转录组数据分析笔记（2）——使用Hisat2比对参考基因组

在进行clean data与参考基因组比对之前，我们需要先建立参考基因组索引。进入下载好参考基因组的文件目录下，运行命令：

$ hisat2-build Arabidopsis_thaliana.dna.genome.fa genome

就可以在当前目录建立参考基因组索引文件，hisat2固定会生成8个以.ht2做后缀名的索引文件，如下所示：

需要注意的一点，比对软件除了hiasat2以外，还有subjunk、bwa、bowtie2等等，各个比对软件生成的索引文件是不同的，不能相互混用，命名的时候注意区分各种比对工具。

比对的意思是将每一个read与参考基因组序列进行对比，目的是得到每一个read在参考基因组上的位置信息，有了这个基础的位置信息才可以进行后续基因或者转录本的定量，最终由定量结果做差异表达矩阵，分析上调或者下调的基因数量。

新建一个shell脚本输入下面的代码

1
2
3
4
5
6
7
8

#!/bin/bash
for i in `seq 194 209`
do
    hisat2 -p 4 -x /media/sf_example/data/ref/genome \  #索引文件绝对路径
    -1 /media/sf_example/data/clean_data/ERR1698"$i"_1.fq.gz \
    -2 /media/sf_example/data/clean_data/ERR1698"$i"_2.fq.gz \
    -S /media/sf_example/data/hisat2_sam/ERR1698"$i".sam    #注意大写的S
done

参数解释：

输出到屏幕的结果如下，我们选取其中一个进行解读：

共有6166091对测序数据，都是双侧测序数据，其中：

最后一块是对两条链拆开比对的结果，这个一般用不到，本来测序的两条reads就应该比对到同一个染色体同一个基因附近，拆开比对到不同染色体没有意义。我们要看的是最后一句话，总比对率为99.97%，通常比对率大于90%说明比对情况较好，与参考基因组基本吻合。

比对结果除了有屏幕上输出的总体报告外，还有记录详细比对结果的sam文件。双端测序的比对会将两个测序文件进行整合和比较，最后只生成一个sam文件，因此这个sam文件非常大，hisat2比对生成的sam文件可以直接打开。我们可以选取一部分进行解读。

VN是格式版本；SO表示比对排序的类型，有unknown，unsorted，queryname和coordinate几种。samtools软件在进行行排序后不能自动更新sam文件的SO值。

参考序列名，这些参考序列决定了比对结果sort的顺序，SN是参考序列名；LN是参考序列长度；每个参考序列为一行。这里表示拟南芥有5条染色体，对应长度都在后面，Mt是线粒体基因，Pt是叶绿体基因。

这里包括了路径，方法，以及我质控后的序列长度（50-100）等详细信息。

后面提供额外的信息，一般不重要，了解一下就行。因为sam文件太大（往往有10G以上），也不适合电脑进行后续处理，所以我们会用到samtools，将sam文件转化为更适合电脑处理的二进制bam文件。这个后面会讲。

以下4种软件均用于序列比对，用法稍有不同，做个记录

$ hisat2 -p 4 -x 索引目录 -1 单端测序数据文件 -2 另一端测序数据文件 -S 输出文件

$ subjunk -T 4 -i 索引目录 -r 单端测序数据文件 -R 另一端测序数据文件 -o 输出文件

$ bowtie2 -p 4 -x 索引目录 -1 单端测序数据文件 -2 另一端测序数据文件 -S 输出文件

$ bwa mem -t 4 -M 索引目录 单端测序数据文件 另一端测序数据文件 > 输出文件