比對軟件很多侵贵，首先大家去收集一下，因為我們是帶大家入門，請統(tǒng)一用hisat2(官網(wǎng)https://ccb.jhu.edu/software/hisat2/index.shtml）叫搁，并且搞懂它的用法甩鳄。直接去hisat2的主頁下載index文件即可逞度，然后把fastq格式的reads比對上去得到sam文件。 接著用samtools把它轉(zhuǎn)為bam文件妙啃，并且排序(注意N和P兩種排序區(qū)別)索引好档泽，載入IGV，再截圖幾個基因看看揖赴！* 順便對bam文件進(jìn)行簡單QC馆匿，參考直播我的基因組系列。
來源于生信技能樹：http://www.biotrainee.com/forum.php?mod=viewthread&tid=1750#lastpost

寫在前面比對問題1：為什么要比對

在比對之前燥滑，我們得了解比對的目的是什么渐北？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這類工具用于找到剪切位點(diǎn)澎嚣。因為RNA-Seq不同于DNA-Seq疏尿，DNA在轉(zhuǎn)錄成mRNA的時候會把內(nèi)含子部分去掉。所以mRNA反轉(zhuǎn)的cDNA如果比對不到參考序列易桃，會被分開褥琐，重新比對一次，判斷中間是否有內(nèi)含子晤郑。
鏈接：http://www.reibang.com/p/681e02e7f9af

寫在前面比對問題2：如何比對敌呈？那些軟件

最近的Nature Communication發(fā)表了一篇題為的Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis的文章--被稱之為史上最全RNA-Seq數(shù)據(jù)分析流程，也是我一直以來想做的事情造寝，只不過他們做的超乎我的想象磕洪。文章中在基于參考基因組的轉(zhuǎn)錄本分析中所用的工具，是TopHat,HISAT2和STAR,結(jié)論就是HISAT2找到j(luò)unction正確率最高诫龙，但是在總數(shù)上卻比TopHat和STAR少析显。從這里可以看出HISAT2的二類錯誤(納偽）比較少，但是一類錯誤（棄真）就高起來签赃。
就唯一比對而言谷异，STAR是三者最佳的分尸，主要是因為它不會像TopHat和HISAT2一樣在PE比對不上的情況還強(qiáng)行把SE也比對到基因組上。而且在處理較長的read和較短read的不同情況歹嘹，STAR的穩(wěn)定性也是最佳的箩绍。
就速度而言，HISAT2比STAR和TopHat2平均快上2.5~100倍荞下。
鏈接：http://www.reibang.com/p/681e02e7f9af

我們比對先采用HISAT2

寫在前面index問題1：為什么要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苟跪，花的時間讓我懷疑人生。
鏈接：http://www.reibang.com/p/681e02e7f9af

來自官網(wǎng)：為了用整個index代表整個基因組蔓涧，HISAT2 用小的index覆蓋了整個基因組件已，每個index覆蓋了56 Kbp的范圍，覆蓋整個人類基因組需要55,000 indexes元暴，這些index結(jié)合其他策略可以快速準(zhǔn)確的比對序列篷扩。

寫在前面index問題2：如何獲得index

1 HISAT2官網(wǎng)下載
人類和小鼠的索引有現(xiàn)成的，HISAT2官網(wǎng)可以直接下載進(jìn)行序列比對茉盏。如下圖所示：選擇hg19和mm10的index鉴未，文章中RNA-Seq測序數(shù)據(jù)，可以包括人類和小鼠的數(shù)據(jù)鸠姨，因此需要小鼠和人類的索引歼狼。

$ cd /mnt/f/rna_seq/data
$ mkdir -p reference/index
$ cd reference/index
$ wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz
$ wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/mm10.tar.gz
# 解壓得到兩個目錄，hg19和mm10
$ tar -zxvf *.tar.gz
# 刪除壓縮包
$ rm -rf *.tar.gz

備注：實際我用迅雷（會員）進(jìn)行下載享怀，每個用時大概1個多小時。

2 自己制作

有時候沒有現(xiàn)成的index趟咆，我們就需要自己用HISAT2重新構(gòu)建索引添瓷；包括外顯子梅屉、剪切位點(diǎn)及SNP索引的建立。
參考網(wǎng)站：http://blog.biochen.com/archives/337

參考代碼

# 其實hisat2-buld在運(yùn)行的時候也會自己尋找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坯汤， 必須選項是基因組所在文件路徑和輸出的前綴
hisat2-build --ss hg19.splice_sites.gtf --exon hg19.exons.gtf genome/hg19/hg19.fa hg19

1 開始比對：用hisat2，得到SAM文件

首先啟動miniconda3環(huán)境

$ source ~/miniconda3/bin/activate

先看下hisat2的用法

(base) kelly@DESKTOP-MRA1M1F:/mnt/f/rna_seq/data$ hisat2 -h
HISAT2 version 2.1.0 by Daehwan Kim (infphilo@gmail.com, www.ccb.jhu.edu/people/infphilo)
Usage:
  hisat2 [options]* -x <ht2-idx> {-1 <m1> -2 <m2> | -U <r> | --sra-acc <SRA accession number>} [-S <sam>]

  <ht2-idx>  Index filename prefix (minus trailing .X.ht2).
  <m1>       Files with #1 mates, paired with files in <m2>.
             Could be gzip'ed (extension: .gz) or bzip2'ed (extension: .bz2).
  <m2>       Files with #2 mates, paired with files in <m1>.
             Could be gzip'ed (extension: .gz) or bzip2'ed (extension: .bz2).
  <r>        Files with unpaired reads.
             Could be gzip'ed (extension: .gz) or bzip2'ed (extension: .bz2).
  <SRA accession number>        Comma-separated list of SRA accession numbers, e.g. --sra-acc SRR353653,SRR353654.
  <sam>      File for SAM output (default: stdout)

  <m1>, <m2>, <r> can be comma-separated lists (no whitespace) and can be
  specified many times.  E.g. '-U file1.fq,file2.fq -U file3.fq'.

主要參數(shù)：

• -x <hisat2-idx>
  參考基因組索引文件的前綴搀愧。
• -1 <m1>
  雙端測序結(jié)果的第一個文件惰聂。若有多組數(shù)據(jù)，使用逗號將文件分隔咱筛。Reads的長度可以不一致搓幌。
• -2 <m2>
  雙端測序結(jié)果的第二個文件。若有多組數(shù)據(jù)迅箩，使用逗號將文件分隔溉愁，并且文件順序要和-1參數(shù)對應(yīng)。Reads的長度可以不一致饲趋。
• -U <r>
  單端數(shù)據(jù)文件拐揭。若有多組數(shù)據(jù)，使用逗號將文件分隔奕塑√梦郏可以和-1、-2參數(shù)同時使用龄砰。Reads的長度可以不一致盟猖。
• –sra-acc <SRA accession number>
  輸入SRA登錄號，比如SRR353653寝贡，SRR353654扒披。多組數(shù)據(jù)之間使用逗號分隔。HISAT將自動下載并識別數(shù)據(jù)類型圃泡，進(jìn)行比對碟案。
• -S <hit>
  指定輸出的SAM文件。

• 我的fastq文件在/mnt/f/rna_seq/data
• 我的index在mnt/f/rna_seq/data/reference/index/hg19
• 比對后得到的bam文件會存放在/mnt/f/rna_seq/aligned

#啟動miniconda3環(huán)境（hisat2所在的環(huán)境）
$ source ~/miniconda3/bin/activate
#進(jìn)入data目錄
$cd /mnt/f/rna_seq/aligned
(base) kelly@DESKTOP-MRA1M1F:/mnt/f/rna_seq/aligned$

# 小鼠和人是分開各自比對自己的index
# 人的比對
$ for ((i=56;i<=58;i++));do hisat2 -t -x /mnt/f/rna_seq/data/reference/index/hg19/genome -1 /mnt/f/rna_seq/data/SRR35899${i}.sra_1.fastq.gz -2 /mnt/f/rna_seq/data/SRR35899${i}.sra_2.fastq.gz -S SRR35899${i}.sam ;done
# 小鼠比對
$ for ((i=59;i<=62;i++));do hisat2 -t -x /mnt/f/rna_seq/data/reference/index/mm10/genome -1 /mnt/f/rna_seq/data/SRR35899${i}.sra_1.fastq.gz -2 /mnt/f/rna_seq/data/SRR35899${i}.sra_2.fastq.gz -S SRR35899${i}.sam; done
#結(jié)果一共得到7個sam文件

• SRR3589956比對結(jié)果顯示如下颇蜡，用時32mins

(base) kelly@DESKTOP-MRA1M1F:/mnt/f/rna_seq/aligned$  for ((i=56;i<=58;i++));do hisat2 -t -x /mnt/f/rna_seq/data/reference/index/hg19/genome -1 /mnt/f/rna_seq/data/SRR35899${i}.sra_1.fastq.gz -2 /mnt/f/rna_seq/data/SRR35899${i}.sra_2.fastq.gz -S SRR35899${i}.sam ;done
Time loading forward index: 00:00:05
Time loading reference: 00:00:01
Multiseed full-index search: 00:32:28
28856780 reads; of these:
  28856780 (100.00%) were paired; of these:
    1838758 (6.37%) aligned concordantly 0 times
    24733251 (85.71%) aligned concordantly exactly 1 time
    2284771 (7.92%) aligned concordantly >1 times
    ----
    1838758 pairs aligned concordantly 0 times; of these:
      90903 (4.94%) aligned discordantly 1 time
    ----
    1747855 pairs aligned 0 times concordantly or discordantly; of these:
      3495710 mates make up the pairs; of these:
        2034757 (58.21%) aligned 0 times
        1221304 (34.94%) aligned exactly 1 time
        239649 (6.86%) aligned >1 times
96.47% overall alignment rate

SRR3589959(小鼠）比對結(jié)果顯示如下价说，用時38mins

(base) kelly@DESKTOP-MRA1M1F:/mnt/f/rna_seq/aligned$ for ((i=59;i<=62;i++));do hisat2 -t -x /mnt/f/rna_seq/data/reference/index/mm10/genome -1 /mnt/f/rna_seq/data/SRR35899${i}.sra_1.fastq.gz -2 /mnt/f/rna_seq/data/SRR35899${i}.sra_2.fastq.gz -S SRR35899${i}.sam; done
Time loading forward index: 00:01:03
Time loading reference: 00:00:07
Multiseed full-index search: 00:37:01
30468155 reads; of these:
  30468155 (100.00%) were paired; of these:
    2722534 (8.94%) aligned concordantly 0 times
    24300948 (79.76%) aligned concordantly exactly 1 time
    3444673 (11.31%) aligned concordantly >1 times
    ----
    2722534 pairs aligned concordantly 0 times; of these:
      156866 (5.76%) aligned discordantly 1 time
    ----
    2565668 pairs aligned 0 times concordantly or discordantly; of these:
      5131336 mates make up the pairs; of these:
        3276533 (63.85%) aligned 0 times
        1334397 (26.00%) aligned exactly 1 time
        520406 (10.14%) aligned >1 times
94.62% overall alignment rate
Time searching: 00:37:09
Overall time: 00:38:13

2 SAM文件轉(zhuǎn)換為BAM文件：SAMtools

為什么要轉(zhuǎn)換格式？為了讓計算機(jī)好處理风秤。

BAM and SAM formats are designed to contain the same information. The SAM format is more human readable, and easier to process by conventional text based processing programs, such as awk, sed, python, cut and so on. The BAM format provides binary versions of most of the same data, and is designed to compress reasonably well.
這個PDF文檔講了SAM和BAM格式的區(qū)別與聯(lián)系鳖目，點(diǎn)擊下載
這個pdf更加具體講了兩種格式之間的轉(zhuǎn)換和對照關(guān)系，點(diǎn)此下載

先介紹下SAMtools

SAMtools的wiki介紹：SAMtools Wiki
官方手冊：Manual page from samtools-1.9
必看：詳細(xì)了解SAMtools的用法和來龍去脈
以下引用hoptop

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)階教程就是看文檔提高

今天我們主要是將sam文件轉(zhuǎn)換為bam文件，并對bam文件進(jìn)行sorted（其中有兩種排序方式N和P）缔赠，最后建立索引衍锚。

• 第一種方式

# 首先將比對后的sam文件轉(zhuǎn)換成bam文件
# 利用的是samtools的view選項，參數(shù)-S 輸入sam文件橡淑；參數(shù)-b 指定輸出的文件為bam构拳；最后重定向?qū)懭隻am文件
$ cd mnt/f/rna_seq/aligned
$ for ((i=56;i<=62;i++));do samtools view -S SRR35899${i}.sam -b > SRR35899${i}.bam;done
# 將所有的bam文件按默認(rèn)的染色體位置進(jìn)行排序
$ for ((i=56;i<=62;i++));do samtools sort SRR35899${i}.bam -o SRR35899${i}_sorted.bam;done
# 將所有的排序文件建立索引，索引文件.bai后綴
$ for ((i=56;i<=62;i++));do samtools index SRR35899${i}_sorted.bam;done

速度有點(diǎn)急人

• 第二種方式

for i in `seq 56 62`
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

為什么 BAM 文件 sort 之后體積會變小因為 BAM 文件是壓縮的二進(jìn)制文件梁棠，對文件內(nèi)容排序之后相似的內(nèi)容排在一起置森，使得文件壓縮比提高了，因此排序之后的 BAM 文件變小了符糊，相對應(yīng)的 SAM 文件就是純文本文件凫海，對SAM 文件進(jìn)行排序就不會改變文件大小。而且由于 RNA-seq 中由于基因表達(dá)量的關(guān)系男娄，RNA-seq 的數(shù)據(jù)比對結(jié)果 BAM 文件使用 samtools 進(jìn)行 sort 之后文件壓縮比例變化會比DNA-seq 更甚行贪。另外，samtools 對 BAM 文件進(jìn)行排序之后那些沒有比對上的 reads 會被放在文件的末尾模闲。

3 對BAM進(jìn)行質(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/

(1)對BAM文件進(jìn)行統(tǒng)計分析

RSeQC的安裝，需要先安裝gcc实夹；numpy橄浓；R；Python2.7亮航，這里比較難裝的就是numpy——可以直接利用anaconda安裝

#啟動python2.7環(huán)境
$ source activate python2
#安裝numpy
$ conda install numpy
$ sudo apt-get install python2.7-dev
#安裝gcc
$ sudo apt  install gcc
#查看是否安裝
$ gcc --version
# 用pip命令安裝
$ pip install RSeQC
# 對bam文件進(jìn)行質(zhì)控荸实，其余都同樣的進(jìn)行
$ bam_stat.py  -i SRR3589956_sorted.bam

##############################3加圖

（2）查看基因組覆蓋率

需要提供bed文件，可以直接RSeQC的網(wǎng)站下載缴淋，或者可以用gtf轉(zhuǎn)換

hg19_refseq.png

#下載hg19_RefSeq.bed文件
$ cd /mnt/f/rna_seq/data/reference/genome/hg19
$ wget https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19_RefSeq.bed.gz/download
#查看基因組覆蓋率
$ read_distribution.py -i /mnt/f/rna_seq/aligned/SRR3589956.sorted.bam -r /mnt/f/rna_seq/data/reference/genome/hg19/hg19_RefSeq.bed

4 加載IGV准给，可視化比對結(jié)果

載入?yún)⒖夹蛄行蛊樱⑨尯虰AM文件