高通量测序质控及可视化工具包RSeQC

0.65字数 3232阅读 6435

本期介绍一个测序质控的工具包:RSeQC包,它提供了一系列有用的小工具能够评估高通量测序尤其是RNA-seq数据.比如一些基本模块;检查序列质量,核酸组分偏性,PCR偏性,GC含量偏性,还有RNA-seq特异性模块:评估测序饱和度映射读数分布覆盖均匀性链特异性转录水平RNA完整性等。

安装:

RSeQC是依赖于python的,直接使用pip进行安装:

pip install RSeQC
输入数据格式:

RSeQC接受4种文件格式:

  • BED格式: Tab分割,12列的表示基因模型的纯文本文件,例如:
bed12格式
  • SAMBAM格式: 用来存储reads比对结果信息.SAM是可读的纯文本文件,然而BAMSAM的二进制文本,一个压缩的可索引的reads比对文件. 例如:
SAM格式
人hg19染色体信息
  • Fasta文件.
使用方法:

最新版本的RSeQC(2.6.4)包含以下一些模块,每个模块都可以单独调用进行分析:

RSeQC包含的模块

我们一个一个来看:

bam2fq.py:

BAMSAM格式的文件转为fastq格式.(这个感觉一般用不到)


bam2wig.py:

将BAM文件转为wig/bigwig格式的文件.(这个在作图尤其是信号图的时候很有用.)如果需要转为bigwig格式,则需要UCSC的wigToBigWig工具,下载地址:
http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/wigToBigWig


bam_stat.py:

对比对结果文件BAMSAM文件进行统计.其实samtools里也有类似工具.结果如下所示:

bam_stat统计结果

统计结果包括:总比对记录,PCR重复数,Non Primary Hits表示多匹配位点.不匹配的reads数,比对到+链的reads,比对到-链的reads,有剪切位点的reads等.


clipping_profile.py:

这个模块用于评估RNA-seqBAMSAM文件中的有切除核苷酸的reads情况.
  clipped reads有两种,一种是soft-clipped,即reads5'或3'不能比对到参考基因组;另一种是hard-clipped,即reads5'或3'不能比对到参考基因组并且被剪切.
  这个模块会生成.r格式的作图脚本以及.pdf格式的报告文件以及.xls的数据文件.
例如双端测序clipping profile图如下:

Read-1
Read-2

deletion_profile.py:

也就是reads deletion位点的分布.

read deletion distribution

divide_bam.py:

随机分割BAM文件(m个比对结果)为n个文件,每个文件包含m/n个比对结果.


FPKM_count.py:

根据read countgene注释文件(bed12格式)计算每个基因的FPM (fragment per million)或者FPKM(fragment per million mapped reads per kilobase exon).

结果类似下面的:

FPKM/FPM结果

geneBody_coverage.py:

计算RNA-seq reads在基因上的覆盖度.

方法示意图

结果如下:

多个样本的RNA-seq reads在基因上的覆盖度

如果是大于3个的样本,则还会生成热图:

reads在基因上的覆盖度热图
geneBody_coverage2.py:

功能和上面的geneBody_coverage.py一样,但输入的是bigwig格式文件

infer_experiment.py:
  1. 这个模块用来"猜"RNA-seq的相关配置信息,针对链特异性测序,通过reads的链型转录本的链型来评估reads是哪一条链的.
  2. reads的链型是通过比对结果得到的,转录本的链型是铜鼓注释文件得到的.
  3. 对于非链特异性测序,reads的链型转录本的链型是没有关系的.
  4. 对于链特异性测序,reads的链型转录本的链型是有很大关系的.通过下面的3个例子说明.
  5. 在测序前你并不需要知道RNA-seq是不是链特异性的,就当他们是非链特异性的,这个模块可以"猜"到"reads"是哪条链的.
对于双端RNA-seq,有两种方法来确定reads在哪条链(如illumina ScriptSeq protocol):

(1). 1++,1–,2+-,2-+ :
说明:12表示read1read2,第一个+/-表示read map 到哪条链,第二个+/-表示这个read 所match的基因在哪条链.

那这个1++,1–,2+-,2-+就表示reads所match的链和其所在gene的"+/-"是一样的,也就是reads的链型与其基因的链型一样,是不独立的.

(2). 1+-,1-+,2++,2– :
这个就表示reads所match的链和其所在gene的"+/-"是不一样的,也就是reads的链型与其基因的链型是不独立的.

对于单端测序:

(1). ++,– : 表示read链型与其所match的gene链型一致.
(2). +-,-+ : 表示read链型与其所match的gene链型不一致.


举三个例子说明:

例一:
"猜"是非链特异性测序

解释:reads数的1.72%被映射到基因组区域,我们无法确定这样的gene的链型(比如这个区域两条链都转录)。 对于剩余的98.28%(1 - 0.0172 = 0.9828)的reads,一半是“1 ++,1-,2 + - ,2- +” readsgene链型一致的,而另一半是“1 +1 - +2++,2-”不一致的。 我们得出结论,这不是一个链特异性的测序,因为reads链型不依赖于gene链型.

例二:
"猜"是链特异性测序

解释: 0.72%是不能确定链型的.94.41%readsgene链型一致的.仅有4.87%是不一致的.我们得出结论,这是一个链特异性的测序,因为reads链型依赖于gene链型.

例三,这是个单端测序的例子:
"猜"是链特异性测序

解释: reads链型依赖于gene链型,这个是链特异性测序.


inner_distance.py:

针对双端测序,计算read pairs内部距离或者插入距离,关于插入距离是什么,看下图:

插入距离或内部距离
插入距离D_size计算公式:

D_size = read2_start - read1_end

但是不同条件下计算方法是不一样的:
  • paried reads 比对到同一个外显子: 插入距离 = D_size
  • paried reads比对到不同外显子:插入距离 = D_size - intron_size
  • paried reads比对到非外显子区域(如内含子或者基因外区域): 插入距离 = D_size
  • 如果两个fragments重叠则插入距离可能为负.
结果举例:
插入距离或内部距离

insertion_profile.py:

计算reads上被插入核苷酸的分布.

结果举例:
Read-1 insertion profile
Read-2 insertion profile

junction_annotation.py:

输入一个BAMSAM文件和一个bed12格式的参考基因文件,这个模块将根据参考基因模型计算剪切融合(splice junctions)事件.

  • splice read: 一个RNA read,能够被剪切一次或多次,所以100个spliced reads能够产生>=100个剪切事件.
  • splice junction:多个跨越同一个内含子的剪切事件能够合并为一个splicing junction.
junction 有三种:
  1. 已经被注释的.5'剪切位点3'剪切位点已经被注释.
  2. 全新的.
  3. 部分是新的.
结果示例:
splicing junctions

junction_saturation.py:

在剪切位点分析时首先检查当前的测序深度是否是足够深的.对于一个已经注释的物种,在某个组织中基因的数量是一定的,所以剪切位点数量也是一定的.可以从参考基因模型(bed12格式的文件)可以提前看出splice junctions的数量. 一个饱和的RNA-seq能够发现所有已经被注释的splice junctioons,否则下游剪切位点分析将会出现问题因为将有较少的splice juncitons将发现不了.这个模块从这个测序结果(SAMBAM文件)中进行重抽样,从5%,10%,15%,...,到95%的饱和度来检查每个阶段的splice junctions并与参考基因组进行比较.

结果示例:
不同测序饱和度下的splicing junctions数量

解释: 在这个例子中,每个饱和度下的测序深度的known junctions基本都是一致的(红线).因为大部分的剪切位点我们基本都发现了.即便加深测序深度也不会发现跟多的known junctions,仅会加深junciton的覆盖度(如junction被更多的reads覆盖.).然而当前测序深度(100%的reads)对于发现新的juncitons是不够的(绿线)


mismatch_profile.py:

计算reads的不匹配位点的分布.

结果示例:
reads中的不匹配位点分布图

上图可以看出5'和3'的不匹配位点最多,这是由于测序本身所决定的.


normalize_bigwig.py:

可视化RNA-seq数据结果是最直接并且高效的质控方式.但在可视化之前我们要保证所有样本的数据是可比较的,这就需要进行归一化.信号值文件'wig'或bigwig文件主要由两个因素决定:(1)总reads数,(2)read长度.因此,如果两个样本的read长度不一样但仅对"总reads数"归一化是有问题的.这里我们将每个bigwig文件归一化到相同的wigsum值.wigsum是对基因组信号值的汇总,例如:wigsum=100,000,000等价于1百万个100nt的reads或2百万个50nt的reads的覆盖度. 结果生成wig格式的文件.


overlay_bigwig.py

这个模块让我们操作两个bigwig文件.可以采取的操作有: 信号值相加,取均值,相除,每个信号值+1,求最大值,求最小值,相乘,相减,求几何平均数.


read_distribution.py:

这个模块根据提供的BAM/SAM文件和bed12格式的gene模型文件就按比对上去的reads在基因组上的分布情况,比如在CDS exon,5'UTR exon,intro,基因间区域reads分布.

结果示例:
reads在基因组不同结构上的分布情况

read_duplication.py:

两种用于计算重复率的策略:(1) 基于序列的,完全相同序列的reads被视为重复的reads. (2) 正好map到同一个基因组位置的reads被视为重复reads. 对于splice reads,map到同一位置并且以相同方式剪切的视为是重复reads.

reads重复

read_GC.py:

计算reads的GC含量分布.

reads的GC含量分布

read_hexamer.py:

计算6mer的频率.


read_NVC.py:

这个模块有用来检查核苷酸的碱基组成偏好性.由于随机引物的影响,reads 5'端开始会有某些模式过表达.这种偏好性能够被NVC(Nucleotide versus cycle)画出.理想状态下, A%=C%=G%=T%=25%.

结果示例:
reads每个位置的碱基偏好性

read_quality.py:

可视化reads每个位置的测序质量.

结果示例:

RNA_fragment_size.py:

在map后计算每个gene上的fragment的大小,包括:每个gene上所有的fragment的均值,中位数,方差.

结果示例:
每个gene上的fragment的统计

RPKM_count.py:

这个模块在最新版里已经被废弃,如果想使用可以翻看以前的版本.


RPKM_saturation.py

任何样本统计(RPKM)的精度受样本大小(测序深度)的影响;重抽样切片是使用部分数据来评估样本统计量的精度的方法. 这个模块从总的RNA reads中重抽样并计算每次的RPKM值.通过这样我们就能检测当前测序深度是不是够的(如果测序深度不够RPKM的值将不稳定,如果测序深度足够则RPKM值将稳定).默认情况下,这个模块将计算20个RPKM值(分别是对个转录本使用5%,10%,...,95%的总reads).

在结果图中,Y轴表示 “Percent Relative Error”“Percent Error”.用来表示当前样本量下的RPKM与实际表达量的偏差.计算公式如下:

计算当前样本量下的RPKM与实际RPKM的偏差
结果示例:
不同表达量等级下不同样本量的RPKM偏差

说明:Q1,Q2,Q3,Q4是按照转录本表达量4分位分开的.Q1表示的是表达量低于25%的转录本,以此类推.可以看出:随着样本量升高,RPKM与实际值的偏差也在降低.而且转录本表达量越高这种趋势越明显(Q4最明显).

样本量与RPKM偏差

可以看出,样本量50%之后线条已经趋于平缓,也就是说对于转录本定量来说,当前测序深度是足够的.


spilt_bam.py:

根据提供的bed12注释文件和BAM文件拆分为以下三个文件:

  1. XXX.in.bam: 包含map到外显子趋于的reads.
  2. XXX.ex.bam:包含map不到外显子趋于的reads.
  3. XXX.junk.bam:质控失败或者没有map上去的reads.

split_paired_bam.py:

将一个双端测序BAM文件拆分为两个单端测序BAM文件.


tin.py:

这个模块用来在转录本级别计算RNA完整性TIN (transcript integrity number)值.

结果示例:

更多原创精彩视频敬请关注生信杂谈:

推荐阅读更多精彩内容