作業(yè)要求
- 比對軟件很多，首先大家去收集一下，因為我們是帶大家入門，請統(tǒng)一用hisat2，并且搞懂它的用法。
- 直接去hisat2的主頁下載index文件即可，然后把fastq格式的reads比對上去得到sam文件。
- 接著用samtools把它轉(zhuǎn)為bam文件，并且排序(注意N和P兩種排序區(qū)別)索引好，載入IGV，再截圖幾個基因看看！
- 順便對bam文件進行簡單QC，參考直播我的基因組系列。
來源于生信技能樹：http://www.biotrainee.com/forum.php?mod=viewthread&tid=1750#lastpost
實驗過程
1. 比對軟件
- HISAT2：http://ccb.jhu.edu/software/hisat2/index.shtml
參考資料：http://blog.biochen.com/archives/337 - STAR：https://codeload.github.com/alexdobin/STAR/zip/master
參考資料：http://www.bio-info-trainee.com/727.html - TopHat：http://ccb.jhu.edu/software/tophat/index.shtml
參考資料：http://blog.sina.com.cn/s/blog_8808cae20101amqp.html - RapMap：https://github.com/COMBINE-lab/RapMap
參考文獻：https://academic.oup.com/bioinformatics/article/32/12/i192/2288985/RapMap-a-rapid-sensitive-and-accurate-tool-for - CIDANE：http://ccb.jhu.edu/software/cidane/
參考文獻：https://genomebiology.biomedcentral.com/articles/10.1186/s13059-015-0865-0 - CLASS2 ：https://sourceforge.net/projects/splicebox/files/?source=navbar
參考文獻：https://academic.oup.com/nar/article/44/10/e98/2516329/CLASS2-accurate-and-efficient-splice-variant
內(nèi)容主要參考：http://www.360doc.com/content/16/1223/13/29483982_617058719.shtml
2. HISAT2的使用
人類和小鼠的索引有現(xiàn)成的，HISAT2官網(wǎng)可以直接下載進行序列比對。如下圖所示：選擇hg19和mm10的index，文章中RNA-Seq測序數(shù)據(jù)，可以包括人類的3個數(shù)據(jù)和小鼠的4個數(shù)據(jù)，因此需要小鼠和人類的索引。
- （1）人類和小鼠index下載
$ mkdir -p ~/disk2/data/reference/index
$ cd ~/disk2/data/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
有時候沒有現(xiàn)成的index，我們就需要自己用HISAT2重新構(gòu)建索引；包括外顯子、剪切位點及SNP索引的建立。
參考網(wǎng)站：http://blog.biochen.com/archives/337
- （2）比對序列，得到sam文件
Usage：hisat2 [options]* -x <hisat2-idx> {-1 <m1> -2 <m2> | -U <r> | --sra-acc <SRA accession number>} [-S <hit>]
參數(shù)：
-x 指定index文件
-1 <m1> 雙端測序第一個文件
-2 <m2> 雙端測序第二個文件
-U 單端測序文件
--sra-acc SRA accession number
-S 指定輸出的格式，一般指定為sam
# 小鼠和人是分開各自比對自己的index
# 人的比對
$ for ((i=56;i<=58;i++));do hisat2 -t -x ~/disk2/data/reference/index/hg19/genome -1 ~/disk2/data/rna-seq/SRR35899${i}_1.fastq.gz -2 ~/disk2/data/rna-seq/SRR35899${i}_2.fastq.gz -S SRR35899${i}.sam;done
# 小鼠比對
$ for ((i=59;i<=62;i++));do hisat2 -t -x ~/disk2/data/reference/index/mm10/genome -1 ~/disk2/data/rna-seq/SRR35899${i}_1.fastq.gz -2 ~/disk2/data/rna-seq/SRR35899${i}_2.fastq.gz -S SRR35899${i}.sam;done
#結(jié)果一共得到7個sam文件
比對結(jié)束的標準輸出
3.SAMtools 的使用
samtools功能眾多，在本次作業(yè)中，我們主要學會將sam文件轉(zhuǎn)換為bam文件，并對bam文件進行sorted（其中有兩種排序方式N和P），最后建立索引。
# 首先將比對后的sam文件轉(zhuǎn)換成bam文件
# 利用的是samtools的view選項，參數(shù)-S 輸入sam文件；參數(shù)-b 指定輸出的文件為bam；最后重定向?qū)懭隻am文件
$ cd ~/disk2/data/reference/rna-seq/aligned
$ for ((i=56;i<=62;i++));do samtools view -S SRR35899${i}.sam -b > SRR35899${i}.bam;done
# 將所有的bam文件進行排序
$ 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
sort過程中，會出現(xiàn)很多tmp文件，最后會融合所有的tmp文件，形成一個sorted文件。如圖所示：
tmp文件
融合所有tmp文件
samtools還有其他功能，包括flagstat命令（查看比對的大致情況）；mileup命令用于生成bcf文件，再使用bcftools進行SNP和Indel的分析；faidx命令對fasta文件建立索引,生成的索引文件以.fai后綴結(jié)尾等。
參考資料：http://blog.csdn.net/sinat_38163598/article/details/72910115
4.比對結(jié)果QC
QC的軟件有很多，下面3種工具都可以用來質(zhì)控比對結(jié)果（參考了hoptop的簡書文章），這里我們使用RSeQC來對我們的比對結(jié)果進行質(zhì)控。
必須吐槽度娘，搜出來的都是什么東西，完全不能看，一臉懵逼。
- RSeQC——http://rseqc.sourceforge.net/
- Qualimap——http://qualimap.bioinfo.cipf.es/
- Picard——http://broadinstitute.github.io/picard/
# RSeQC的安裝，需要先安裝gcc；numpy；R；Python2.7，這里比較難裝的就是numpy——可以直接利用anaconda安裝（http://www.itdecent.cn/p/14fd4de54402）
# 我的環(huán)境已經(jīng)配置好了，所以直接可以用pip命令安裝
$ pip install RSeQC
# 對bam文件進行質(zhì)控，其余都同樣的進行
$ bam_stat.py -i SRR3589956_sorted.bam
RSeQC包括了十多個Python腳本，實現(xiàn)很多功能，具體每個腳本的參數(shù)用法，都可以在官網(wǎng)學習這里暫時不多說（之后有時間再詳細介紹）。
5.IGV查看比對結(jié)果
載入?yún)⒖蓟蚪M，基因組注釋文件，很bam文件，看一些基因。
真正的數(shù)據(jù)分析開始，果然還是很吃力，不過依舊要加油，馬上就要學會了。
