RNA-seq测序方法及数据分析

RNA高通量测序（RNA-sequencing，缩写为RNA-seq）是目前高通量测序技术中被用得最广的一种技术。

RNA-seq可以帮助我们了解：各种比较条件下，所有基因的表达情况的差异。

RNA-seq可以检测的差异有：正常组织和肿瘤组织的之间的差异，药物治疗前后基因表达的差异，发育过程中不同的发育阶段不同的组织之间的基因表达差异，等等。

在所有检测的差异类型中，最常见的就是检测所有mRNA的表达量的差异。

同时，还可以检测 RNA 的结构上的差异。例如：mRNA的剪接方式的差异，也就是我们一般说的“可变剪接”，还可以检测“融合基因”，同时还可以检测基因单点突变导致的SNP（Single Nucleotide Polymorphisom)。

RNA-seq测序方法

去除核糖体RNA、建库

在测mRNA的过程当中，首先要解决的问题，就是如何去除核糖体RNA，即rRNA”(Ribosomal RNA)。

在通常抽提到的总RNA中，绝大部分都是核糖体RNA（rRNA）。以人类的细胞或组织为例，一般抽提到的总RNA当中，95%都是核糖体RNA。剩下的2%到3%是mRNA。还有2%到3%是Long non-coding RNA、或者tRNA、microRNA这些RNA。也就是说，mRNA只占了所有RNA中的一小部分。

如果把所有的RNA都拿来测序的话，测到的绝大部分的序列数据都是核糖体RNA。而且这当中rRNA的比例会高达95%左右，但是，核糖体RNA在整个人类当中都是非常保守的，而且在人的各个组织、器官当中也是极度稳定的。也就是说，测rRNA，它得到的数据，并不能为实验者提供什么有用的信息，而mRNA才是RNA当中信息含量最丰富的那个部分。

我们一般的RNA-seq要测的，也是mRNA的各种变化，所以，在实验过程当中，我们一般要把核糖体RNA先去掉。然后再进行建库测序。

去除核糖体RNA，并进行建库的方法有许多种。目前应用最广泛的是illumina公司的TruseqRNA建库方法。

上图是mRNA测序的建库过程图。

首先，利用高等生物的mRNA都有Poly(A)尾巴这个特点，用带有Poly(T)探针的磁珠与总RNA进行杂交。然后，Poly(T)探针就和带Poly(A)尾巴的mRNA结合在一起，接下来就回收磁珠，然后就把这些带Poly(A)的mRNA从磁珠上洗脱下来。

然后，再把这些洗脱下来的mRNA用镁离子溶液进行处理。镁离子溶液会把mRNA打断。被打断的这些mRNA片段，再用随机引物进行逆转录。

逆转录成（第一链）cDNA后，再合成出第二链（cDNA)。这样就成为双链的cDNA。

接下来，再在双链的cDNA的两端加“A”碱基，并连上“Y”型的接头。经过PCR扩增，成为标准的测序文库，然后，这个标准的测序文库就可以拿到HiSeq测序仪上进行测序了。

样本质量要求

上述建库方法对RNA的完整度有较高的要求。只有在mRNA大部分是完整的状态下，才能得到比较好的效果。

这是因为带Poly(T)的磁珠，它所吸附的是Poly(A)序列。如果mRNA发生了降解，即mRNA断掉了，那么磁珠所吸附下来的片段都是靠近3'端的断片，而靠近5'端的断片是吸附不下来的，会在富集过程中被洗脱掉。这样在接下来的数据分析当中，就会发生一定的数据偏差。

为了保证能够测到尽可能完整的mRNA序列，Illumina公司是这样建议的：先对总RNA进行一次质量检测，一般是用Agilent公司出品的Bioanalyzer 2100毛细管电泳仪，对总RNA样本进行一次电泳质检。Bioanalyzer会根据18S和28S这两个核糖体RNA的电泳峰是否高、是否尖，来判断RNA的质量。并且会自动打分。如下：

这两个峰越高、越尖，也就说明RNA的降解就越少，完整度就越高，打分也会越高。反之，打分就会低。这个分值叫“RIN”值（“RNA Integrity Number”），即RNA的完整度评分值。RIN值最高是10分，最低是0分。

Illumina公司推荐用RIN值在8.0以上的RNA进行建库和测序。测序完成之后，就可以进行数据分析了。

RNA-seq数据分析

判断测序的质量

分析的第一步，一般是先把测到的RNA片段，先mapping(比对)到基因组上。在比对完后，可以先看一下，有多少RNA片段是在靠近基因的5'端位置，又有多少片段在是靠近基因的3'端位置。

上图就是把所有的基因，都按其外显子的长度拉直，然后归一化到“0 - 100”的长度。看比对上的片段有多少落在0到100这一个轴的哪个位置上。

比对的结果可以让我们看到前面Poly(T)磁珠在抓mRNA的时侯，捕获下来的这些mRNA是不是完整的。

如果捕获下来的这些mRNA大部分是完整的话，那么这个图形靠近5'端的曲线就会显得比较饱满。它的高度会和3'端的高度差不多。反之，如果这条曲线的3'端很高，而5'端比较低，我们就可以初步判断，该RNA有一定程度的降解。

因此，我们可以推断在捕获过程当中，有相当一部分mRNA的5'端片段因为与3'端片段的Poly(A)片段的尾巴断开了，所以没有被捕获下来。所以该RNA是有一定程度降解的。

在知道了测序的质量之后，接下要关注的就是不同样本之间、各个基因的mRNA的表达量的差异。

数据标准化（RPKM、FKPM、TPM、CPM）

统计比对到基因上的reads数即为counts，也就是测序原始表达矩阵，rawdata。但由于两大原因（不同样本的测序深度，不同基因的长度），直接用counts比较没有意义，因此需要进行标准化。RPKM、FKPM、TPM是实现消除上述二者影响的方法，三者差异如下。

RPKM：

Reads Per Kilobase of exon model per Million mapped reads，代表每一百万条可以比对到基因组上的Read当中，有几条是可以比对到某个特定基因的，然后这数值再除以该基因的外显子的长度，得到的这样一个最终的比值。即某一基因的counts先除以测序深度（总reads数），再除以基因长度。公式如下：

公式的理解：

①去除测序深度的影响：比对到某个基因的外显子上的Read数，除以这次所测到的、全部可以比对到基因组上的Read数。由于总reads太大了，直接除以这个数字就会使得标准化出来的Read数出现太多的小数，所以为了美观，一般都是除以以百万为单位的总Read数。

②去除基因长度的影响：除以这个外显子的长度，按KB为单位，即1000。目的是修正这个mRNA长度所引起的mRNA的Read数的偏差。这是因为建库过程当中，这个RNA是用镁离子溶液来处理，然后打断（并逆录）成若干个180-200BP左右的小片段，如果一个基因的外显子越长，那么它所产生的mRNA就越长，被打出来的小片段就越多。

注意顺序：是先除以总reads数；再除以基因长度。

FPKM：同RPKM是一样的，只是RPKM用于单末端测序，而FPKM用于双末端测序。

TPM：TPM的计算方法其实也同RPKM很类似，同样的对基因长度和测序深度进行标准化；即counts先除基因长度，再除总reads数。这样每个样本最后的“结果和”都相等，不同样本间差异更清楚。这就意味着TPM数值能体现出比对上某个基因的reads的比例，使得该数值可以直接进行样本间的比较。

事实也证明TPM的标准化方法更有优势，目前都已经推荐进行TPM标准化，不再使用了RPKM、FPKM了。

CPM：Counts per million (CPM) mapped reads，只对测序文库（每个样本总reads数）标准化，而不对长度标准化。这是因为，差异分析往往是同一基因在两组或多组样本量的差异，因此不必在计算单位长度基因的表达量。

RNA表达量差异分析（火山图、聚类分析图、GO分析、KEGG分析）

火山图：

作为一种针对全转录组的分析，我们希望是一次看到一个整体的样本表达差异的情况，而不仅仅是看少数几个基因的表达差异。

科学家做了一种叫“火山图”的一个图形，就像火山喷发的样子，来比较形象地来说明2个样本之间的表达差异，即RNA表达量的对比。火山图就为我们提供了一个形象的、直观的、整体表达差异信息。

横轴表示某个基因的表达是上升了，还是下降了。如果这个基因的表达上调了，那么这个点就往右移动。反之，如果这个基因的表达量下调了，那么这个点就往原点的左移动。

纵轴表示这种变化差异的置信程度，置信程度越高，那么这个点的纵轴位置也越高。

这其中的每个点，就是两个样本当中同一个基因的mRNA表达量的变化。

我们在纵轴上划这样一条水平线，超过这个水平线以上的点，其差异水平的置信程度是很高的，我们就把它标示成红颜色。如果低于这条水平线的点，它的置信程度也相对低一些，我们把它标成蓝颜色。

为什么差异程度是相同的情况下，它们的差异置信程度是不一样的。比如说同样是差了2的5次方，也就是32倍，它的差异置信程度会不一样，有些是蓝点，有些是红点。举例如下：

A基因在甲样本中，被测到了3200条，而在乙样本中被测到了100条；B基因在甲样本中，被测到了320条，而在乙样本中被测到了10条。它们同样是差了31倍，但是因为A基因的样本统计数远大于B基因的样本统计数，也就是说，它们的Reads数有很大的差距。所以，A基因表达差异的置信程度会比B基因表达差异的置信程度要高许多。

聚类分析图：

聚类分析是RNA分析中非常常用的一个手段。它是通过多个样本的全基因表达谱对比，找到它们之间的相似性和相近关系。

上图是一张聚类分析的图，横轴是样本，纵轴是基因。通过聚类分析可以发现：在这个群体中，样本被分成了3个群体。每个群体的内部，都有相似的表达特征。同时，我们还可以看到，基因的表达也是成簇的，大体上分成3个基因群。这3个基因群，各自有着相似的表达量。

聚类分析有很多应用，比如说：我们可以分析疾病的亚型，还可以通过对多个基因在特定疾病当中的表达倾向性来找出可能的、新的、诊断用的Biomark。

GO分析：

GO分析是RNA-seq分析中非常常用的一种分析。GO是Gene Ontology的缩写，Gene Ontology是一个国际化的、基因功能分类体系。这个体系用一整套动态更新的标准词汇和严格定义的概念，来全面地概括任何生物中基因和基因产物的属性。

GO主要描述基因的三个属性：基因参与的生物过程（biological process, BP），基因产物的功能（molecular function, MF），基因产物在细胞器内的空间定位（cellular component, CC）。

差异基因GO富集柱状图：可以直观的反映出在生物过程、细胞组分、和分子功能富集的差异基因的个数分布情况。

有向无环图，是差异基因GO富集分析的图形化展示方式，从上到下所定义的功能范围越来越小、越来越精准。它的分支，表示包含关系。圆圈的颜色越深，表示这个富集关系程度越高。

Pathway分析（KEGG分析）：

通路分析：通路（Pathway）是指在系统水平上完成生物的某一功能的基本单元、或者局部子网络。

KEGG（Kyoto Encyclopaedia of Genes and Genomes，《京都基因和基因组百科全书》）是目前公认的、最权威的基因功能数据库。其中的Pathway(通路）是KEGG的核心内容。目前针对Pathway的分析、注释，大多数是基于KEGG Pathway来做的。

散点图是KEGG富集分析结果的图形化展示方式。

图中，KEGG富集程度通过Rich factor、Qvalue和富集到此通路上的基因个数来衡量。

点的面积越大，则富集的基因数越多。富集的因子越大，则表示富集的程度越大。qValue是校正之后的pValue，越接近0，表示富集程度越显著。

RNA结构变异分析（可变剪接、融合基因、点突变）

结构上的变异，也就是RNA序列的变异。主要是3种：可变剪接、融合基因、点突变（SNP）。

结构分析需要较深的测序深度，一般建议测10G以上的数据量。原因是二代测序目前的测长还不是很长，每一个Read只有大约100到125个Bp左右。如果测序深度不够，那么读到的这些read在整个的mRNA上的分布，是一种比较零碎的一种状态。在这种比较零碎的、不完整的覆盖情况下，要去分析哪里有一个剪接点、断点、SNP，不是很准确。