基于SLURM HPC批量运行GATK4 HaplotypeCaller的技术方案咨询

背景 #

已完成黑麦草全基因组测序（PE 150bp）数据预处理：通过trim galore去除低质量reads，bwa比对至参考基因组，完成BAM排序、添加RG标签，并使用以下SLURM脚本完成重复序列标记。现需编写类似脚本，实现多样本批量运行GATK4 HaplotypeCaller进行SNP检测。

参考的重复序列标记脚本 #

#!/bin/bash
#SBATCH --job-name=markdup
#SBATCH --array=1-11
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=4
#SBATCH --mem=30G
#SBATCH --time=02:00:00
#SBATCH --output=logs/markdup_%A_%a.out
#SBATCH --error=logs/markdup_%A_%a.err


module load apps/samtools/1.9/gcc-14.1.0
module load apps/picard/3.0.0/bin

# 获取样本信息
LINE=$(sed -n "${SLURM_ARRAY_TASK_ID}p" samples.tsv)
SAMPLE=$(echo $LINE | cut -d' ' -f1)

IN_BAM=${SAMPLE}.rg.bam
OUT_BAM=${SAMPLE}.dedup.bam
METRICS=${SAMPLE}.dedup.metrics.txt

if [[ ! -f "${IN_BAM}" ]]; then
  echo "缺失输入BAM文件: ${IN_BAM}" >&2
  exit 2
fi

# 运行Picard MarkDuplicates
picard MarkDuplicates \
    I=${IN_BAM} \
    O=${OUT_BAM} \
    M=${METRICS} \
    CREATE_INDEX=true \
    VALIDATION_STRINGENCY=SILENT \
    REMOVE_DUPLICATES=false

# 结果快速检查
echo "${SAMPLE}的重复序列统计信息:"
head -n 20 ${METRICS} || echo "未生成统计文件"

samtools flagstat ${OUT_BAM} || echo "samtools flagstat处理${OUT_BAM}失败"
samtools view -H ${OUT_BAM} | grep '^@RG' || echo "${OUT_BAM}中未找到@RG头信息" >&2

批量运行GATK4 HaplotypeCaller的SLURM脚本 #

#!/bin/bash
#SBATCH --job-name=gatk_hc
#SBATCH --array=1-11  # 对应samples.tsv中的样本数量，按需修改
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=4  # 与--native-pair-hmm-threads参数对应，按需调整
#SBATCH --mem=32G
#SBATCH --time=04:00:00  # 根据样本基因组大小调整，大基因组可延长
#SBATCH --output=logs/gatk_hc_%A_%a.out  # 输出日志路径，确保logs目录存在
#SBATCH --error=logs/gatk_hc_%A_%a.err
#SBATCH --mail-type=END,FAIL  # 任务结束/失败时发送邮件，按需开启


# 加载依赖模块，根据集群环境修改
module load apps/java/18.0.1.1/noarch
module load apps/gatk/4.2.2.0/noarch
module load apps/samtools/1.9/gcc-14.1.0


# 从samples.tsv中获取当前样本名（文件每行一个样本名）
LINE=$(sed -n "${SLURM_ARRAY_TASK_ID}p" samples.tsv)
SAMPLE=$(echo $LINE | cut -d' ' -f1)

# 定义输入输出文件路径
BAM=${SAMPLE}.dedup.bam  # 输入为去重后的BAM文件，需确保已生成对应的索引文件(.bai)
GVCF=${SAMPLE}.g.vcf.gz  # 输出的GVCF文件，适合后续联合基因分型
REFERENCE=GCF_019359855.2_Kyuss_2.0_genomic.fna  # 参考基因组路径，需确保已生成索引(.fai)和dict文件

# 检查输入BAM文件是否存在
if [[ ! -f "${BAM}" ]]; then
  echo "缺失输入BAM文件: ${BAM}" >&2
  exit 2
fi

# 运行GATK4 HaplotypeCaller
gatk --java-options "-Xmx28G" HaplotypeCaller \
  -R ${REFERENCE} \
  -I ${BAM} \
  -O ${GVCF} \
  -ERC GVCF \  # 生成GVCF格式，支持后续多样本联合分型
  --native-pair-hmm-threads ${SLURM_CPUS_PER_TASK}  # 利用分配的CPU核心加速

关键注意事项 #

• samples.tsv格式：文件每行存放一个样本名称，确保无多余空格，数量与--array=1-11中的数字一致。
• 模块适配：根据自己集群的软件环境修改module load命令，确保Java、GATK、Samtools版本兼容。
• 资源参数调整：--mem设置需预留足够内存给Java（-Xmx建议比总内存少4G左右），--time根据黑麦草基因组大小和样本数据量调整，若数据量大可延长至8-12小时。
• 参考基因组准备：需提前为参考基因组生成索引文件（samtools faidx）和序列字典（gatk CreateSequenceDictionary），否则HaplotypeCaller会报错。
• 输入BAM要求：去重后的BAM文件必须带有RG标签且已生成索引（.bai），否则会影响变异检测结果。

内容的提问来源于stack exchange，提问作者bioinfonext