本期介紹一個(gè)測(cè)序質(zhì)控的工具包:RSeQC包,它提供了一系列有用的小工具能夠評(píng)估高通量測(cè)序尤其是RNA-seq數(shù)據(jù).比如一些基本模塊;檢查序列質(zhì)量,核酸組分偏性,PCR偏性,GC含量偏性,還有RNA-seq特異性模塊:評(píng)估測(cè)序飽和度，映射讀數(shù)分布，覆蓋均勻性，鏈特異性亥鸠，轉(zhuǎn)錄水平RNA完整性等审丘。

安裝:

RSeQC是依賴于python的,直接使用pip進(jìn)行安裝:

pip install RSeQC

輸入數(shù)據(jù)格式:

RSeQC接受4種文件格式:

• BED格式: Tab分割,12列的表示基因模型的純文本文件,例如:

bed12格式

• SAM或BAM格式: 用來存儲(chǔ)reads比對(duì)結(jié)果信息.SAM是可讀的純文本文件,然而BAM是SAM的二進(jìn)制文本,一個(gè)壓縮的可索引的reads比對(duì)文件. 例如:

SAM格式

• 染色體大小文件: 只有兩列的純文本文件,這個(gè)之前在生物信息學(xué)文本處理大雜燴(一)里已經(jīng)講過.hg19.chrom_24.sizes是人基因組hg19版本的size文件，是使用UCSC 的fetchChromSizes（從這里下載：http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/） 下載的.

人hg19染色體信息

• Fasta文件.

使用方法:

最新版本的RSeQC(2.6.4)包含以下一些模塊,每個(gè)模塊都可以單獨(dú)調(diào)用進(jìn)行分析:

RSeQC包含的模塊

我們一個(gè)一個(gè)來看:

bam2fq.py:

將BAM或SAM格式的文件轉(zhuǎn)為fastq格式.(這個(gè)感覺一般用不到)

bam2wig.py:

將BAM文件轉(zhuǎn)為wig/bigwig格式的文件.(這個(gè)在作圖尤其是信號(hào)圖的時(shí)候很有用.)如果需要轉(zhuǎn)為bigwig格式,則需要UCSC的wigToBigWig工具,下載地址:
http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/wigToBigWig

bam_stat.py:

對(duì)比對(duì)結(jié)果文件BAM或SAM文件進(jìn)行統(tǒng)計(jì).其實(shí)samtools里也有類似工具.結(jié)果如下所示:

bam_stat統(tǒng)計(jì)結(jié)果

統(tǒng)計(jì)結(jié)果包括:總比對(duì)記錄,PCR重復(fù)數(shù),Non Primary Hits表示多匹配位點(diǎn).不匹配的reads數(shù),比對(duì)到+鏈的reads,比對(duì)到-鏈的reads,有剪切位點(diǎn)的reads等.

clipping_profile.py:

這個(gè)模塊用于評(píng)估RNA-seq的BAM或SAM文件中的有切除核苷酸的reads情況.
??clipped reads有兩種,一種是soft-clipped,即reads5'或3'不能比對(duì)到參考基因組;另一種是hard-clipped,即reads5'或3'不能比對(duì)到參考基因組并且被剪切.
??這個(gè)模塊會(huì)生成.r格式的作圖腳本以及.pdf格式的報(bào)告文件以及.xls的數(shù)據(jù)文件.
例如雙端測(cè)序clipping profile圖如下:

Read-1

Read-2

deletion_profile.py:

也就是reads deletion位點(diǎn)的分布.

read deletion distribution

divide_bam.py:

隨機(jī)分割BAM文件(m個(gè)比對(duì)結(jié)果)為n個(gè)文件,每個(gè)文件包含m/n個(gè)比對(duì)結(jié)果.

FPKM_count.py:

根據(jù)read count和gene注釋文件(bed12格式)計(jì)算每個(gè)基因的FPM (fragment per million)或者FPKM(fragment per million mapped reads per kilobase exon).

結(jié)果類似下面的:

FPKM/FPM結(jié)果

geneBody_coverage.py:

計(jì)算RNA-seq reads在基因上的覆蓋度.

方法示意圖

結(jié)果如下:

多個(gè)樣本的RNA-seq reads在基因上的覆蓋度

如果是大于3個(gè)的樣本,則還會(huì)生成熱圖:

reads在基因上的覆蓋度熱圖

geneBody_coverage2.py:

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

infer_experiment.py:

1. 這個(gè)模塊用來"猜"RNA-seq的相關(guān)配置信息,針對(duì)鏈特異性測(cè)序,通過reads的鏈型與轉(zhuǎn)錄本的鏈型來評(píng)估reads是哪一條鏈的.
2. reads的鏈型是通過比對(duì)結(jié)果得到的,轉(zhuǎn)錄本的鏈型是銅鼓注釋文件得到的.
3. 對(duì)于非鏈特異性測(cè)序,reads的鏈型與轉(zhuǎn)錄本的鏈型是沒有關(guān)系的.
4. 對(duì)于鏈特異性測(cè)序,reads的鏈型與轉(zhuǎn)錄本的鏈型是有很大關(guān)系的.通過下面的3個(gè)例子說明.
5. 在測(cè)序前你并不需要知道RNA-seq是不是鏈特異性的,就當(dāng)他們是非鏈特異性的,這個(gè)模塊可以"猜"到"reads"是哪條鏈的.

對(duì)于雙端RNA-seq,有兩種方法來確定reads在哪條鏈(如illumina ScriptSeq protocol):

(1). 1++,1–,2+-,2-+ :
說明:1和2表示read1和read2,第一個(gè)+/-表示read map 到哪條鏈,第二個(gè)+/-表示這個(gè)read 所match的基因在哪條鏈.

那這個(gè)1++,1–,2+-,2-+就表示reads所match的鏈和其所在gene的"+/-"是一樣的,也就是reads的鏈型與其基因的鏈型一樣,是不獨(dú)立的.

(2). 1+-,1-+,2++,2– :
這個(gè)就表示reads所match的鏈和其所在gene的"+/-"是不一樣的,也就是reads的鏈型與其基因的鏈型是不獨(dú)立的.

對(duì)于單端測(cè)序:

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

舉三個(gè)例子說明:

例一:

"猜"是非鏈特異性測(cè)序

解釋: 總reads數(shù)的1.72％被映射到基因組區(qū)域灿巧，我們無法確定這樣的gene的鏈型（比如這個(gè)區(qū)域兩條鏈都轉(zhuǎn)錄）赂摆。 對(duì)于剩余的98.28％（1 - 0.0172 = 0.9828）的reads，一半是“1 ++狸窘，1-墩朦，2 + - ，2- +” reads與gene鏈型一致的翻擒，而另一半是“1 +1 - +2++氓涣，2-”不一致的。 我們得出結(jié)論韭寸，這不是一個(gè)鏈特異性的測(cè)序，因?yàn)?code>reads鏈型不依賴于gene鏈型.

例二:

"猜"是鏈特異性測(cè)序

解釋: 0.72%是不能確定鏈型的.94.41%的reads與gene鏈型一致的.僅有4.87%是不一致的.我們得出結(jié)論荆隘，這是一個(gè)鏈特異性的測(cè)序恩伺，因?yàn)?code>reads鏈型依賴于gene鏈型.

例三,這是個(gè)單端測(cè)序的例子:

"猜"是鏈特異性測(cè)序

解釋: reads鏈型依賴于gene鏈型,這個(gè)是鏈特異性測(cè)序.

inner_distance.py:

針對(duì)雙端測(cè)序,計(jì)算read pairs的內(nèi)部距離或者插入距離,關(guān)于插入距離是什么,看下圖:

插入距離或內(nèi)部距離

插入距離D_size計(jì)算公式:

D_size = read2_start - read1_end

但是不同條件下計(jì)算方法是不一樣的:

• paried reads 比對(duì)到同一個(gè)外顯子: 插入距離 = D_size
• paried reads比對(duì)到不同外顯子:插入距離 = D_size - intron_size
• paried reads比對(duì)到非外顯子區(qū)域(如內(nèi)含子或者基因外區(qū)域): 插入距離 = D_size
• 如果兩個(gè)fragments重疊則插入距離可能為負(fù).

結(jié)果舉例:

插入距離或內(nèi)部距離

insertion_profile.py:

計(jì)算reads上被插入核苷酸的分布.

結(jié)果舉例:

Read-1 insertion profile

Read-2 insertion profile

junction_annotation.py:

輸入一個(gè)BAM或SAM文件和一個(gè)bed12格式的參考基因文件,這個(gè)模塊將根據(jù)參考基因模型計(jì)算剪切融合(splice junctions)事件.

• splice read: 一個(gè)RNA read,能夠被剪切一次或多次,所以100個(gè)spliced reads能夠產(chǎn)生>=100個(gè)剪切事件.
• splice junction:多個(gè)跨越同一個(gè)內(nèi)含子的剪切事件能夠合并為一個(gè)splicing junction.

junction 有三種:

1. 已經(jīng)被注釋的.5'剪切位點(diǎn)和3'剪切位點(diǎn)已經(jīng)被注釋.
2. 全新的.
3. 部分是新的.

結(jié)果示例:

splicing junctions

junction_saturation.py:

在剪切位點(diǎn)分析時(shí)首先檢查當(dāng)前的測(cè)序深度是否是足夠深的.對(duì)于一個(gè)已經(jīng)注釋的物種,在某個(gè)組織中基因的數(shù)量是一定的,所以剪切位點(diǎn)數(shù)量也是一定的.可以從參考基因模型(bed12格式的文件)可以提前看出splice junctions的數(shù)量. 一個(gè)飽和的RNA-seq能夠發(fā)現(xiàn)所有已經(jīng)被注釋的splice junctioons,否則下游剪切位點(diǎn)分析將會(huì)出現(xiàn)問題因?yàn)閷⒂休^少的splice juncitons將發(fā)現(xiàn)不了.這個(gè)模塊從這個(gè)測(cè)序結(jié)果(SAM或BAM文件)中進(jìn)行重抽樣,從5%,10%,15%,...,到95%的飽和度來檢查每個(gè)階段的splice junctions并與參考基因組進(jìn)行比較.

結(jié)果示例:

不同測(cè)序飽和度下的splicing junctions數(shù)量

解釋: 在這個(gè)例子中,每個(gè)飽和度下的測(cè)序深度的known junctions基本都是一致的(紅線).因?yàn)榇蟛糠值募羟形稽c(diǎn)我們基本都發(fā)現(xiàn)了.即便加深測(cè)序深度也不會(huì)發(fā)現(xiàn)跟多的known junctions,僅會(huì)加深junciton的覆蓋度(如junction被更多的reads覆蓋.).然而當(dāng)前測(cè)序深度(100%的reads)對(duì)于發(fā)現(xiàn)新的juncitons是不夠的(綠線)

mismatch_profile.py:

計(jì)算reads的不匹配位點(diǎn)的分布.

結(jié)果示例:

reads中的不匹配位點(diǎn)分布圖

上圖可以看出5'和3'的不匹配位點(diǎn)最多,這是由于測(cè)序本身所決定的.

normalize_bigwig.py:

可視化RNA-seq數(shù)據(jù)結(jié)果是最直接并且高效的質(zhì)控方式.但在可視化之前我們要保證所有樣本的數(shù)據(jù)是可比較的,這就需要進(jìn)行歸一化.信號(hào)值文件'wig'或bigwig文件主要由兩個(gè)因素決定:(1)總reads數(shù),(2)read長度.因此,如果兩個(gè)樣本的read長度不一樣但僅對(duì)"總reads數(shù)"歸一化是有問題的.這里我們將每個(gè)bigwig文件歸一化到相同的wigsum值.wigsum是對(duì)基因組信號(hào)值的匯總,例如:wigsum=100,000,000等價(jià)于1百萬個(gè)100nt的reads或2百萬個(gè)50nt的reads的覆蓋度. 結(jié)果生成wig格式的文件.

overlay_bigwig.py

這個(gè)模塊讓我們操作兩個(gè)bigwig文件.可以采取的操作有: 信號(hào)值相加,取均值,相除,每個(gè)信號(hào)值+1,求最大值,求最小值,相乘,相減,求幾何平均數(shù).

read_distribution.py:

這個(gè)模塊根據(jù)提供的BAM/SAM文件和bed12格式的gene模型文件就按比對(duì)上去的reads在基因組上的分布情況,比如在CDS exon,5'UTR exon,intro,基因間區(qū)域的reads分布.

結(jié)果示例:

reads在基因組不同結(jié)構(gòu)上的分布情況

read_duplication.py:

兩種用于計(jì)算重復(fù)率的策略:(1) 基于序列的,完全相同序列的reads被視為重復(fù)的reads. (2) 正好map到同一個(gè)基因組位置的reads被視為重復(fù)reads. 對(duì)于splice reads,map到同一位置并且以相同方式剪切的視為是重復(fù)reads.

reads重復(fù)

read_GC.py:

計(jì)算reads的GC含量分布.

reads的GC含量分布

read_hexamer.py:

計(jì)算6mer的頻率.

read_NVC.py:

這個(gè)模塊有用來檢查核苷酸的堿基組成偏好性.由于隨機(jī)引物的影響,reads 5'端開始會(huì)有某些模式過表達(dá).這種偏好性能夠被NVC(Nucleotide versus cycle)畫出.理想狀態(tài)下, A%=C%=G%=T%=25%.

結(jié)果示例:

reads每個(gè)位置的堿基偏好性

read_quality.py:

可視化reads每個(gè)位置的測(cè)序質(zhì)量.

結(jié)果示例:

RNA_fragment_size.py:

在map后計(jì)算每個(gè)gene上的fragment的大小,包括:每個(gè)gene上所有的fragment的均值,中位數(shù),方差.

結(jié)果示例:

每個(gè)gene上的fragment的統(tǒng)計(jì)

RPKM_count.py:

這個(gè)模塊在最新版里已經(jīng)被廢棄,如果想使用可以翻看以前的版本.

RPKM_saturation.py

任何樣本統(tǒng)計(jì)（RPKM）的精度受樣本大小（測(cè)序深度）的影響;重抽樣或切片是使用部分?jǐn)?shù)據(jù)來評(píng)估樣本統(tǒng)計(jì)量的精度的方法. 這個(gè)模塊從總的RNA reads中重抽樣并計(jì)算每次的RPKM值.通過這樣我們就能檢測(cè)當(dāng)前測(cè)序深度是不是夠的(如果測(cè)序深度不夠RPKM的值將不穩(wěn)定,如果測(cè)序深度足夠則RPKM值將穩(wěn)定).默認(rèn)情況下,這個(gè)模塊將計(jì)算20個(gè)RPKM值(分別是對(duì)個(gè)轉(zhuǎn)錄本使用5%,10%,...,95%的總reads).

在結(jié)果圖中,Y軸表示 “Percent Relative Error” 或 “Percent Error”.用來表示當(dāng)前樣本量下的RPKM與實(shí)際表達(dá)量的偏差.計(jì)算公式如下:

計(jì)算當(dāng)前樣本量下的RPKM與實(shí)際RPKM的偏差

結(jié)果示例:

不同表達(dá)量等級(jí)下不同樣本量的RPKM偏差

說明:Q1,Q2,Q3,Q4是按照轉(zhuǎn)錄本表達(dá)量4分位分開的.Q1表示的是表達(dá)量低于25%的轉(zhuǎn)錄本,以此類推.可以看出:隨著樣本量升高,RPKM與實(shí)際值的偏差也在降低.而且轉(zhuǎn)錄本表達(dá)量越高這種趨勢(shì)越明顯(Q4最明顯).

樣本量與RPKM偏差

可以看出,樣本量50%之后線條已經(jīng)趨于平緩,也就是說對(duì)于轉(zhuǎn)錄本定量來說,當(dāng)前測(cè)序深度是足夠的.

spilt_bam.py:

根據(jù)提供的bed12注釋文件和BAM文件拆分為以下三個(gè)文件:

1. XXX.in.bam: 包含map到外顯子趨于的reads.
2. XXX.ex.bam:包含map不到外顯子趨于的reads.
3. XXX.junk.bam:質(zhì)控失敗或者沒有map上去的reads.

split_paired_bam.py:

將一個(gè)雙端測(cè)序的BAM文件拆分為兩個(gè)單端測(cè)序的BAM文件.

tin.py:

這個(gè)模塊用來在轉(zhuǎn)錄本級(jí)別計(jì)算RNA完整性TIN (transcript integrity number)值.

結(jié)果示例:

更多原創(chuàng)精彩視頻敬請(qǐng)關(guān)注生信雜談: