轉錄組入門（5）： 序列比對

歡迎來GitHub上fork，一起進步： https://github.com/xuzhougeng

比對軟體很多，首先大家去收集一下，因為我們是帶大家入門，請統一用hisat2，并且搞懂它的用法。

直接去hisat2的首頁下載下傳index檔案即可，然後把fastq格式的reads比對上去得到sam檔案。

接着用samtools把它轉為bam檔案，并且排序(注意N和P兩種排序差別)索引好，載入IGV，再截圖幾個基因看看！

順便對bam檔案進行簡單QC，參考直播我的基因組系列。

前面四篇基本都算是準備工作，從這一篇開始才算進入了RNA-Seq資料分析的核心部分。

比對

比對還是不比對

在比對之前，我們得了解比對的目的是什麼？RNA-Seq資料比對和DNA-Seq資料比對有什麼差異？

RNA-Seq資料分析分為很多種，比如說找差異表達基因或尋找新的可變剪切。如果找差異表達基因單純隻需要确定不同的read計數就行的話，我們可以用bowtie, bwa這類比對工具，或者是salmon這類align-free工具，并且後者的速度更快。

但是如果你需要找到新的isoform，或者RNA的可變剪切，看看外顯子使用差異的話，你就需要TopHat, HISAT2或者是STAR這類工具用于找到剪切位點。因為RNA-Seq不同于DNA-Seq，DNA在轉錄成mRNA的時候會把内含子部分去掉。是以mRNA反轉的cDNA如果比對不到參考序列，會被分開，重新比對一次，判斷中間是否有内含子。

工具抉擇

在2016年的一篇綜述A survey of best practices for RNA-seq data analysis，提到目前有三種RNA資料分析的政策。那個時候的工具也主要用的是TopHat,STAR和Bowtie.其中TopHat目前已經被它的作者推薦改用HISAT進行替代。

最近的Nature Communication發表了一篇題為的Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis的文章--被稱之為史上最全RNA-Seq資料分析流程，也是我一直以來想做的事情，隻不過他們做的超乎我的想象。文章中在基于參考基因組的轉錄本分析中所用的工具，是TopHat,HISAT2和STAR，結論就是HISAT2找到junction正确率最高，但是在總數上卻比TopHat和STAR少。從這裡可以看出HISAT2的二類錯誤(納僞）比較少，但是一類錯誤（棄真）就高起來。

就唯一比對而言，STAR是三者最佳的，主要是因為它不會像TopHat和HISAT2一樣在PE比對不上的情況還強行把SE也比對到基因組上。而且在處理較長的read和較短read的不同情況，STAR的穩定性也是最佳的。

就速度而言，HISAT2比STAR和TopHat2平均快上2.5~100倍。

如果學習RNA-Seq資料分析，上面提到的兩篇文獻是必須要看上3遍以上的，而且建議每隔一段時間回顧一下。但是如果就比對工具而言，基本上就是HISAT2和STAR選一個就行。

下載下傳index

首先，問自己一個問題，為什麼比對的時候需要用到index？這裡強烈建議大家去看Jimmy寫的

bowtie算法原理探究bowtie算法原理探究

。但是隻是建議，你不需要真的去看，反正你也看不懂。

高通量測序遇到的第一個問題就是，成千上萬甚至上幾億條read如果在合理的時間内比對到參考基因組上，并且保證錯誤率在接受範圍内。為了提高比對速度，就需要根據參考基因組序列，經過BWT算法轉換成index，而我們比對的序列其實是index的一個子集。當然轉錄組比對還要考慮到可變剪切的情況，是以更加複雜。

是以我門不是直接把read回貼到基因組上，而是把read和index進行比較。人類的index一般都是有現成的，我建議大家下載下傳現成的，我曾經嘗試過用伺服器自己建立index，花的時間讓我懷疑人生。

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

覺得電腦組態還行的，或者是沒有現成index的，可以通過HISAT2的方法進行建立

# 其實hisat2-buld在運作的時候也會自己尋找exons和splice_sites，但是先做的目的是為了提高運作效率
extract_exons.py gencode.v26lift37.annotation.sorted.gtf > hg19.exons.gtf &
extract_splice_sites.py gencode.v26lift37.annotation.gtf > hg19.splice_sites.gtf &
# 建立index， 必須選項是基因組所在檔案路徑和輸出的字首
hisat2-build --ss hg19.splice_sites.gtf --exon hg19.exons.gtf genome/hg19/hg19.fa hg19
           

我的是i7-7700處理器，記憶體是64G，運作的資源效率如下：

正式比對

hisat2基本用法就是

hisat2 [options]* -x <hisat2-idx> {-1 <m1> -2 <m2> | -U <r> } [-S <hit>]

，基本就是提供index的位置，PE資料或者是SE資料存放位置。然而其他可選參數卻是進階的一大名堂。新手就用預設參數呗。

因為RNA-Seq資料是從 SRR3589957 ~ SRR3589962，6個樣本的PE資料，也就是有6次循環，可以寫腳本，也可以直接在指令行裡運作

如下指令運作所在目錄為

/mnt/f/Data/

，我的參考基因組的index資料存放在

/mnt/f/Data/reference/index/hg19/

，而RNA-seq資料存放在

/mnt/f/Data/RNA-Seq

下。比對結果會存放在

/mnt/f/Data/RNA-Seq/aligned

mkdir -p RNA-Seq/aligned
for i in `seq 57 62`
do
    hisat2 -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
           

&

會把任務丢到背景，是以會同時執行這3個比對程式，如果CPU和記憶體承受不住，去掉

&

一個個來。比對這一步是非常消耗記憶體資源的，這是比對工具要将索引資料放入記憶體引起的。我有64G記憶體，理論上可以同時處理20個PE資料。在我的電腦組態下，大緻花了2個小時同時才完成這一步.

基本參數說明

在資料比對的時候，可以安靜一下讀讀HISAT2的額外選項，主要分為如下幾塊

• 主要參數，一定要填寫的内容
• 輸入選項， 對結果影響不大
• 比對選項，主要是

  --n-ceil

  決定模糊字元的數量
• 得分選項， 當一個read比對到不同部位時，确定那個才是最優的。基于mismatch, soft-cliping, gap得分。
• 可變剪切比對選項， 你要決定exon，intron的長度，GT/AG的得分，還可以提供已知的可變剪切和外顯子gtf檔案，
• 報告選項，确定要找多少的位置
• PE選項， 與gap有關的參數
• 輸出選項，建議加上-t記錄時間，其他就是壓縮格式，不影響比對
• SAM選項， 主要是決定SAM的header應該添加哪些内容
• 性能選項和其他選項不考慮

注: soft clipping 指的是比對的read隻有部分比對到參考序列上，還有部分沒有比對上。也就是一個100bp的read，就比對上前面20 bp或者是後面20bp，或者是後面20bp比對的效果不太好。

是以影響比對結果就是比對選項，得分選項，可變剪切選項和PE選項，在有生之年我應該會寫一片文章介紹這些選項對結果的影響。

HISAT2輸出結果

比對之後會輸出如下結果，解讀一下就是全部資料都是100%的，96.68%的配對資料一次都沒有比對，1.23%的資料比是唯一比對，2.09%是多個比對。然後96.68%一次都沒有比對的資料，如果不按照順序來，有0.05%的比對。之後把剩下的部分用單端資料進行比對的話，95.20%資料沒比對上，3.60%的資料比對一次，1.20%比對超過一次。零零總總的加起來是8%的比對！！！

這個總體比對率讓我開始懷疑人生，怎麼可能呀，我翻了翻輸出記錄，發現有幾個結果的比對率超過90%呀。我思索了片刻，驚醒這個實驗好像是用人類和小鼠都做了一遍。于是又去GEO上查了一下記錄，恍然大悟，差點翻車。

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

同時我反思了一下出錯的原因，我預設這個實驗是KO和非KO各3個重複，其實文章的實驗設計并不是如此，可見了解實驗設計很重要，于是我把資料下載下傳這一部分進行了完善。

mkdir -p RNA-Seq/aligned
for i in `seq 56 58`
do
    hisat2 -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)資料格式是目前高通量測序中存放比對資料的标準格式，當然他可以用于存放未比對的資料。是以，

SAM

的格式說明

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

• view: BAM-SAM/SAM-BAM 轉換和提取部分比對
• sort: 比對排序
• merge: 聚合多個排序比對
• index: 索引排序比對
• faidx: 建立FASTA索引，提取部分序列
• tview : 文本格式檢視序列
• pileup : 産生基于位置的結果和 consensus/indel calling

最常用的三闆斧就是格式轉換，排序，索引。而進階教程就是看文檔提高。

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

注

• -S是最新版samtools為了相容以前版本寫的，是以可以省去
• 0.1.19版本和最新版有比較大差别，請注意版本

Jimmy說樣我們仔細判斷sam排序兩種方式的不同，是以我截取前面100行資料，分别排序然後檢視結果。

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

預設排序是根據染色體位置

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
           

也就說說預設按照染色體位置進行排序，而

-n

參數則是根據read名進行排序。當然還有一個

-t

根據TAG進行排序。

說說samtools view

三闆斧的

view

是一個非常實用的子指令，除了之前的格式轉換以外，還能進行資料提取和提取。

比如說提取1号染色體1234-123456區域的比對read

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

在比如搭配

flag

(0.1.19版本沒有）和

flagstat

，使用

-f

或

-F

參數提取不同比對情況的read。

flag是一種描述read比對情況的标記，一種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
           

也就是說如果我想用samtools篩選恰好配對的read,就需要用0x10

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

我應該會在有生之年寫一篇文章好好介紹samtools。

比對質控(QC)

還是在A survey of best practices for RNA-seq data analysis裡面，提到了人類基因組應該有70%~90%的比對率，并且多比對read（multi-mapping reads）數量要少。另外比對在外顯子和所比對鍊（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環境下
pip install RSeQC
           

一共有如下幾個檔案，根據命名就知道功能是啥了。

• 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

先對bam檔案進行統計分析， 從結果上看是符合70~90的比對率要求。

bam_stat.py -i SRR3589956_sorted.bam
           

基因組覆寫率的QC需要提供bed檔案，可以直接RSeQC的網站下載下傳，或者可以用gtf轉換

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

IGV檢視

載入參考序列，注釋和BAM檔案，随便看看吧。

最後請允許我給自己的小密圈打個廣告，如果你非常支援我的創作，想和我保持聯系的話。