一篇文章讲清RNA-seq原理与应用 | 生信笔记07

一篇文章讲清RNA-seq原理与应用 | 生信笔记07

引用

CSDN

1.

https://blog.csdn.net/dxs18459111694/article/details/139136824

RNA-seq（RNA测序）技术是现代生物医学研究中的一项重要工具，广泛应用于基因表达分析、转录组研究等领域。本文将系统地介绍RNA-seq数据分析的完整流程，从数据获取到最终的功能注释和整合分析，帮助读者全面了解这一复杂而重要的生物信息学过程。

前言

接下来我们要介绍的是RNA-seq数据的处理分析流程，根据RNA-seq测序技术的不同，可以分为三种：

1. short-read
2. long-read
3. direct RNA-seq

而我们一般的RNA-seq测序数据分析流程算法，基本上都是基于short-read（短读长）技术所产生的数据文件。目前，我们可以从Short Read Archive(SRA)数据库获取的RNA-seq数据中，有超过95%的数据是由Illumina公司的short read测序技术所产生的。

其分析过程可以用下面的路线图表示：

该路线图大致分为三个部分：

1. 数据获取：

• 包括实验设计、测序设计以及数据下机后的raw reads数据的质控

2. 数据分析

• 在获取到干净的数据之后，可以进行reads的比对，然后进行基因表达的量化、差异表达分析、功能富集分析等

3. 高级分析

• 包括数据的可视化，其他小分子RNA分析、融合分析以及与其他类型的数据进行整合分析等

而我们分析的起始点，是从原始数据开始的，也就是获取raw reads数据。通常这种高通量测序数据会保存为FASTQ格式的文件。FASTQ格式是一种以ASCII码字符的形式保存生物序列及其对应的每个碱基的质量的文本文件。

FASTQ文件中每条序列（通常是一条read）是由4行组成，其中：

• 第一行以@字符开头，之后的字符为序列的标识符和描述信息
• 第二行为具体的序列
• 第三行以+符号开头，之后可以可选地加上与第一行一样的序列标识或描述信息
• 第四行为碱基质量分数（Phred），其字符数量与第二行相等，每个字符表示对应碱基的质量得分，例如

@SEQ_ID
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65

其中，碱基质量值的编码方式为

1. 先将碱基错误率P进行负对数转换，得到Q值
   Q = − 10 log10(P)
2. 然后将Q值加上33或64得到的值所对应的ASCII码即为碱基质量分数

例如，错误率P = 0.01，则Q = 20，如果是Phred33则对应的质量为字符5（53），如果是Phred64则对应的字符为T(84)

分析流程

1. 数据获取

一般情况下，如果自己有送样检测数据的话，测序公司会提供原始的FASTQ格式的数据。如果我们要使用别人文章中发表的公开数据，还需要从数据库中下载对应的数据

例如，我们从SRA数据中下载的原始测序文件是sra格式，我们需要先使用工具将其转换为FASTQ格式

2. 质量控制

主要在三个地方需要对数据的质量进行监控

• 获取原始数据之后
• 比对完之后
• 表达定量之后

2.1 raw read

对raw reads数据进行质量控制，需要分析序列的质量、GC含量、是否存在接头、短重复序列的分布、测序错误以及PCR重复和污染

质控软件：

• FastQC：用于分析Illumina测序平台的数据
• NGSQC：可应用于所有平台

一般来说，reads的质量会朝着3'端递减，如果碱基的质量太低，我们需要删除它以提高比对率

FASTX-Toolkit和Trimmomatic两个软件可以用于切除低质量的碱基和接头序列

2.2 比对后

reads通常需要比对到一个参考基因组或转录组，而比对的质量是评估测序准确率和是否存在DNA污染的一个重要指标

比对质量通常为比对到的reads数占总reads数的比例。例如，比对到人类参考基因组的比对质量通常需要在70-90%，且有大量的reads映射到一个相同的区间内。如果是比对到转录本上，由于可变剪切的影响，可以适当放宽比对质量

在外显子和比对方向上的read覆盖率的均一性，也是评估质量的重要指标。如果 reads 主要聚集在转录本的3'端，可能表明原始样本的RNA质量较低

比对上的reads的GC含量，可能揭示了PCR的错误率

主要软件有：Picard、RSeQC和Qualimap

2.3 定量后

在计算完表达的量化值之后，可以计算GC含量和基因长度的误差，在必要时可以使用标准化方法来进行校正

如果参考转录组注释得很好，则可以分析样本的生物构成，来评估RNA纯化步骤的质量。例如，rRNA和small RNA不能出现在polyA longRNA的制备中

NOISeq和EDASeq等R包可以使用图形来展示count数据的质量控制

2.4 可重复性

上面的质量控制都只是针对单个样本的，此外，不同样本之间的可重复性评估，对于评价整个数据集的质量也是至关重要的

技术重复样本的可重现性一般很高（spearman R2 > 0.9），但是生物学重复样本之间并没有明确的标准，取决于实验系统的异质性。如果不同实验系统之间存在差异基因，则同一条件下的生物学重复在主成分分析（PCA）中会被聚类在一起。

3. 序列比对

在对样本的raw reads进行质控之后，就可以进行序列比对了，序列比对主要有三种策略，如下图

如果有参考序列，根据参考序列的不同，可以分为

• 比对到基因组：使用间隔比对算法，如TopHat、STAR等，然后根据是否提供了注释文件（GFF格式文件，包含转录本位置信息），又可以分为转录本识别和转录本发现并进行定量分析
• 比对到转录组：使用非间隔比对算法，如Bowtie等，然后使用RSEM或Kallisto方法识别转录本并计算定量信息
  如果没有参考序列，则需要先把序列组装成转录本，再将reads比对到组装后的参考转录本上，然后使用HTseq-count等算法对转录本进行定量

3.1 转录本发现

使用Illumina技术检测的short reads来发现新的转录本是RNA-seq分析中的一个挑战。通常来说，短reads很少会跨越多个剪切位点，这就很难直接推断出一个转录本的整体长度。

此外，转录的起始和终止位置也比较难识别，一些像GRIT的工具，通过合并5'端的信息可以提高异构体识别的准确性。其他如Cufflinks、iReckon、SLIDE和StringTie等方法，通过结合现有的注释信息，作为一个可能的异构体列表

一些寻找基因的工具，如Augustus，结合RNA-seq数据，可以更好的注释蛋白编码转录本，但是对非编码转录本的性能更差。

3.2 De novo 转录本重构

在没有转录本或转录本不全的情况下，可以对reads进行组装来重构一份转录本。可选的方法很多，如SOAPdenovoTrans、Oases、Trans-ABySS或Trinity

通常来说，使用双端链特异性测序和long reads测序包含更多的信息，会有更好的效果

虽然，对于低表达的转录本进行组装的可靠性较低，但是reads太多也会导致潜在的组装错误和较长的时间消耗等问题。因此，在深度测序的样本中，可以适当减少reads的数量

对于多样本的比较，可以将所有样本作为一个输入来构建参考转录本，然后分别对每个样本的reads进行比对

无论是使用参考序列还是从头开始组装，使用短reads的Illumina技术来完全重构转录组仍然是一个具有挑战性的问题

4. 转录组定量

RNA-seq最广泛的应用就是用来评估基因和转录本的表达，这一应用主要是基于比对到转录组区间内的reads的数量

最简单的方法是，使用HTSeq-count或featureCounts计算区间内的reads数来量化基因的表达。这种基因水平的（不是转录本水平）的量化方法使用的是GTF文件，这种文件包含外显子和基因在基因组上的坐标。

但一般不能直接使用read count来比较基因的表达水平，因为该值会受到转录本长度、reads总数以及测序偏差等因素的影响。所以需要先进行标准化，标准化方法有

1. RPKM/FPKM : 每百万reads每一千碱基对中包含的reads数

该方法先计算测序深度系数，即总reads数除以 一百万，然后计算基因或转录本的长度（单位为kb），标准化顺序为先消除测序深度的影响，再消除长度的影响：

RPKM(x)=Reads per transcripttotal reads106×transcript length1000=Reads per transcriptmillion reads×transcript length (kb)=109×CxR×Lx

其中

• x表示一个基因或转录本，或基因组上一段特定的区域
• Cx表示比对到x外显子区域的reads数；
• R表示当前样本中包含的全部reads数
• Lx表示x外显子区域包含的碱基数（长度，bp）

FPKM与RPKM的计算公式一样，只是RPKM用于单端测序，FPKM用于双端测序

1. TPM : 其与RPKM最大的区别是，标准化顺序为先消除基因长度的影响，再消除测序深度的影响

首先，将reads count除以基因或转录本的长度（kb）得到RPK(reads per kilobase)，然后将样本中所有的RPK加起来除以106，得到标准化系数，最后使用RPK除以标注化系数

TPM(x)=Cx/Lx×106∑i=1NCi/Li

其中

• x表示一个基因或转录本，或基因组上一段特定的区域
• Cx表示比对到x外显子区域的reads数
• Lx表示x外显子区域包含的碱基数（kp）
• N表示基因或转录本总数

这样，每个样本的TPM总和是一样的，便于比较样本间的差异

目前，也有许多复杂的算法通过解决相关转录本共享reads的问题来评估转录本水平的表达，例如，Cufflinks使用TopHat的比对结果，应用期望最大化算法来评估转录本的丰度。这一方法考虑到长度不同的基因的reads分布并不均匀等因素的影响。

还有其他算法也可以量化转录组的表达，例如RSEM、eXpress、Sailfish和kallisto等。这些方法允许转录本之间存在多比对的reads，并输出经测序偏差校正的样本内归一化值。

5. 差异表达分析

差异表达分析是对样本间基因的表达值进行比较，虽然RPKM、FPKM和TPM标准化方法消除了测序深度和基因或转录本的长度因素的影响，但这些方法依赖于总的或有效的reads数，当样本的具有异质性转录本分布或当高表达或差异表达的特征扭曲了count分布时，表现欠佳

而像TMM、DESeq、PoissonSeq和UpperQuartile等方法会忽略高变异或高表达的特征。

干扰样本内比较的其他因素包括不同样本的转录本长度变化、转录本覆盖位置的偏差、平均片段大小以及基因的GC含量等

NOISeq这个R包提供了多种绘图，来识别RNA-seq数据中的误差来源，并应用相应的方法来标准化这些误差

除了这些样本内特异的标准化方法，还需要解决数据集之间的批次效应（不同实验条件下产生的数据之间存在的差异），批次矫正方法有COMBAT和ARSyN等，虽然这些方法是针对芯片数据设计的，但是在RNA-seq数据中也有很好的效果

计算差异表达的方法有很多，有些方法，如edgeR将原始的read counts作为输入，并在统计模型中加入了标准化，另一些方法，需要先对数据进行标准化，如DESeq2使用的是负二项分布作为参考分布，并提供了自己的标准化方法。

baySeq和EBSeq是贝叶斯方法，还有一些基于线性模型的方法。最后，一些非参数方法，如NOISeq和SAMseq对于小样本量的研究，负二项分布会存在噪音污染，这种情况下，一些简单点方法，如基于Poisson分布的DEGseq，或者基于经验分布的NOISeq可能会更好些。

但是需要强调的是，在没有足够生物学重复的情况下，无法进行总体的推断，因此任何p值计算都是无效的。

许多独立的研究都已经证实，选择不同的方法会对结果有一定的影响，而且没有哪一种方法能够适用于所有的数据，所以，推荐在分析的时候使用多个软件进行相互验证。

6. 可变剪切分析

可变剪接(Alternative Splicing) 是指转录形成的前体RNA通过去除内含子、连接外显子而形成成熟RNA的过程，从而实现一个基因同时编码多种蛋白质，实现生物功能多样性

在不同组织或者发育的不同阶段，可变剪接不是一成不变的，在特定的组织或条件下，通过连接不同的外显子，会产生特定的剪接异构体（isoform）。有大量的研究发现，可变剪接的变化与癌症等多种疾病相关，所以研究可变剪接在不同组织中的作用是非常有意义的。

转录本水平的差异表达分析可以潜在地检测同一基因的转录异构体表达的变化，已经有一些算法应用于RNA-seq数据的中进行可变剪切分析

这些方法主要分为两大类：

1. 异构体表达估计与差异表达相结合，来揭示总基因表达中每种异构体的比例变化

例如，BASIS方法使用分层贝叶斯模型来直接推断转录异构体的差异表达；CuffDiff2方法先评估异构体的表达，然后比较它们之间的差异；rSeqDiff方法使用分层似然率检验同时检测无剪接变化的差异基因表达和差异异构体表达。所有这些方法通常都受限于短读长测序的内在局限性，无法在异构体水平上进行准确识别

1. 一种所谓的exon-based的方法，它跳过了对异构体表达的估计，通过比较样本之间基因外显子和连接点上的reads分布来检测可变剪接的信号

其基本假设为：可以在外显子及其连接点的信号中追踪异构体表达的差异。DEXseq和DSGSeq采用类似的思路，通过检测基因的外显子（和连接点）上read counts的差异显著性来识别不同的异构体。rMATS是通过比较用连接点的reads定义的外显子inclusion levels表达水平的差异

7. 融合分析

基因融合是指两个基因的全部或一部分的序列相互融合为一个新的基因的过程。其有可能是染色体易位、中间缺失或染色体倒置所导致的，可在DNA或RNA层面上表达。

融合基因通过基因失调、融合产生嵌合体蛋白这两种机制引发癌症的发生。

目前，RNA-seq融合算法100多种，有人对常用的15中融合检测算法进行了比较

没有哪一个算法具有明显的优势，整体来看，SOAPfuse可能会好一些，FusionCatcher和JAFFA其次。

8. 功能注释

标准的转录组分析的最后一步，是使用差异表达基因来进行功能或通路的注释。最常用的两类方法是：

• 基于超几何分布的过表达富集分析
• GSEA富集分析

一些工具，如GOseq考虑了基因长度等因素对差异表达结果的影响，并使用超几何分布进行富集分析，GSVA和SeqGSEA使用类似GSEA的方法进行功能富集

功能富集需要预先定义的基因集合或通路，包括GO、KEGG、Reactome等数据库。

通过在蛋白质数据库（例如SwissProt）和包含保守蛋白质结构域（例如Pfam和InterPro）的数据库中搜索相似序列，使用直系同源分析对蛋白质编码的转录本进行功能注释。而Rfam数据库包含许多特征明确的RNA家族，例如rRNA或tRNA，而mirBase和Miranda是专门研究miRNA的

9. 整合分析

将RNA-seq数据与其他类型的全基因组数据进行整合分析，使我们能够将基因表达的调控与分子生理学和功能基因组学的特定方面联系起来。

1. DNA测序
   将RNA测序和DNA测序相结合，可以进行SNP、RNA编辑和表达数量性状基因座（eQTL）比对等分析。

2. DNA甲基化
   将DNA甲基化和RNA-seq整合，可以分析差异表达基因和甲基化模式之间的相关性。使用的算法包括：广义线性模型、logistic回归模型和经验贝叶斯模型等

3. 染色体特征
   通过整合RNA-seq和Chip-seq数据，可以降低Chip-seq分析的假阳性，并展示TF对其靶基因的激活或抑制作用。

4. MicroRNA整合
   RNA-seq和miRNA-seq有可能揭示miRNAs对转录稳态水平的调节作用

5. 蛋白质组和甲基化组
   RNA-seq与蛋白质组学的整合是有争议的，因为这两种测量结果通常显示出很低的相关性（~0.4）。尽管如此，蛋白质组学和RNA-seq的配对分析可用于识别新的异构体

转录组与代谢组数据的整合已被用于识别在基因表达和代谢物水平上受调控的通路，并且可以使用工具来可视化通路上下文的结果，如MassTRIX、Paintomics、VANTED v2和SteinerNet