RNA-seq(5):序列比对：Hisat2
比对软件很多，首先大家去收集一下，因为我们是带大家入门，请统一用hisat2(官网https://ccb.jhu.edu/software/hisat2/index.shtml），并且搞懂它的用法。直接去hisat2的主页下载index文件即可，然后把fastq格式的reads比对上去得到sam文件。 接着用samtools把它转为bam文件，并且排序(注意N和P两种排序区别)索引好，载入IGV，再截图几个基因看看！* 顺便对bam文件进行简单QC，参考直播我的基因组系列。 来源于生信技能树：http://www.biotrainee.com/forum.php?mod=viewthread&tid=1750#lastpost
写在前面比对问题1：为什么要比对
在比对之前，我们得了解比对的目的是什么？RNA-Seq数据比对和DNA-Seq数据比对有什么差异？ RNA-Seq数据分析分为很多种，比如说找差异表达基因或寻找新的可变剪切。如果找差异表达基因单纯只需要确定不同的read计数就行的话，我们可以用bowtie, bwa这类比对工具，或者是salmon这类align-free工具，并且后者的速度更快。 但是如果你需要找到新的isoform，或者RNA的可变剪切，看看外显子使用差异的话，你就需要TopHat, HISAT2或者是STAR这类工具用于找到剪切位点。因为RNA-Seq不同于DNA-Seq，DNA在转录成mRNA的时候会把内含子部分去掉。所以mRNA反转的cDNA如果比对不到参考序列，会被分开，重新比对一次，判断中间是否有内含子。 链接：https://www.jianshu.com/p/681e02e7f9af
写在前面比对问题2：如何比对？那些软件
最近的Nature Communication发表了一篇题为的Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis的文章--被称之为史上最全RNA-Seq数据分析流程，也是我一直以来想做的事情，只不过他们做的超乎我的想象。文章中在基于参考基因组的转录本分析中所用的工具，是TopHat,HISAT2和STAR,结论就是HISAT2找到junction正确率最高，但是在总数上却比TopHat和STAR少。从这里可以看出HISAT2的二类错误(纳伪）比较少，但是一类错误（弃真）就高起来。 就唯一比对而言，STAR是三者最佳的，主要是因为它不会像TopHat和HISAT2一样在PE比对不上的情况还强行把SE也比对到基因组上。而且在处理较长的read和较短read的不同情况，STAR的稳定性也是最佳的。 就速度而言，HISAT2比STAR和TopHat2平均快上2.5~100倍。 链接：https://www.jianshu.com/p/681e02e7f9af
我们比对先采用HISAT2
写在前面index问题1：为什么要index
为什么比对的时候需要用到index？这里强烈建议大家去看Jimmy写的bowtie算法原理探究bowtie算法原理探究。 高通量测序遇到的第一个问题就是，成千上万甚至上几亿条read如果在合理的时间内比对到参考基因组上，并且保证错误率在接受范围内。为了提高比对速度，就需要根据参考基因组序列，经过BWT算法转换成index，而我们比对的序列其实是index的一个子集。当然转录组比对还要考虑到可变剪切的情况，所以更加复杂。 因此我门不是直接把read回贴到基因组上，而是把read和index进行比较。人类的index一般都是有现成的，我建议大家下载现成的，我曾经尝试过用服务器自己创建index，花的时间让我怀疑人生。 链接：https://www.jianshu.com/p/681e02e7f9af
来自官网：为了用整个index代表整个基因组，HISAT2 用小的index覆盖了整个基因组，每个index覆盖了56 Kbp的范围，覆盖整个人类基因组需要55,000 indexes，这些index结合其他策略可以快速准确的比对序列。
写在前面index问题2：如何获得index
1 HISAT2官网下载 人类和小鼠的索引有现成的，HISAT2官网可以直接下载进行序列比对。如下图所示：选择hg19和mm10的index，文章中RNA-Seq测序数据，可以包括人类和小鼠的数据，因此需要小鼠和人类的索引。
$ 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备注：实际我用迅雷（会员）进行下载，每个用时大概1个多小时。
2 自己制作
有时候没有现成的index，我们就需要自己用HISAT2重新构建索引；包括外显子、剪切位点及SNP索引的建立。 参考网站：http://blog.biochen.com/archives/337
# 其实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 hg191 开始比对：用hisat2，得到SAM文件
首先启动miniconda3环境
$ 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'.主要参数：
- -x <hisat2-idx> 参考基因组索引文件的前缀。
- -1 <m1> 双端测序结果的第一个文件。若有多组数据，使用逗号将文件分隔。Reads的长度可以不一致。
- -2 <m2> 双端测序结果的第二个文件。若有多组数据，使用逗号将文件分隔，并且文件顺序要和-1参数对应。Reads的长度可以不一致。
- -U <r> 单端数据文件。若有多组数据，使用逗号将文件分隔。可以和-1、-2参数同时使用。Reads的长度可以不一致。
- –sra-acc <SRA accession number> 输入SRA登录号，比如SRR353653，SRR353654。多组数据之间使用逗号分隔。HISAT将自动下载并识别数据类型，进行比对。
- -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环境（hisat2所在的环境）
$ source ~/miniconda3/bin/activate
#进入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
#结果一共得到7个sam文件- SRR3589956比对结果显示如下，用时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 rateSRR3589959(小鼠）比对结果显示如下，用时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:132 SAM文件转换为BAM文件：SAMtools
为什么要转换格式？为了让计算机好处理。
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格式的区别与联系，点击下载 这个pdf更加具体讲了两种格式之间的转换和对照关系，点此下载
先介绍下SAMtools
SAMtools的wiki介绍：SAMtools Wiki 官方手册：Manual page from samtools-1.9 必看：详细了解SAMtools的用法和来龙去脉 以下引用hoptop
SAM（sequence Alignment/mapping)数据格式是目前高通量测序中存放比对数据的标准格式，当然他可以用于存放未比对的数据。所以，SAM的格式说明 而目前处理SAM格式的工具主要是SAMTools，这是Heng Li大神写的.除了C语言版本，还有Java的Picard，Python的Pysam，Common lisp的cl-sam等其他版本。SAMTools的主要功能如下：
- view: BAM-SAM/SAM-BAM 转换和提取部分比对
- sort: 比对排序
- merge: 聚合多个排序比对
- index: 索引排序比对
- faidx: 建立FASTA索引，提取部分序列
- tview: 文本格式查看序列
- pileup: 产生基于位置的结果和 consensus/indel calling 最常用的是格式转换，排序，索引。而进阶教程就是看文档提高
今天我们主要是将sam文件转换为bam文件，并对bam文件进行sorted（其中有两种排序方式N和P），最后建立索引。
- 第一种方式
# 首先将比对后的sam文件转换成bam文件
# 利用的是samtools的view选项，参数-S 输入sam文件；参数-b 指定输出的文件为bam；最后重定向写入bam文件
$ 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文件按默认的染色体位置进行排序
$ 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速度有点急人
- 第二种方式
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 文件是压缩的二进制文件，对文件内容排序之后相似的内容排在一起，使得文件压缩比提高了，因此排序之后的 BAM 文件变小了，相对应的 SAM 文件就是纯文本文件，对SAM 文件进行排序就不会改变文件大小。而且由于 RNA-seq 中由于基因表达量的关系，RNA-seq 的数据比对结果 BAM 文件使用 samtools 进行 sort 之后文件压缩比例变化会比DNA-seq 更甚。另外，samtools 对 BAM 文件进行排序之后那些没有比对上的 reads 会被放在文件的末尾。
3 对BAM进行质控QC
在A survey of best practices for RNA-seq data analysis里面，提到了人类基因组应该有70%~90%的比对率，并且多比对read（multi-mapping reads）数量要少。另外比对在外显子和所比对链（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文件进行统计分析
RSeQC的安装，需要先安装gcc；numpy；R；Python2.7，这里比较难装的就是numpy——可以直接利用anaconda安装
#启动python2.7环境
$ 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文件进行质控，其余都同样的进行
$ bam_stat.py -i SRR3589956_sorted.bam##############################3加图
（2）查看基因组覆盖率
需要提供bed文件，可以直接RSeQC的网站下载，或者可以用gtf转换
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.bed4 加载IGV，可视化比对结果
载入参考序列，注释和BAM文件