GATK4 实用技巧
前言
本篇笔记围绕GATK4流程中注意事项进行介绍,包括如何选择合适的线程和内存?如何按步骤进行数据预处理。
- GATK实战工作流程
- 脚本优化方法
- 参数设置技巧
GATK 全称“基因组分析工具箱” Genome Analysis Toolkit,它是一组命令行工具,用于分析高通量测序数据,主要侧重于变异发现。这些工具可以单独使用,也可以链接在一起形成完整的工作流程。
所需数据与软件
需要有简单的Linux知识储备和高通量基因组学基础,流程中用到的软件如下:
GATK 4.3.0.0 fastp 0.20.1 bwa 0.7.17 samtools 1.11 mosdepth 0.3.0
示例数据来自ERP001960项目,数据结构为illumina测序原始数据,包含三个样本,如果想重复测试需要大概700GB空间。
数据预处理
预处理的步骤包括质控、去接头、排序等,需要用到FATSP等工具,为了提高效率建议直接用fastp标准输出到BAM文件,避免不必要的磁盘读写,用管道串联起来能提高工作效率。
无论是否进行排序,都不要直接将未压缩的sam文件写入磁盘!最好用管道直接把两步合为一起生成排序去重之后的bam文件
fastp去除接头
使用fastp程序进行操作,为了提高数据吞吐效率,将默认输出信息定向到/dev/null中,避免磁盘读写。
fastp -i /lscratch/${SLURM_JOB_ID}/$(basename $r1) \ -I /lscratch/${SLURM_JOB_ID}/$(basename $r2) \ --stdout --thread $threads \ -j "${logs}/fastp-${SLURM_JOB_ID}.json" \ -h "${logs}/fastp-${SLURM_JOB_ID}.html" \ > /dev/null
下图展示了不同线程与处理速度的关系,很显然啊!这里如果超过6个线程就没多大意义了,设置6个线程就能达到最佳计算效率。
bwa mem 序列比对
比对过程实际上是拿着参考基因组,然后把每个小片段进行对齐,举个比较抽象的例子,就好比是拼图的过程,目的是让每一个小碎片都找到匹配的位置。
使用bwa mem程序进行比对,注意这里每个样品需要两条序列(因为二代双端测序两个方向是一对儿)。
bwa mem -M -t ${threads} \ -R "@RG\tID:$id\tPL:ILLUMINA\tLB:$lb\tSM:$sm" \ $genome \ /lscratch/$SLURM_JOBID/$(basename ${r1[small]}) \ /lscratch/$SLURM_JOBID/$(basename ${r2[small]}) \ > /dev/null
通过下图可以看出线程数与计算效率的关系,很明显在32线程以内基本呈正相关,因此建议此处设置的线程数稍微高一点,理论上线程数越多的话,拼图的速度也越快,就好比你找了几个小伙伴帮你同时拼图。
但是,也需要注意不要设置的太高,否侧多个线程的并行过程可能出现干扰,反而影响速度。就好比你为了拼一张地图把全校的人都叫过来了,一间屋里全是人,每个人移动起来都费劲,CPU也一样,线程过多相互打架。
samtools sort 排序
SAM文件输出后直接管道传递给sort进行排序,这样可以避免磁盘读写浪费时间(就好比买了一份外卖我直接吃就行了,没必要拿到厨房加工一下再取出来,多此一举浪费时间)。可以再理解一下磁盘吞吐和溢出概念,这对提高性能非常重要。
使用重定向把输出信息丢弃,基准测试过程中BAM文件大概25GB左右,临时文件储存在lscratch中,需要格外注意设置线程数和每个线程的内存大小时综合考虑,需要根据机器配置和计算需求灵活修改。
samtools sort -T /lscratch/$SLURM_JOBID/part \ -m ${mem_mb_per_thread}M \ -@ ${threads} \ -O BAM \ /lscratch/$SLURM_JOBID/$(basename $bam) \ > /dev/null
通过对比可以发现,在设置10个线程时,计算效率最高,盲目的增加线程数反而导致计算效率降低了。另外,启用10个线程之后,设置内存参数也对速度有影响,例如16线程状态下,采用10GB内存的速度比30GB明显要快。为了方便理解,举个不太抽象的例子:假如有很多乒乓球,每个上面有一个数字编号,打乱后放在一个屋子里,想让你按顺序排好放到另一个屋子。
如果人数增加很明显是有助于提高排序速度(线程数增加),正常情况下每个人一次只能拿一个球(单个线程内存比较小),如果让每个人每次可以多拿几个球就可以避免来回跑(增加内存也能提高效率)。
但是存在某个临界情况,此时整体效率最大,也就是线程数和内存最匹配(上面图中的最高点处)。
samtools view转换bam文件
上步生成的sam文件要转换成bam文件,使用如下命令,此步骤占用的内存相对较少。
samtools view -@ ${threads} \ -O BAM \ /lscratch/\$SLURM_JOBID/$(basename $sam) \ > /dev/null
通过下图发现这一步设置的线程数和计算速度呈正比关系,测试最大采用16线程。
脚本的优化过程
有没有一种更厉害的方法,能达到比较高的计算效率呢?
下述代码中,使用两个管道命令将中间结果依次传递,采用32核心进行计算,完成了fastp、bwa、samtools等步骤。
fastp -i $read1 -I $read2 \ --stdout --thread 2 \ -j "${logs}/fastp-${SLURM_JOB_ID}.json" \ -h "${logs}/fastp-${SLURM_JOB_ID}.html" \ 2> "${logs}/fastp-${SLURM_JOB_ID}.log" \ | bwa mem -v 2 -M -t 32 -p \ -R "@RG\tID:$id\tPL:ILLUMINA\tLB:$lb\tSM:$sm" \ $genome - 2> ${logs}/bwa-${SLURM_JOB_ID}.log \ | samtools sort -T /lscratch/$SLURM_JOBID/part \ -m 1706M \ -@ 12 \ -O BAM \ -o /lscratch/$SLURM_JOBID/$bam \ 2> ${logs}/samtools-${SLURM_JOB_ID}.log
运行中占用的最大线程数为45左右,最大内存为60GB左右,大部分时间都是32线程进行运算,比对结果逐步记录到缓存区,当缓存区满了后再进行排序,这样极大地提高了效率。
https://hpc.nih.gov/training/gatk_tutorial
END
© 素材来源于网络,侵权请联系后台删除
往期推荐: