知识中心 - 北京概普生物科技有限公司(GapTech)

浅谈转录组Read在基因组分布和饱和度分析

生信干货生信人 ·2018年8月11日 17:29

1 介绍

一旦reads 被mapped 到参考基因组，他们的位置就可以匹配到基因组上的注释信息。这样，我们可以通过计算基因，转录本或者外显子的reads 数来定量基因的表达。使用被比对到基因组的reads对测序饱和度，不同基因组特征的read 分布，转录本的均匀覆盖度进行评估。

本章前半部分介绍基于注释的定量，并且提供一些软件检验它们。下半部分讨论基因表达的定量，这是大多数RNA-seq 研究的组成部分。原则上，计算定位read 的数量提供了一个直接估计转录本丰度的方法，但在实践中需要考虑一些并发状况。真核基因通过选择性剪接和启动子的使用，通常产生一个基因的多种剪接异构体。但是，用短reads在转录本水平定量是没有意义的，因为转录本不同异构体之间往往有共同或重叠的外显子。而且，由于比对问题和文库偏差会导致转录本覆盖度不均匀。因为这些并发状况，表达经常在基因水平或者外显子水平进行分析。但是，由于不同条件下采用不同的可变剪接体，而不同的可变剪接体长度也不同，因此在进行差异表达分析时，基于基因水平计数通常不是最佳的。这个问题在第八章差异表达分析的部分有详细论述。

2 基于注释的定量质控

正如在第3章中讨论的，生产RNA-seq数据的实验室步骤还不完美，但幸运的是，很多reads质量问题例如低测序质量的碱基和核苷酸组成上的偏差已经可以在原始数据水平检测出来。但是一些重要的质量，仅仅在reads被mapped 到参考基因组以及它们的位置匹配到注释信息的时候，才能被测量出来。这些质量包括：

1）测序饱和度。基因表达的丰度分析，可变剪切分析和转录本重建分析的可靠性都依赖测序深度。因为测序花费高，所以检查数据的丰度是很重要的，也即是说，要发现新的基因或者剪切点是否要增加测序量。理想情况下，会预先确定合适的深度，但是这需要来自同一物种和组织的数据，因为丰度依赖转录组的复杂性。

2）不同基因组特征的read分布。这个可以在多个水平来做，例如,reads可以外显子，内含子和基因间区为单位进行计数。外显子可以分布在编码区之间，3′UTR和5′UTR。如果很高比例的reads映射到内含子和基因间区，应该有新的异构体和基因，但这也可能是污染基因

组DNA的标志。reads可能更多地分布在一些生物分子如蛋白基因、假基因，核糖体rRNA、miRNA等。rRNA含量尤其重要，因为实验室移除rRNA的步骤可能会导致样本不可靠或者样本不一致。如果reads大量地映射到rRNA，你可以移除它们，例如，用Bowtie2(如第4章所描述的)将它们映射到rRNA序列，并保留未比对上的reads。

3）转录本的覆盖均匀度。不同的实验室操作方法会引入不同的偏差。例如包含有polyA帽子捕获步骤的会导致reads主要来自于转录本3′端。3′端的偏差在不同的样本里表现多种多样，因此，有必要重新估计它的丰度。

2.1 基于注释的质量控制工具

许多质控RNA-seq数据的工具是可以用的，RSeQC【1】, RNA-SeQC【2】, Qualimap【3】,和Picard’s CollectRNASeqMetrics tool【4】。它们给出了许多overlap 的质量统计，并且各有优点。他们都提供命令行接口，RNAseQC 和Qualimap 有它们自己的GUI，并且RseQC在Chipster软件中可用。注释信息通常是在GTF[5]或BED[6]文件中给出的，这些文件需要与BAM文件有相同的染色体命名。

RNA-SeQC用Java实现，它可以从GTF获取注释。它还需要一个带有索引（.fai）的参考和序列字典文件(. dict)的引用FASTA文件。RNA-SeQC提供了一个特别详细的覆盖深度报告，并且可以比较不同的样本。覆盖率统计报告包括平均覆盖率,转录本两端区域的覆盖度,3′,5′的偏差,和gap的数量,累计长度和占的百分比。所有的值都分别对低、中、高表达基因计算。覆盖率值对三个级别的GC内容进行统计。另外，除了覆盖均匀分布图,RNA-SeQC也给出离3′端距离与覆盖深度分布图。输出结果由HTML报告和TAB分隔的文本文件组成。

Qualimap是一个Java 程序，它在内部使用R和某些Bioconduct 包。它可以从GTF/BED格式获取注释信息，而且它也需要一个独立的生物分子文件。Qualimap可以对丰度和遗传丰度绘制很好的图。丰度图可以展示在不同的测序深度中检测到的不同特征的数量，而且它也很方便地指出每增加一百万测序深度有多少新的特征会被检测出来。分子类型分布图显示了在蛋白质编码基因、假基因、rRNA、miRNAs 等上面的reads分布，以及基因组中这些特征上有多大的百分比。

RseQC由几个Python程序组成，并且从BED格式的基因组文件中获取注释信息。注意，使用R来绘制图形。RseQC有几个其他程序没有的功能:(a)计算不同基因组特征的分布时，同时也包含了在转录本的上游和下游的几个bins。(b)更重要的是，除了基因外，它还可以对剪切位点计算丰度状态，并且(c)可以注释剪切位点为已知、新的和部分新的。BED文件有3个必须的列，9个可选[6]。RseQC需要12列的BED文件，因为每个基因的外显子信息包含在最后的三列中(blockCount、blocksizes和blockstart)。你可以用UCSC 的table Brower 获得不同组织的BED文件，在分组菜单中选择 ¨ Genes and gene predictions.¨ ,从菜单中选择一个基因集。例如，Refseq 基因或者Ensemble 基因，并且输出文件设置为BED格式，来自UCSC的BED文件中的染色体有前¨ chr¨,用Ensemble则没有前缀，你可以使用下面的命令删除¨ chr¨ :

Sed ‘s/ ∧ chr// ‘hg19 Ensembl chr.bed > hg19 Ensembl.bed

下面RseQc例子中使用的双端比对文件，accepted hits.bam，选自第四章.

read distribution.py用来计算不同基因组特征型的reads 分布。

python read distribution.py −r hg19< E nsembl.bed −i accepted h its.bam

下面的结果表中给出了总的reads 数（excluding nonprimary hits）和标签（separate splice fragments of a read）。所有的标签指明了有多少标签可以明确地分配给下面列出的是个种类。

使用geneBody coverage.py来产生覆盖率图，用来检查沿着转录本的覆盖是不是均匀的。参数−o用来给出结果文件的前缀。

python geneBody coverage.py −r hg19 Ensembl.bed −i accepted h its.bam −o file

使用RPKM saturation.py来估计在当前深度下基因表达丰度估计的精确性。从reads 中随机取样，以RPKM 为单位对每一个子集计算丰度。检查他们是否稳定。这个在四个水平下分别来做。

pytho RPKM saturation.py −r hg19 Ensembl.bed −i accepted hits.bam −o file

使用junction annotation.py来将剪切点分为新的，部分新的（有一个是新的）或者注释到基因模型（两个都在参考基因模型上），用一个饼图指出结果。

python junction annotation.py −r hg19 Ensembl.bed −i accepted hits.bam −o file

可以使用junction saturation.py来检查剪切位点的测序丰度。从reads 中随机取样，检测每一个子集，与参考基因组注释比较，结果是分别指出新的剪切点和已知的剪切点。

python junction saturation.py −r hg19 Ensembl.bed −i accepted hits.bam −o file

Chipster中基于注释的质量控制选取BAM文件，包含注释的BED文件和使用RseQC进行质量控制或者RNA-seq质量定量。按照参数表，确定文件被正确分配。

3 总结

借助于参考注释信息，对比对后reads在基因组上的位置分布特征进行分析，方便我们研究一些重要的质量情况，例如测序深度丰度，转录本的覆盖均匀度和不同基因组特征的read分布。一些工具可以用来基于注释的质量控制，并且它们各有优势。当reads被mapped到参考基因组之后，我们可以通过计算基因，转录本或者外显子的reads数来定量基因表达。定量和差异表达分析是天然联系在一起的，最好的方法仍在讨论中。reads可以使用工具如HTSeq来计算，但是对于不同剪接体转换的基因的差异表达分析，基因水平的计数不是最优的（因为长转录本有更多的计数）。

如有问题请扫描二维码，去WeGAP论坛提问，WeGAP各板块招募版主，欢迎留言。

欢迎关注生信人

TCGA | 小工具 | 数据库 |组装| 注释 | 基因家族 | Pvalue

基因预测 |bestorf | sci | NAR | 在线工具 | 生存分析 | 热图

舞台|基因组 | 黄金测序 | 套路 | 杂谈组装 | 进化 | 测序简史