【生信小课堂】如何直观的理解RNA-seq分析流程

【生信小课堂】如何直观的理解RNA-seq分析流程

【小撰同学的惯例废话前言】上个星期在我的入门群里(原本是为了让考研群的小伙伴进复试不要一问三不知才开的群,不是孟孟的三千人交流大群)布置了一个任务,跑通RNA-seq的流程,结果跑通的人寥寥无几,并且即使跑通也不理解究竟发生了什么事,仅仅就是傻乎乎的做完。这让我开始反思RNA-seq入门究竟要怎么讲才能最简单的理解,不像我之前开的直播课一样讲得有些晦涩难懂。故而前天我就在群里用最简单的灵魂作图来告诉群里的小伙伴,RNA-seq究竟发生了什么,接着就有你现在看到的这篇废话连篇又很可能有些表述很不科学的文章。

本文适合刚入坑的小伙伴观看,并且在跑流程的时候对着看,可以对数据在干嘛有个粗略的理解,而后再补充到位,而不会说跑完整个流程也不知道在干嘛

要理解整个流程,个人觉得可以按数据的四个流程来拆分:高通量测序准备工作上游分析下游分析

【什么是高通量转录组测序】

所谓高通量测序技术是什么?

顾名思义,就是通量很高(对比sanger测序)的测序,一次性可以获得海量的数据,所以叫高通量测序。

转录组是什么?

转录组,一般指的就是某一时空条件下细胞所产生的所有转录产物,说人话就是,你的样品经过了某种处理,然后拿去提了总RNA,这个总RNA就是一个转录组

理解高通量转录组测序的关键在哪?

首先是建库,我们建的文库用的是什么,rna吗?不是,那用的是什么?

cDNA,即rna拿去逆转录的产物,为什么要用DNA而不是RNA?

除了单链RNA不稳定外,还有一小部分原因是DNA的建库流程已经确定了,只要把RNA变成DNA后面流程完全一样,可以偷个懒,不过为了节约时间可以在一二链合成的时候直接加好接头,后面就连接头都不用加了,缩短建库的时间,一天可以轻松完成建库

其次就是什么是桥式pcr

上面就是桥式pcr的流程,简而言之就是序列接头(adapter)一端被固定,然后另一端跟反应槽里的互补序列配对,呈现桥状,然后再进行pcr,故而称桥式pcr。

经过n轮桥式pcr之后,一个序列可以扩增到一个叹为观止的水平,故而通量就非常高了~

最后是测序信号是怎么得到的?

荧光基团,在所有的碱基上我都接了荧光基团

想更直观理解的这里放个illumina的官方介绍

【准备工作】

首先,我们拿到的原始序列文件就是fastq,那么怎么去理解fastq文件呢?

ST-E00126:128:HJFLHCCXX:2:1101:7405:1133
TTGCAAAAAATTTCTCTCATTCTGTAGGTTGCCTGTTCACTCTGATGATAGTTTGTTTTGG
+
FFKKKFKKFKF<KK<F,AFKKKKK7FFK77<FKK,<F7K,,7AF<FF7FKK7AA,7<FA,,

上面这里是一个fastq文件的格式,每一行代表什么呢?

第一行就是测序的坐标信息,即告诉你这条reads的名字是什么

第二列就是我们测到的序列

第三列就一个加号,没卵用

第四列质量信息,对应着上面各个碱基,测得有多臭,具体多臭下面说怎么直观的看

需要了解的就这么多,如果要仔细了解,看下面这个帖子

怎么对测序质量这些东西进行直观化

有个东西叫做fastqc的软件,可以对fq文件进行质检,具体怎么看呢?看这个贴子

在明白了自己的测序数据有多臭之后,我们就要将数据中低质量的部分全部剔除掉,剔除的软件有很多,类似Trimmomatic,fastp,cutadapter

一般给定的标准就是清除存在的所有接头序列,过滤掉q小于20的碱基,去除N碱基大于5%的序列,去除A与T或者C与G含量相差10%的序列,去除切除碱基后过短的序列,这个标准一般通用,具体可以根据自己数据去筛选。

想详细了解的看这个贴子

准备工作完成之后,我们就得到了一份高质量的原始数据(clean data),从而正式进入分析工作

【上游分析】

无论是以前的bowtie2+samtools+cufflinks+deseq2,还是现在转录组的当红炸子鸡流程hisat2+stringtie+ballgown,其本质的工作流程其实是一样的,只不过使用的算法不同而已。

第一步叫做回帖,这一步是干嘛的呢?

首先,我们的fastq文件存储的数据是一个零散的状态,那要怎么样把它恢复到打断前的状态?

这里我们就需要一个模板,按照模板,我们把序列排序,大概就长这个样子

这里的ref就是模板,即参考基因组,而我们的fastq文件本质就是一条一条的小序列,在模板的指引下,我们得到了他们原本在基因组上应该在位置,这一步就是回帖的含义。

也即是bowtie2跟hisat2所干的事。

而关于回帖的细节,可以看这两篇

回帖完之后,我们的回帖信息会被输入到一个文本文件:SAM文件(二进制位bam文件)

sam文件有个头文件,即你看到这张图前面那样,存储着染色体的信息,还有你之前比对的指令,但这些不是我们需要了解的重点,我们需要来看看下面存储着什么?

第一列是什么?就刚刚fastq文件的第一列,就是这条reads的名字

第二列是什么?flag?太复杂了,不记了

第三列染色体

第四列染色体的起始位置

第五列回帖的可信度,即回帖质量

第六列?第七列?看不懂,不管了

后面还有回帖上的序列

总的来看,所谓的sam/bam文件就是记录回帖的序列是什么,回帖上多少,回帖的质量行不行,回帖到什么位置

而后就是用cufflinks或者stringtie结合注释文件gff/gtf,将转录本构建出来。

那么gff/gtf是什么?简而言之,gff就是记录了这个物种在哪个位置有功能,是gene还是调控因子

而cufflinks要做的事情就是将bam文件的比对信息跟gff的信息结合起来,拼出一条转录本

bowtie2做的事情是

而cufflinks做的事情则是这个

【下游分析】

当使用cufflinks构建得到raw count之后,我们就想比较不同处理间的差异在哪,那么这个时候我们可以直接比较吗

当看到我这么问的时候,肯定就是说不可以

那么,为什么不可以?

拿孟孟之前举的例子

问题1: 比如我有gene3,有1000条测序reads,gene4有2000条测序reads,那么我能否说gene4就一定比gene3的表达量高?

问题2: 比如我有gene1,有1000条测序reads,我的另一个处理条件下gene2有2000条测序reads,我能否就说geneA在处理条件下表达量降低了?

图1 ( Manuel Garber et al., Nature Methods, 2011 )

很明显,第一个问题,如果两个基因的长度不一致,那是无法直接比较的;而第二个问题,我们就需要考虑如何矫正了,而这个矫正值就是所谓的RPKM/FPKM/TPM,关于这些是什么

请看这个贴子

当我们将所有的东西放同一个标准下,就可以进行比较了,而比较的时候,即肯定存在两个组才能进行比较,也就是我们的control跟treatment

以control为标准,比较treatment,我们就知道了差异究竟在哪些基因,即所谓的差异表达

现实计算肯定要复杂的多,但作为粗略理解,这样会比较容易理解

而当我们得到一堆差异基因之后,就通常要做所谓的富集分析,常见的有GO,KEGG。

以常见基于ORA算法的富集举例

本质其实就是一个超几何分布,常见的就是Fisher extract test

算出一个p值即可,然后自选标准,大于多少认为是显著的,认为某某通路上存在差异表达。

想更具体了解,可以看看这个视频

本期内容就到这里,还望各路大神轻喷,同时欢迎各位大神指点一下哪里可以写得很通俗而不失谨慎,方便新入门的小伙伴更好的理解整个分析流程~

关于群号:孟孟的3k群需要在腾讯课堂(推荐)或者知乎live上购买课程,而后在课程互动里面找到群号跟加群的验证(知乎live),或者进入购课群后找群主要入群验证(腾讯课堂)

我的小白互勉群群号:933196696,只推荐没入门或者刚入门的同学加群,因为面向新手,我创建本意也只是为了让考研群的小伙伴进复试不会被一问三不知(很多人考了生信却没有真正实践),大神其实也欢迎加群帮忙解答问题,但您的问题可能我们解答不了,毕竟太菜了

编辑于 2020-05-12 10:52