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

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

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

比對

比對還是不比對

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

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

工具抉擇

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

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

下載index

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

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

因此我門不是直接把read回貼到基因組上，而是把read和index進(jìn)行比較。人類的index一般都是有現(xiàn)成的，我建議大家下載現(xiàn)成的，我曾經(jīng)嘗試過用服務(wù)器自己創(chuàng)建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

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

# 其實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處理器，內(nèi)存是64G，運行的資源效率如下：

正式比對

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ù)唄。

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

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

基本參數(shù)說明

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

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

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

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

HISAT2輸出結(jié)果

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

這個總體比對率讓我開始懷疑人生，怎么可能呀，我翻了翻輸出記錄，發(fā)現(xiàn)有幾個結(jié)果的比對率超過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).

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

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)數(shù)據(jù)格式是目前高通量測序中存放比對數(shù)據(jù)的標(biāo)準(zhǔn)格式，當(dāng)然他可以用于存放未比對的數(shù)據(jù)。所以，SAM的格式說明

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

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

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

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說樣我們仔細(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

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

說說samtools view

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

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

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

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

比對質(zhì)控(QC)

還是在A survey of best practices for RNA-seq data analysis里面，提到了人類基因組應(yīng)該有70%~90%的比對率，并且多比對read（multi-mapping reads）數(shù)量要少。另外比對在外顯子和所比對鏈（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

一共有如下幾個文件，根據(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

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

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)⒖夹蛄?，注釋和BAM文件，隨便看看吧。