转录组数据比对hisat2-6

生信菜鸟团

发布于 2024-07-10 16:56:45

1040

发布于 2024-07-10 16:56:45

文章被收录于专栏：生信菜鸟团

生信技能树学习笔记

比对过程：

• 1.建索引

• 2.比对参考基因组

• 3.sam转bam

用到的软件——Hisat2

Hisat2主要是用来进行转录组数据的比对。使用--help查看选项和参数。

hisat2主要参数：

## ----构建索引# 进入参考基因组目录cd $HOME/database/GRCh38.105# Hisat2构建索引，构建索引时间比较长，建议提交后台运行，一般会运行20多分钟左右vim Hisat2Index.shmkdir Hisat2Indexfa=Homo_sapiens.GRCh38.dna.primary_assembly.fafa_baseName=GRCh38.dnahisat2-build -p 12 -f ${fa} Hisat2Index/${fa_baseName}

下一步时间较长，要提前预估服务器空间是否够用：du -sh ./

# 运行nohup sh Hisat2Index.sh >Hisat2Index.sh.log &

## ----比对# 进入比对文件夹cd $HOME/project/Human-16-Asthma-Trans/Mapping/Hisat2

定义三个文件夹：index索引、inputdir输入文件夹、out输出文件夹。这一步需要激活conda环境。

单个样本比对

## 单个样本比对，步骤分解index=/home/t_rna/database/GRCh38.104/Hisat2Index/GRCh38.dnainputdir=$HOME/project/Human-16-Asthma-Trans/data/cleandata/trim_galore/outdir=$HOME/project/Human-16-Asthma-Trans/Mapping/Hisat2#-p 10代表线程，根据自己的情况改hisat2 -p 10 -x ${index} \ -1 ${inputdir}/SRR1039510_1_val_1.fq.gz \ -2 ${inputdir}/SRR1039510_2_val_2.fq.gz \ -S ${outdir}/SRR1039510.Hisat_aln.sam

SRR1039510_1_val_1.fq.gz和SRR1039510_2_val_2.fq.gz是之前过滤后的文件，过滤方式在之前的笔记中。SRR1039510.Hisat_aln.sam是比对结果。

上述操作的结果

Sort指对sam文件进行排序，-@ 15是指使用线程，-o指生成的bam文件名字，后面接转之前的sam文件

# sam转bamsamtools sort -@ 15 -o SRR1039510.Hisat_aln.sorted.bam SRR1039510.Hisat_aln.sam

运行结果

Sam文件bam文件内容一样，主要区别是文件大小。

多个样本比对

这里需要用到管道符|串联比对参考基因组和 sam转bam两个步骤

这里的2代表下面这个程序中输出的过程，并将其重定向到样本对应的log文件中

关注点：

• 总比对率：一般都能在80%以上

• 唯一比对：92.89%，越高越好

具体程序

# 多个样本批量进行比对，排序，建索引# Hisat.sh内容：注意命令中的-，表示占位符，表示|管道符前面的输出。## 此处索引直接使用服务器上已经构建好的进行练习# vim Hisat.shindex=/home/t_rna/database/GRCh38.104/Hisat2Index/GRCh38.dnainputdir=$HOME/project/Human-16-Asthma-Trans/data/cleandata/trim_galore/outdir=$HOME/project/Human-16-Asthma-Trans/Mapping/Hisat2 cat ../../data/cleandata/trim_galore/ID | while read iddohisat2 -p 5 -x ${index} -1 ${inputdir}/${id}_1_val_1.fq.gz -2 ${inputdir}/${id}_2_val_2.fq.gz 2>${id}.log | samtools sort -@ 3 -o ${outdir}/${id}.Hisat_aln.sorted.bam -done # 提交后台运行nohup sh Hisat.sh >Hisat.log & # 统计比对情况multiqc -o ./ SRR*log

结果

可视化结果

比对率过低可能

1.细菌污染 2.核糖体RNA 3.比对文件物种错误

比对结果文件：sam/bam格式

SAM(The Sequence Alignment/Map format)格式，即序列比对文件格式，详细

介绍见：http://samtools.github.io/hts-specs/SAMv1.pdf

BAM是SAM的二进制文件(B源自binary)

sam/bam头部

sam/bam主体区

比对结果部分（alignment section）

1.每一行表示一个read的比对信息。

2.每行包括11个必须字段和1个可选字段，字段之间用制表符分割。

sam/bam格式-FLAG

FLAG详解：http://broadinstitute.github.io/picard/explain-flags.html

FLAG: Combination of bitwise FLAGs，表示比对的详细情况

例如：我要看FLAG 99是什么意思：samtools flags 99

0x63 99 PAIRED,PROPER_PAIR,MREVERSE,READ1

0x63 只能由0x1，0x2，0x20和0x40组成，也就是这个read是双端测序；这个序列和参考序列完全匹配，没有插入缺失；这个序列对应的另一端序列比对到参考序列的负链上；代表这个序列是R1端序列， read1

sam/bam格式-CIGAR

CIGAR详解

CIGAR string，简要比对信息表达式（Compact Idiosyncratic Gapped Alignment Report），其以参考序列为基础，使用数字加字母表示比对结果，比如3S6M1P1I4M，前三个碱基被剪切去除了，然后6个比对上了，然后打开了一个缺口，有一个碱基插入，最后是4个比对上了，是按照顺序的，字母的含义如下