本期介绍一个测序质控的工具包:RSeQC包,它提供了一系列有用的小工具能够评估高通量测序尤其是RNA-seq数据.比如一些基本模块;
检查序列质量,核酸组分偏性,PCR偏性,GC含量偏性,还有RNA-seq特异性模块:评估测序饱和度，映射读数分布，覆盖均匀性，链特异性，转录水平RNA完整性等。
安装:
RSeQC是依赖于python的,直接使用pip进行安装:
pip install RSeQC
输入数据格式:
RSeQC接受4种文件格式:
-
BED格式:Tab分割,12列的表示基因模型的纯文本文件,例如:
-
SAM或BAM格式: 用来存储reads比对结果信息.SAM是可读的纯文本文件,然而BAM是SAM的二进制文本,一个压缩的可索引的reads比对文件. 例如:
- 染色体大小文件: 只有两列的纯文本文件,这个之前在生物信息学文本处理大杂烩(一)里已经讲过.
hg19.chrom_24.sizes是人基因组hg19版本的size文件，是使用UCSC 的fetchChromSizes（从这里下载：http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/） 下载的.
- Fasta文件.
使用方法:
最新版本的RSeQC(2.6.4)包含以下一些模块,每个模块都可以单独调用进行分析:
我们一个一个来看:
bam2fq.py:
将BAM或SAM格式的文件转为fastq格式.(这个感觉一般用不到)
bam2wig.py:
将BAM文件转为wig/bigwig格式的文件.(这个在作图尤其是信号图的时候很有用.)如果需要转为bigwig格式,则需要UCSC的wigToBigWig工具,下载地址:
http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/wigToBigWig
bam_stat.py:
对比对结果文件BAM或SAM文件进行统计.其实samtools里也有类似工具.结果如下所示:
统计结果包括:总比对记录,PCR重复数,Non Primary Hits表示多匹配位点.不匹配的reads数,比对到+链的reads,比对到-链的reads,有剪切位点的reads等.
clipping_profile.py:
这个模块用于评估RNA-seq的BAM或SAM文件中的有切除核苷酸的reads情况.
clipped reads有两种,一种是soft-clipped,即reads5'或3'不能比对到参考基因组;另一种是hard-clipped,即reads5'或3'不能比对到参考基因组并且被剪切.
这个模块会生成.r格式的作图脚本以及.pdf格式的报告文件以及.xls的数据文件.
例如双端测序clipping profile图如下:
deletion_profile.py:
也就是reads deletion位点的分布.
divide_bam.py:
随机分割BAM文件(m个比对结果)为n个文件,每个文件包含m/n个比对结果.
FPKM_count.py:
根据read count和gene注释文件(bed12格式)计算每个基因的FPM (fragment per million)或者FPKM(fragment per million mapped reads per kilobase exon).
结果类似下面的:
geneBody_coverage.py:
计算RNA-seq reads在基因上的覆盖度.
结果如下:
如果是大于3个的样本,则还会生成热图:
geneBody_coverage2.py:
功能和上面的geneBody_coverage.py一样,但输入的是bigwig格式文件
infer_experiment.py:
- 这个模块用来"猜"RNA-seq的相关配置信息,针对链特异性测序,通过
reads的链型与转录本的链型来评估reads是哪一条链的. -
reads的链型是通过比对结果得到的,转录本的链型是铜鼓注释文件得到的. - 对于非链特异性测序,
reads的链型与转录本的链型是没有关系的. - 对于链特异性测序,
reads的链型与转录本的链型是有很大关系的.通过下面的3个例子说明. - 在测序前你并不需要知道RNA-seq是不是链特异性的,就当他们是非链特异性的,这个模块可以"猜"到"reads"是哪条链的.
对于双端RNA-seq,有两种方法来确定reads在哪条链(如illumina ScriptSeq protocol):
(1). 1++,1–,2+-,2-+ :
说明:1和2表示read1和read2,第一个+/-表示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- +” reads与gene链型一致的，而另一半是“1 +1 - +2++，2-”不一致的。 我们得出结论，这不是一个链特异性的测序，因为reads链型不依赖于gene链型.
例二:
解释: 0.72%是不能确定链型的.94.41%的reads与gene链型一致的.仅有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上被插入核苷酸的分布.
结果举例:
junction_annotation.py:
输入一个BAM或SAM文件和一个bed12格式的参考基因文件,这个模块将根据参考基因模型计算剪切融合(splice junctions)事件.
-
splice read: 一个RNA read,能够被剪切一次或多次,所以100个
spliced reads能够产生>=100个剪切事件. -
splice junction:多个跨越同一个内含子的剪切事件能够合并为一个
splicing junction.
junction 有三种:
- 已经被注释的.
5'剪切位点和3'剪切位点已经被注释. - 全新的.
- 部分是新的.
结果示例:
junction_saturation.py:
在剪切位点分析时首先检查当前的测序深度是否是足够深的.对于一个已经注释的物种,在某个组织中基因的数量是一定的,所以剪切位点数量也是一定的.可以从参考基因模型(bed12格式的文件)可以提前看出splice junctions的数量. 一个饱和的RNA-seq能够发现所有已经被注释的splice junctioons,否则下游剪切位点分析将会出现问题因为将有较少的splice juncitons将发现不了.这个模块从这个测序结果(SAM或BAM文件)中进行重抽样,从5%,10%,15%,...,到95%的饱和度来检查每个阶段的splice junctions并与参考基因组进行比较.
结果示例:
解释: 在这个例子中,每个饱和度下的测序深度的known junctions基本都是一致的(红线).因为大部分的剪切位点我们基本都发现了.即便加深测序深度也不会发现跟多的known junctions,仅会加深junciton的覆盖度(如junction被更多的reads覆盖.).然而当前测序深度(100%的reads)对于发现新的juncitons是不够的(绿线)
mismatch_profile.py:
计算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分布.
结果示例:
read_duplication.py:
两种用于计算重复率的策略:(1) 基于序列的,完全相同序列的reads被视为重复的reads. (2) 正好map到同一个基因组位置的reads被视为重复reads. 对于splice reads,map到同一位置并且以相同方式剪切的视为是重复reads.
read_GC.py:
计算reads的GC含量分布.
read_hexamer.py:
计算6mer的频率.
read_NVC.py:
这个模块有用来检查核苷酸的碱基组成偏好性.由于随机引物的影响,reads 5'端开始会有某些模式过表达.这种偏好性能够被NVC(Nucleotide versus cycle)画出.理想状态下, A%=C%=G%=T%=25%.
结果示例:
read_quality.py:
可视化reads每个位置的测序质量.
结果示例:
RNA_fragment_size.py:
在map后计算每个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与实际表达量的偏差.计算公式如下:
结果示例:
说明:Q1,Q2,Q3,Q4是按照转录本表达量4分位分开的.Q1表示的是表达量低于25%的转录本,以此类推.可以看出:随着样本量升高,RPKM与实际值的偏差也在降低.而且转录本表达量越高这种趋势越明显(Q4最明显).
可以看出,样本量50%之后线条已经趋于平缓,也就是说对于转录本定量来说,当前测序深度是足够的.
spilt_bam.py:
根据提供的bed12注释文件和BAM文件拆分为以下三个文件:
- XXX.in.bam: 包含map到外显子趋于的reads.
- XXX.ex.bam:包含map不到外显子趋于的reads.
- XXX.junk.bam:质控失败或者没有map上去的reads.
split_paired_bam.py:
将一个双端测序的BAM文件拆分为两个单端测序的BAM文件.
tin.py:
这个模块用来在转录本级别计算RNA完整性TIN (transcript integrity number)值.
结果示例:
更多原创精彩视频敬请关注生信杂谈:
