序列比对 hisat2
转录组入门(5)：序列比对
任务列表
比对软件
hisat2的用法
下载index文件
比对、排序、索引
质量控制
载入IGV，截图几个基因
hisat2的用法
本作业是比对到基因组，所以使用gapped or splices mapper，此流程已经更新。TopHat首次被发表已经是7年前，STAR的比对速度是TopHat的50倍，HISAT更是STAR的1.2倍。HISAT2是TopHat2/Bowti2的继任者，使用改进的BWT算法，实现了更快的速度和更少的资源占用，作者推荐TopHat2/Bowti2和HISAT的用户转换到HISAT2。
官网：https://ccb.jhu.edu/software/hisat2/index.shtml（学习一个软件最好的方法就是结合现有中文资料，加上阅读官方说明书和HELP文档，一般刚开始学习的时候先使用默认参数，不要乱调参数）
下载index文件
1
2
3
4
5
6
cd ~/reference
mkdir -p index/hisat && cd index/hisat
wget -c ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz
wget -c ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/mm10.tar.gz
tar zxvf hg19.tar.gz
tar xvzf mm10.tar.gz
-c：断点续传
比对、排序、索引
把fastq格式的reads比对上去得到sam文件，接着用samtools把它转为bam文件，并且排序（注意N和P两种排序区别）索引好（可以使用管道实现，省去中间SAM保存的过程，直接输出BAM文件）
编写bash脚本：map.sh
1
2
3
4
5
6
7
8
9
10
11
12
#! usr/bin/bash
set -u
set -e
set -o pipefail
hg19_ref=/mnt/hgfs/2017/reference/index/hisat/hg19/genome
mm10_ref=/mnt/hgfs/2017/reference/index/hisat/mm10/genome
data_path=/mnt/hgfs/2017/rna_seq/data
NUM_THREADS=25
ls --color=never Homo1.fastq.gz | while read id;do(~/biosoft/hisat2-2.1.0/hisat2 -t -p $NUM_THREADS -x $hg19_ref -1
d
a
t
a
p
a
t
h
/
data_path/
datapath/{id%_}1.fastq.gz -2
d
a
t
a
p
a
t
h
/
data_path/
datapath/{id%}2.fastq.gz 2 > ${id%}map.log | samtools view -Sb - > ${id%}.bam);done
ls --color=never Mus1.fastq.gz | while read id;do(~/biosoft/hisat2-2.1.0/hisat2 -t -p $NUM_THREADS -x $mm10_ref -1
d
a
t
a
p
a
t
h
/
data_path/
datapath/{id%_}1.fastq.gz -2
d
a
t
a
p
a
t
h
/
data_path/
datapath/{id%}2.fastq.gz 2 > ${id%}map.log | samtools view -Sb - > ${id%}.bam);done
ls --color=never .bam | while read id;do(samtools sort --threads $NUM_THREADS $id -o ${id%.}_sorted.bam);done
ls --color=never *_sorted.bam | while read id;do(samtools index $id);done
运行脚本：
1
bash map.sh
质量控制
对bam文件进行简单QC
Reads比对后的质量控制(评估比对质量的指标)
比对上的reads占总reads的百分比；
Reads比对到外显子和参考链上的覆盖度是否一致；
比对到基因组序列，多重比对reads；
相关质控软件除了Picard，RSeQC，Qualimap还有一大堆