转录组入门（5）：序列比对
作业要求
- 比对软件很多，首先大家去收集一下，因为我们是带大家入门，请统一用hisat2，并且搞懂它的用法。
- 直接去hisat2的主页下载index文件即可，然后把fastq格式的reads比对上去得到sam文件。
- 接着用samtools把它转为bam文件，并且排序(注意N和P两种排序区别)索引好，载入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
内容主要参考：http://www.360doc.com/content/16/1223/13/29483982_617058719.shtml
2. HISAT2的使用
人类和小鼠的索引有现成的，HISAT2官网可以直接下载进行序列比对。如下图所示：选择hg19和mm10的index，文章中RNA-Seq测序数据，可以包括人类的3个数据和小鼠的4个数据，因此需要小鼠和人类的索引。
index下载链接
- （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
有时候没有现成的index，我们就需要自己用HISAT2重新构建索引；包括外显子、剪切位点及SNP索引的建立。
参考网站：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>]
参数：
-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
#结果一共得到7个sam文件
比对结束的标准输出
比对结束输出
3.SAMtools 的使用
samtools功能众多，在本次作业中，我们主要学会将sam文件转换为bam文件，并对bam文件进行sorted（其中有两种排序方式N和P），最后建立索引。
# 首先将比对后的sam文件转换成bam文件
# 利用的是samtools的view选项，参数-S 输入sam文件；参数-b 指定输出的文件为bam；最后重定向写入bam文件
$ 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过程中，会出现很多tmp文件，最后会融合所有的tmp文件，形成一个sorted文件。如图所示：
tmp文件
融合所有tmp文件
samtools还有其他功能，包括flagstat命令（查看比对的大致情况）；mileup命令用于生成bcf文件，再使用bcftools进行SNP和Indel的分析；faidx命令对fasta文件建立索引,生成的索引文件以.fai后缀结尾等。
比对大致情况
参考资料：http://blog.csdn.net/sinat_38163598/article/details/72910115
4.比对结果QC
QC的软件有很多，下面3种工具都可以用来质控比对结果（参考了hoptop的简书文章），这里我们使用RSeQC来对我们的比对结果进行质控。
必须吐槽度娘，搜出来的都是什么东西，完全不能看，一脸懵逼。
- 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.jianshu.com/p/14fd4de54402）
# 我的环境已经配置好了，所以直接可以用pip命令安装
$ pip install RSeQC
# 对bam文件进行质控，其余都同样的进行
$ bam_stat.py -i SRR3589956_sorted.bam
RSeQC包括了十多个Python脚本，实现很多功能，具体每个脚本的参数用法，都可以在官网学习这里暂时不多说（之后有时间再详细介绍）。
5.IGV查看比对结果
载入参考基因组，基因组注释文件，很bam文件，看一些基因。
选了一条染色体随便一看
真正的数据分析开始，果然还是很吃力，不过依旧要加油，马上就要学会了。
