4、RNAseq(4)--Hisat2进行序列比对及Samtools格式转化

概况：使用处理后的fastq文件和基因组与转录组比对，确定在转录组或者基因组中的关系。在转录组和基因组的比对采取的方案不同。分别是ungapped alignment to transcriptome和Gapped aligenment to genome。软件：hisat2和STAR在比对上都有比较好的表现。有文献显示，hisat2在纳伪较少但是弃真较多，但是速度比较快。STAR就比对而言综合质量比较好，在长短reads回帖上都有良好发挥。由于hisat2的速度优势，选择hisat2作为本次比对的软件。在比对之前首先要先进行索引文件的获取或者制作。

比对还是不比对

在比对之前，我们得了解比对的目的是什么？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如果比对不到参考序列，会被分开，重新比对一次，判断中间是否有内含子。

比对工具抉择

在2016年的一篇综述A survey of best practices for RNA-seq data analysis，提到目前有三种RNA数据分析的策略。那个时候的工具也主要用的是TopHat,STAR和Bowtie.其中TopHat目前已经被它的作者推荐改用HISAT进行替代。

最近的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倍。

image

如果学习RNA-Seq数据分析，上面提到的两篇文献是必须要看上3遍以上的，而且建议每隔一段时间回顾一下。但是如果就比对工具而言，基本上就是HISAT2和STAR选一个就行。

1. 索引文件的获取

下载index

首先，问自己一个问题，为什么比对的时候需要用到index？这里强烈建议大家去看Jimmy写的bowtie算法原理探究bowtie算法原理探究。但是只是建议，你不需要真的去看，反正你也看不懂。

高通量测序遇到的第一个问题就是，成千上万甚至上几亿条read如果在合理的时间内比对到参考基因组上，并且保证错误率在接受范围内。为了提高比对速度，就需要根据参考基因组序列，经过BWT算法转换成index，而我们比对的序列其实是index的一个子集。当然转录组比对还要考虑到可变剪切的情况，所以更加复杂。

因此我门不是直接把read回贴到基因组上，而是把read和index进行比较。人类的index一般都是有现成的，我建议大家下载现成的，我曾经尝试过用服务器自己创建index，花的时间让我怀疑人生。
所以一般建议的做法就是直接从网站上下载已经index好的数据

hisat2的index

# 下载小鼠基因组的index
mkdir mm10_indexed
cd mm10_indexed/
# 从网站上下载安装包
wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/mm10.tar.gz
#解压
tar zvfx mm10.tar.gz

不同的比对软件构建索引方式不同，所用的索引也不尽相同

索引文件可以去网站下载也可以自己构建。但是索引构建会比较费时间。建立索引文件需要大约一个小时（MAC: 2.6 GHz Intel Core i5/ 8 GB 1600 MHz DDR3）。

网站下载hisat2基因组索引：http://ccb.jhu.edu/software/hisat2/index.shtml

本地索引文件构建参考了CSDN@ Richard_Jolin的构建过程

索引文件的格式如下，是由多个文件构成，要保证索引文件的格式和名称部分一致。

2. hisat2的比对

使用hisat2
公式构建根据hisat2 的使用说明书构建了以下公式：

hisat2  -p 6    -x <dir of index of genome>  -1  seq_val_1.fq.gz   -2  seq_val_2.fq.gz  -S  tem.hisat2.sam

参数说明：
-p #多线程数 -x #参考基因组索引文件目录和前缀 -1 #双端测序中一端测序文件 -2 #同上 -S #输出的sam文件

说明：在比对过程中，hisat会自动将双端测序匹配同一reads并在基因组中比对，最后两个双端测序生成一个sam文件。比对回帖过程需要消耗大量时间和电脑运行速度和硬盘存储空间。5G左右fastq文件比对回帖过程消耗大概一个小时，生成了17G的sam格式文件。回帖完成会生成一个回帖报告。

# 首先利用conda安装版本为2.1.0的hisat软件
conda install hisat2=2.1.0
#因为比对的过程比较长，所以利用nohup命令将程序放在后台执行，并将执行结束的
nohup /home/cenhui2018/QWJ/sequence_data/bio_soft/Hisat2/hisat2-2.1.0/./hisat2 -p 16 -x ~/QWJ/sequence_data/genome/mm10_indexed/mm10/genome  -1 ~/QWJ/sequence_data/20191030_NGS_DATA/Clean_data/19R576_combined_paired_R1.fastq -2 ~/QWJ/sequence_data/20191030_NGS_DATA/Clean_data/19R576_combined_paired_R2.fastq -S ~/QWJ/sequence_data/20191030_NGS_DATA/19R576_paired.hisat2_2.sam > program_2.log 2>&1 &

&会把任务丢到后台，所以会同时执行这3个比对程序，如果CPU和内存承受不住，去掉&一个个来。比对这一步是非常消耗内存资源的，这是比对工具要将索引数据放入内存引起的。

HISAT2输出结果

比对之后会输出如下结果

根据结果，显示的overall的比对率为97.19%

3. SAMtools三板斧

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

最常用的三板斧就是格式转换，排序，索引。而进阶教程就是看文档提高。

# 利用samtools工具将比对得到的sam文件转换为bam文件
samtools view -S 19R576_paired.hisat2_2.sam -b > 19R576_paired.hisat2_2.bam &
# 将得到的sam文件sort一下
samtools sort 19R576_paired.hisat2_2.bam -o 19R576_paired.hisat2_2_sorted.bam
#将得到的sort之后的bam文件排序
samtools index 19R576_paired.hisat2_2_sorted.bam &

4. IGV查看

这个seq_sourted.bam文件可以通过samtools或者IGV( Integrative Genomics Viewer)独立软件进行查看。在IGV软件中载入seq_sourted.bam文件。可以很直观清晰地观察到reads在基因组中的回帖情况和外显子与内含子的关系。

IGV查看结果

参考文献：
https://zhuanlan.zhihu.com/p/61847802
https://www.jianshu.com/p/681e02e7f9af