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

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

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

比對

比對還是不比對

在比對之前沃呢，我們得了解比對的目的是什么年栓？RNA-Seq數(shù)據(jù)比對和DNA-Seq數(shù)據(jù)比對有什么差異？
RNA-Seq數(shù)據(jù)分析分為很多種薄霜，比如說找差異表達基因或?qū)ふ倚碌目勺兗羟心匙ァＨ绻也町惐磉_基因單純只需要確定不同的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進行替代秘通。

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

下載index

首先跪楞，問自己一個問題，為什么比對的時候需要用到index侣灶？這里強烈建議大家去看Jimmy寫的bowtie算法原理探究bowtie算法原理探究甸祭。但是只是建議，你不需要真的去看褥影，反正你也看不懂淋叶。

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

因此我門不是直接把read回貼到基因組上，而是把read和index進行比較坛缕。人類的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的方法進行創(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ù)卻是進階的一大名堂服猪。新手就用默認參數(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ù)量
• 得分選項惦蚊， 當一個read比對到不同部位時器虾，確定那個才是最優(yōu)的”姆妫基于mismatch, 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ù)進行比對的話铁坎，95.20%數(shù)據(jù)沒比對上蜂奸，3.60%的數(shù)據(jù)比對一次，1.20%比對超過一次硬萍。零零總總的加起來是8%的比對Ｎ涯臁！襟铭！

這個總體比對率讓我開始懷疑人生碌奉，怎么可能呀，我翻了翻輸出記錄寒砖，發(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).

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

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ù)的標準格式，當然他可以用于存放未比對的數(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)換，排序做祝，索引砾省。而進階教程就是看文檔提高。

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行數(shù)據(jù)闰渔，分別排序然后查看結(jié)果嘹承。

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

默認排序是根據(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

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

說說samtools view

三板斧的view是一個非常實用的子命令，除了之前的格式轉(zhuǎ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文件進行統(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)⒖夹蛄屑讼妫⑨尯虰AM文件，隨便看看吧乌询。