當(dāng)前位置：首頁(yè) > 编程资源 > 编程问答 >内容正文

编程问答

转录组入门（5）：序列比对

發(fā)布時(shí)間：2023/12/20 编程问答 28 豆豆

生活随笔收集整理的這篇文章主要介紹了转录组入门（5）：序列比对小編覺(jué)得挺不錯(cuò)的,現(xiàn)在分享給大家,幫大家做個(gè)參考.

歡迎來(lái)GitHub上fork，一起進(jìn)步： https://github.com/xuzhougeng

比對(duì)軟件很多，首先大家去收集一下，因?yàn)槲覀兪菐Т蠹胰腴T，請(qǐng)統(tǒng)一用hisat2，并且搞懂它的用法。
直接去hisat2的主頁(yè)下載index文件即可，然后把fastq格式的reads比對(duì)上去得到sam文件。
接著用samtools把它轉(zhuǎn)為bam文件，并且排序(注意N和P兩種排序區(qū)別)索引好，載入IGV，再截圖幾個(gè)基因看看！
順便對(duì)bam文件進(jìn)行簡(jiǎn)單QC，參考直播我的基因組系列。

前面四篇基本都算是準(zhǔn)備工作，從這一篇開(kāi)始才算進(jìn)入了RNA-Seq數(shù)據(jù)分析的核心部分。

比對(duì)

比對(duì)還是不比對(duì)

在比對(duì)之前，我們得了解比對(duì)的目的是什么？RNA-Seq數(shù)據(jù)比對(duì)和DNA-Seq數(shù)據(jù)比對(duì)有什么差異？
RNA-Seq數(shù)據(jù)分析分為很多種，比如說(shuō)找差異表達(dá)基因或?qū)ふ倚碌目勺兗羟小Ｈ绻也町惐磉_(dá)基因單純只需要確定不同的read計(jì)數(shù)就行的話，我們可以用bowtie, bwa這類比對(duì)工具，或者是salmon這類align-free工具，并且后者的速度更快。

但是如果你需要找到新的isoform，或者RNA的可變剪切，看看外顯子使用差異的話，你就需要TopHat, HISAT2或者是STAR這類工具用于找到剪切位點(diǎn)。因?yàn)镽NA-Seq不同于DNA-Seq，DNA在轉(zhuǎn)錄成mRNA的時(shí)候會(huì)把內(nèi)含子部分去掉。所以mRNA反轉(zhuǎn)的cDNA如果比對(duì)不到參考序列，會(huì)被分開(kāi)，重新比對(duì)一次，判斷中間是否有內(nèi)含子。

工具抉擇

在2016年的一篇綜述A survey of best practices for RNA-seq data analysis，提到目前有三種RNA數(shù)據(jù)分析的策略。那個(gè)時(shí)候的工具也主要用的是TopHat,STAR和Bowtie.其中TopHat目前已經(jīng)被它的作者推薦改用HISAT進(jìn)行替代。

最近的Nature Communication發(fā)表了一篇題為的Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis的文章--被稱之為史上最全RNA-Seq數(shù)據(jù)分析流程，也是我一直以來(lái)想做的事情，只不過(guò)他們做的超乎我的想象。文章中在基于參考基因組的轉(zhuǎn)錄本分析中所用的工具，是TopHat,HISAT2和STAR，結(jié)論就是HISAT2找到j(luò)unction正確率最高，但是在總數(shù)上卻比TopHat和STAR少。從這里可以看出HISAT2的二類錯(cuò)誤(納偽）比較少，但是一類錯(cuò)誤（棄真）就高起來(lái)。
就唯一比對(duì)而言，STAR是三者最佳的，主要是因?yàn)樗粫?huì)像TopHat和HISAT2一樣在PE比對(duì)不上的情況還強(qiáng)行把SE也比對(duì)到基因組上。而且在處理較長(zhǎng)的read和較短read的不同情況，STAR的穩(wěn)定性也是最佳的。
就速度而言，HISAT2比STAR和TopHat2平均快上2.5~100倍。

如果學(xué)習(xí)RNA-Seq數(shù)據(jù)分析，上面提到的兩篇文獻(xiàn)是必須要看上3遍以上的，而且建議每隔一段時(shí)間回顧一下。但是如果就比對(duì)工具而言，基本上就是HISAT2和STAR選一個(gè)就行。

下載index

首先，問(wèn)自己一個(gè)問(wèn)題，為什么比對(duì)的時(shí)候需要用到index？這里強(qiáng)烈建議大家去看Jimmy寫的bowtie算法原理探究bowtie算法原理探究。但是只是建議，你不需要真的去看，反正你也看不懂。

高通量測(cè)序遇到的第一個(gè)問(wèn)題就是，成千上萬(wàn)甚至上幾億條read如果在合理的時(shí)間內(nèi)比對(duì)到參考基因組上，并且保證錯(cuò)誤率在接受范圍內(nèi)。為了提高比對(duì)速度，就需要根據(jù)參考基因組序列，經(jīng)過(guò)BWT算法轉(zhuǎn)換成index，而我們比對(duì)的序列其實(shí)是index的一個(gè)子集。當(dāng)然轉(zhuǎn)錄組比對(duì)還要考慮到可變剪切的情況，所以更加復(fù)雜。

因此我門不是直接把read回貼到基因組上，而是把read和index進(jìn)行比較。人類的index一般都是有現(xiàn)成的，我建議大家下載現(xiàn)成的，我曾經(jīng)嘗試過(guò)用服務(wù)器自己創(chuàng)建index，花的時(shí)間讓我懷疑人生。

cd referece && mkdir index && cd index wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz tar -zxvf hg19.tar.gz

覺(jué)得電腦配置還行的，或者是沒(méi)有現(xiàn)成index的，可以通過(guò)HISAT2的方法進(jìn)行創(chuàng)建

# 其實(shí)hisat2-buld在運(yùn)行的時(shí)候也會(huì)自己尋找exons和splice_sites，但是先做的目的是為了提高運(yùn)行效率 extract_exons.py gencode.v26lift37.annotation.sorted.gtf > hg19.exons.gtf & extract_splice_sites.py gencode.v26lift37.annotation.gtf > hg19.splice_sites.gtf & # 建立index，必須選項(xiàng)是基因組所在文件路徑和輸出的前綴 hisat2-build --ss hg19.splice_sites.gtf --exon hg19.exons.gtf genome/hg19/hg19.fa hg19

我的是i7-7700處理器，內(nèi)存是64G，運(yùn)行的資源效率如下：

正式比對(duì)

hisat2基本用法就是hisat2 [options]* -x <hisat2-idx> {-1 <m1> -2 <m2> | -U <r> } [-S <hit>]，基本就是提供index的位置，PE數(shù)據(jù)或者是SE數(shù)據(jù)存放位置。然而其他可選參數(shù)卻是進(jìn)階的一大名堂。新手就用默認(rèn)參數(shù)唄。

因?yàn)镽NA-Seq數(shù)據(jù)是從 SRR3589957 ~ SRR3589962，6個(gè)樣本的PE數(shù)據(jù)，也就是有6次循環(huán)，可以寫腳本，也可以直接在命令行里運(yùn)行
如下命令運(yùn)行所在目錄為/mnt/f/Data/，我的參考基因組的index數(shù)據(jù)存放在/mnt/f/Data/reference/index/hg19/，而RNA-seq數(shù)據(jù)存放在/mnt/f/Data/RNA-Seq下。比對(duì)結(jié)果會(huì)存放在/mnt/f/Data/RNA-Seq/aligned

mkdir -p RNA-Seq/aligned for i in `seq 57 62` dohisat2 -t -x reference/index/hg19/genome -1 RNA-Seq/SRR35899${i}_1.fastq.gz -2 SRR35899${i}_2.fastq.gz -S RNA-Seq/aligned/SRR35899${i}.sam & done

&會(huì)把任務(wù)丟到后臺(tái)，所以會(huì)同時(shí)執(zhí)行這3個(gè)比對(duì)程序，如果CPU和內(nèi)存承受不住，去掉&一個(gè)個(gè)來(lái)。比對(duì)這一步是非常消耗內(nèi)存資源的，這是比對(duì)工具要將索引數(shù)據(jù)放入內(nèi)存引起的。我有64G內(nèi)存，理論上可以同時(shí)處理20個(gè)PE數(shù)據(jù)。在我的電腦配置下，大致花了2個(gè)小時(shí)同時(shí)才完成這一步.

基本參數(shù)說(shuō)明

在數(shù)據(jù)比對(duì)的時(shí)候，可以安靜一下讀讀HISAT2的額外選項(xiàng)，主要分為如下幾塊

主要參數(shù)，一定要填寫的內(nèi)容
輸入選項(xiàng)，對(duì)結(jié)果影響不大
比對(duì)選項(xiàng)，主要是--n-ceil決定模糊字符的數(shù)量
得分選項(xiàng)，當(dāng)一個(gè)read比對(duì)到不同部位時(shí)，確定那個(gè)才是最優(yōu)的。基于mismatch, soft-cliping, gap得分。
可變剪切比對(duì)選項(xiàng)，你要決定exon，intron的長(zhǎng)度，GT/AG的得分，還可以提供已知的可變剪切和外顯子gtf文件，
報(bào)告選項(xiàng)，確定要找多少的位置
PE選項(xiàng)，與gap有關(guān)的參數(shù)
輸出選項(xiàng)，建議加上-t記錄時(shí)間，其他就是壓縮格式，不影響比對(duì)
SAM選項(xiàng)，主要是決定SAM的header應(yīng)該添加哪些內(nèi)容
性能選項(xiàng)和其他選項(xiàng)不考慮

注: soft clipping 指的是比對(duì)的read只有部分匹配到參考序列上，還有部分沒(méi)有匹配上。也就是一個(gè)100bp的read，就匹配上前面20 bp或者是后面20bp，或者是后面20bp比對(duì)的效果不太好。

因此影響比對(duì)結(jié)果就是比對(duì)選項(xiàng)，得分選項(xiàng)，可變剪切選項(xiàng)和PE選項(xiàng)，在有生之年我應(yīng)該會(huì)寫一片文章介紹這些選項(xiàng)對(duì)結(jié)果的影響。

HISAT2輸出結(jié)果

比對(duì)之后會(huì)輸出如下結(jié)果，解讀一下就是全部數(shù)據(jù)都是100%的，96.68%的配對(duì)數(shù)據(jù)一次都沒(méi)有比對(duì)，1.23%的數(shù)據(jù)比是唯一比對(duì)，2.09%是多個(gè)比對(duì)。然后96.68%一次都沒(méi)有比對(duì)的數(shù)據(jù)，如果不按照順序來(lái)，有0.05%的比對(duì)。之后把剩下的部分用單端數(shù)據(jù)進(jìn)行比對(duì)的話，95.20%數(shù)據(jù)沒(méi)比對(duì)上，3.60%的數(shù)據(jù)比對(duì)一次，1.20%比對(duì)超過(guò)一次。零零總總的加起來(lái)是8%的比對(duì)！！！

這個(gè)總體比對(duì)率讓我開(kāi)始懷疑人生，怎么可能呀，我翻了翻輸出記錄，發(fā)現(xiàn)有幾個(gè)結(jié)果的比對(duì)率超過(guò)90%呀。我思索了片刻，驚醒這個(gè)實(shí)驗(yàn)好像是用人類和小鼠都做了一遍。于是又去GEO上查了一下記錄，恍然大悟，差點(diǎn)翻車。

Samples 9-15 are mRNA-seq to determine effect of AKAP95 knockdown in human 293 cells (9-11) or mouse ES cells (12-15).

同時(shí)我反思了一下出錯(cuò)的原因，我默認(rèn)這個(gè)實(shí)驗(yàn)是KO和非KO各3個(gè)重復(fù)，其實(shí)文章的實(shí)驗(yàn)設(shè)計(jì)并不是如此，可見(jiàn)理解實(shí)驗(yàn)設(shè)計(jì)很重要，于是我把數(shù)據(jù)下載這一部分進(jìn)行了完善。

mkdir -p RNA-Seq/aligned for i in `seq 56 58` dohisat2 -t -x reference/index/hg19/genome -1 RNA-Seq/SRR35899${i}_1.fastq.gz -2 SRR35899${i}_2.fastq.gz -S RNA-Seq/SRR35899${i}.sam & done

如上是修改后的代碼

SAMtools三板斧

SAM（sequence Alignment/mapping)數(shù)據(jù)格式是目前高通量測(cè)序中存放比對(duì)數(shù)據(jù)的標(biāo)準(zhǔn)格式，當(dāng)然他可以用于存放未比對(duì)的數(shù)據(jù)。所以，SAM的格式說(shuō)明

而目前處理SAM格式的工具主要是SAMTools，這是Heng Li大神寫的.除了C語(yǔ)言版本，還有Java的Picard，Python的Pysam，Common lisp的cl-sam等其他版本。SAMTools的主要功能如下：

view: BAM-SAM/SAM-BAM 轉(zhuǎn)換和提取部分比對(duì)
sort: 比對(duì)排序
merge: 聚合多個(gè)排序比對(duì)
index: 索引排序比對(duì)
faidx: 建立FASTA索引，提取部分序列
tview: 文本格式查看序列
pileup: 產(chǎn)生基于位置的結(jié)果和 consensus/indel calling

最常用的三板斧就是格式轉(zhuǎn)換，排序，索引。而進(jìn)階教程就是看文檔提高。

for i in `seq 56 58` dosamtools view -S SRR35899${i}.sam -b > SRR35899${i}.bamsamtools sort SRR35899${i}.bam -o SRR35899${i}_sorted.bamsamtools index SRR35899${i}_sorted.bam done

注

-S是最新版samtools為了兼容以前版本寫的，所以可以省去
0.1.19版本和最新版有比較大差別，請(qǐng)注意版本

Jimmy說(shuō)樣我們仔細(xì)判斷sam排序兩種方式的不同，因此我截取前面100行數(shù)據(jù)，分別排序然后查看結(jié)果。

head -1000 SRR3589957.sam > test.sam samtools view -b test.sam > test.bam samtools view test.bam | head

默認(rèn)排序是根據(jù)染色體位置

samtools sort test.bam default samtools view default.bam | head

Sort alignments by leftmost coordinates, or by read name when -n is used

samtools sort -n test.bam sort_left samtools view sort_left.bam | head

也就說(shuō)說(shuō)默認(rèn)按照染色體位置進(jìn)行排序，而-n參數(shù)則是根據(jù)read名進(jìn)行排序。當(dāng)然還有一個(gè)-t 根據(jù)TAG進(jìn)行排序。

說(shuō)說(shuō)samtools view

三板斧的view是一個(gè)非常實(shí)用的子命令，除了之前的格式轉(zhuǎn)換以外，還能進(jìn)行數(shù)據(jù)提取和提取。
比如說(shuō)提取1號(hào)染色體1234-123456區(qū)域的比對(duì)read

samtools view SRR3589957_sorted.bam chr1:1234-123456 | head

在比如搭配flag(0.1.19版本沒(méi)有）和flagstat，使用-f或-F參數(shù)提取不同匹配情況的read。
flag是一種描述read比對(duì)情況的標(biāo)記，一種12種，可以搭配使用。

0x1 PAIRED paired-end (or multiple-segment) sequencing technology 0x2 PROPER_PAIR each segment properly aligned according to the aligner 0x4 UNMAP segment unmapped 0x8 MUNMAP next segment in the template unmapped 0x10 REVERSE SEQ is reverse complemented 0x20 MREVERSE SEQ of the next segment in the template is reverse complemented 0x40 READ1 the first segment in the template 0x80 READ2 the last segment in the template 0x100 SECONDARY secondary alignment 0x200 QCFAIL not passing quality controls 0x400 DUP PCR or optical duplicate 0x800 SUPPLEMENTARY supplementary alignment

可以先用flagstat看下總體情況

samtools flagstat SRR3589957_sorted.bam

也就是說(shuō)如果我想用samtools篩選恰好配對(duì)的read,就需要用0x10

samtools view -b -f 0x10 SRR3589957_sorted.bam chr1:1234-123456 > flag.bam samtools flagstat flag.bam

我應(yīng)該會(huì)在有生之年寫一篇文章好好介紹samtools。

比對(duì)質(zhì)控(QC)

還是在A survey of best practices for RNA-seq data analysis里面，提到了人類基因組應(yīng)該有70%~90%的比對(duì)率，并且多比對(duì)read（multi-mapping reads）數(shù)量要少。另外比對(duì)在外顯子和所比對(duì)鏈（uniformity of read coverage on exons and the mapped strand）的覆蓋度要保持一致。

常用工具有

Picard https://broadinstitute.github.io/picard/
RSeQC http://rseqc.sourceforge.net/
Qualimap http://qualimap.bioinfo.cipf.es/

我們就用RSeQC吧，畢竟使用python寫的工具，天生的親切感，而且安裝非常方便。

# Python2.7環(huán)境下 pip install RSeQC

一共有如下幾個(gè)文件，根據(jù)命名就知道功能是啥了。

bam2fq.py
bam2wig.py
bam_stat.py
clipping_profile.py
deletion_profile.py
divide_bam.py
FPKM_count.py
geneBody_coverage.py
geneBody_coverage2.py
infer_experiment.py
inner_distance.py
insertion_profile.py
junction_annotation.py
junction_saturation.py
mismatch_profile.py
normalize_bigwig.py
overlay_bigwig.py
read_distribution.py
read_duplication.py
read_GC.py
read_hexamer.py
read_NVC.py
read_quality.py
RNA_fragment_size.py
RPKM_count.py
RPKM_saturation.py
spilt_bam.py
split_paired_bam.py
tin.py

先對(duì)bam文件進(jìn)行統(tǒng)計(jì)分析，從結(jié)果上看是符合70~90的比對(duì)率要求。

bam_stat.py -i SRR3589956_sorted.bam

基因組覆蓋率的QC需要提供bed文件，可以直接RSeQC的網(wǎng)站下載，或者可以用gtf轉(zhuǎn)換

read_distribution.py -i RNA-Seq/aligned/SRR3589956_sorted.bam -r reference/hg19_RefSeq.bed

IGV查看

載入?yún)⒖夹蛄?#xff0c;注釋和BAM文件，隨便看看吧。

最后請(qǐng)?jiān)试S我給自己的小密圈打個(gè)廣告，如果你非常支持我的創(chuàng)作，想和我保持聯(lián)系的話。

總結(jié)

以上是生活随笔為你收集整理的转录组入门（5）：序列比对的全部?jī)?nèi)容，希望文章能夠幫你解決所遇到的問(wèn)題。

如果覺(jué)得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。