转录组分析(7) - 比对

Posted 2023-05-05

参考技术A 对比（Mapping）就是通过局部比对将短的reads定位到参考基因组上，以共后续分析使用。序列比对的过程包括：首先对参考基因组建立索引，建立索引之后，利用比对软件将clean data与参考基因组进行比对，得到sam文件，通过samtools view转换为bam文件，bam文件为二进制文件，对bam文件进行排序，得到排序后的bam文件。

subread是个套件,里面有subread aligner, subjunc aligner, featureCounts, exactSNP。subread aligner可以用于DNA-seq和RNA-seq；当用于RNA-seq时,subread只适用于差异分析；对于检测基因组变异如可变剪接之类的,需要reads的完全比对,这时候可以使用subjunc进行比对；在比对RNA-seq数据时,subread不会取检测exon-exon junctions的存在,只会把exon-spanning eads的最大可比对区域作为比对结果.但是,如果只是进行差异分析的话,subread的结果足以进行.subread的比对上reads可能会比subjunc多.

HISAT2利用大量FM索引以覆盖整个基因组，使用改进的BWT算法，实现了更快的速度和更少的资源占用。

TopHat是一个基于Bowtie的RNA-Seq数据分析工具，它可以快速确认exon-exon剪切拼接事件。TopHat2比对有三个主要的过程，转录组的比对、基因组的比对、剪接比对，双端测序的Read首先每端数据要分别单独进行比对，然后考虑比对的片段长度以及方向将单独比对结果合并在一起形成双端比对结果。

转录组入门：序列比对

任务列表

• 比对软件
• 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文件

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

#! 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 Homo*1.fastq.gz | while read id;do(~/biosoft/hisat2-2.1.0/hisat2 -t -p $NUM_THREADS -x $hg19_ref -1 $data_path/${id%_*}_1.fastq.gz -2 $data_path/${id%_*}_2.fastq.gz 2 > ${id%_*}_map.log | samtools view -Sb  - > ${id%_*}.bam);done
ls --color=never Mus*1.fastq.gz | while read id;do(~/biosoft/hisat2-2.1.0/hisat2 -t -p $NUM_THREADS -x $mm10_ref -1 $data_path/${id%_*}_1.fastq.gz -2 $data_path/${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

运行脚本： 

bash map.sh

质量控制

对bam文件进行简单QC

Reads比对后的质量控制(评估比对质量的指标)

比对上的reads占总reads的百分比；

Reads比对到外显子和参考链上的覆盖度是否一致；

比对到基因组序列，多重比对reads；

相关质控软件除了Picard，RSeQC，Qualimap还有一大堆