引言
今天给大家介绍的是基于 Sentieon 软件开发的用于猪全基因组测序数据的自动化流程脚本。该流程实现了从原始测序数据（FASTQ）到变异检测结果（GVCF）以及 joint calling 的完整分析流程，支持多个测序平台和输出格式。
脚本支持原始测序数据（raw_fastq）、过滤后的测序数据（clean_fastq），进行质控、比对、排序、标记重复、生成质量评估指标，最终通过 Haplotyper 算法进行变异检测输出 gVCF 文件，通过 joint calling 输出 VCF 文件。
测试猪样本测序深度 9.98X，从 FASTQ 到 VCF 全流程分析最快用时 14.74 分钟，大幅缩短了猪全基因组 WGS 分析时间，有效加快动物的分子育种进程。
感谢 Ampere Computing LLC 和比亚迪服务器部门对本次测试的大力支持！！！
1. 环境设置与参数解析
1.1 Sentieon 环境配置
1.1.1 下载地址
软件地址链接
模型下载链接（请转到 VX 查看链接）
1.1.2 环境设置变量（点击跳转）
1.2 变量定义与参数解析
#!/bin/bashecho $0 \$SAMPLEID \$WORKDIR \$FASTQ_1 \$FASTQ_2 \$FASTA \$SUFFIX \$DATATYPE \$KEEP_CLEAN \$KEEP_BAM \$PLOIDYset -euxo pipefail
复制代码
echo：打印脚本名及输入参数占位符，方便调试时确认参数顺序。
set -euxo pipefail：设置脚本执行规则，增强健壮性。
export SAMPLEID=$1export WORKDIR=$2export FASTQ_1=$3export FASTQ_2=$4export FASTA=$5BSUFFIX=$6TYPE=${7:-"raw"}KEEP_CLEAN=${8:-keep}KEEP_BAM=${9:-keep}PLOIDY=${10:-2}LOGFILE=$SAMPLEID.run.logexport root=/Path/sentieon/202503/sentieon-genomics-202503/bin/
复制代码
从命令行读取参数，包括样本 ID、工作目录、输入文件路径等。
部分参数设置默认值。
LOGFILE：定义日志文件名（样本 ID+.run.log）。
root：sentieon 工具（基因组分析软件）的安装路径。
更详细的参数表格，如下表所示：
1.3 输入文件有效性检查
if [ "$BSUFFIX" = "bam" ] || [ "$BSUFFIX" = "cram" ] || [ "$BSUFFIX" = " " ]; then echo "$BSUFFIX check"else die "Error: check 6th blank, BSUFFIX must be 'bam' or 'cram' or space"fi
复制代码
检查 BSUFFIX，只能是 bam（二进制比对格式）、cram（压缩比对格式）或空格，否则报错退出。
1.4 测序平台判断
if [ -e $LOGFILE ];then export PLATFORM=$(awk '/Platform/{print $NF;exit}'$LOGFILE)else export count=$(zcat $FASTQ_1|head -n 1|awk '{print NF}') if [ $count -eq 1 ];then export PLATFORM="DNBSEQ" elif [ $count == 2 ];then export PLATFORM="ILLUMINA" elif [ $count == 3 ];then export PLATFORM="ILLUMINA" else echo"Unrecognized platform" export PLATFORM="ILLUMINA" fifi
复制代码
若日志文件已存在，从日志中提取测序平台（Platform 字段）。
若日志不存在，通过 FASTQ 文件首行的字段数判断平台。
1.5 输出文件格式确定
FAI=$FASTA.faiif [ "$BSUFFIX" = "bam" ] || [ "$BSUFFIX" = "cram" ];thenexport SUFFIX=$(awk -v preset=$BSUFFIX 'BEGIN{max=0}{if($2>max)max=$2}END{if(max>536870911){print "cram"}else{print preset}}' $FAI)elseexport SUFFIX=$(awk 'BEGIN{max=0}{if($2>max)max=$2}END{if(max>536870911){print "cram"}else{print "bam"}}' $FAI)fi
复制代码
根据参考基因组索引文件（$FASTA.fai）中最长序列的长度决定输出比对文件格式：
补充说明：对于单个染色体长度>536870911（约 512MB）的物种，Sentieon 软件可以切换至 cram，不用因 BAM 文件索引 (.bai) 的格式限制切割染色体。
1.6 环境与目录配置
export TMPDIR=$WORKDIRexport THREADS=$(nproc)[ -e $WORKDIR ]||mkdir -p $WORKDIRcd $WORKDIRexec >>$LOGFILE 2>&1echo "SampleID:" $SAMPLEIDecho "DataType:" $TYPE……
复制代码
配置临时目录、线程数，确保工作目录存在，并将所有输出写入日志文件。
1.7 设置计时函数(可选)
timer(){ start_time=$(date +%s) eval $2 && touch $3 end_time=$(date +%s) cost_time=$[ $end_time-$start_time ] echo -e "TIMER: $1\t$(($cost_time/60)) min $(($cost_time%60)) s"}
复制代码
用于记录每个分析步骤的开始/结束时间、耗时，并通过创建标记文件（如 qc.ok）标记步骤完成，避免重复执行。
2. 数据质控（Raw2clean）
raw2clean(){ cmd="fastp -w 16 -i $FASTQ_1 -I $FASTQ_2 -o $clean1 -O $clean2 -j $SAMPLEID.qc.json -h $SAMPLEID.qc.html&&rm $FASTQ_1 $FASTQ_2" timer raw2clean "$cmd" qc.ok}
复制代码
raw2clean(){ cmd="fastp -w 16 -i $FASTQ_1 -o $clean1 -j $SAMPLEID.qc.json -h $SAMPLEID.qc.html&&rm $FASTQ_1" timer raw2clean "$cmd" qc.ok}
复制代码
使用 fastp 工具对原始 FASTQ 数据进行质控（过滤低质量 reads、接头等），输出过滤后的 FASTQ 文件及质控报告文件（JSON/HTML）。
成功后创建 qc.ok 标记文件。
3. 序列比对与排序（Alignment）
alignment(){ tag="@RG\tID:rg_$SAMPLEID\tSM:$SAMPLEID\tPL:$PLATFORM" cmd="sentieon bwa mem -R \"$tag\" -t $THREADS -K 10000000 -x $ML_MODEL/bwa.model $FASTA $clean1 $clean2|sentieon util sort --temp_dir $TMPDIR -r $FASTA -o $SAMPLEID.sorted.$SUFFIX -t $THREADS --sam2bam -i -" echo $cmd timer alignment "$cmd" align.ok}
复制代码
alignment(){ tag="@RG\tID:rg_$SAMPLEID\tSM:$SAMPLEID\tPL:$PLATFORM" cmd="$root/sentieon bwa mem -R \"$tag\" -t $THREADS -K 10000000 $FASTA -x $ML_MODEL/bwa.model $clean1 |$root/sentieon util sort --temp_dir $TMPDIR -r $FASTA -o $SAMPLEID.sorted.$SUFFIX -t $THREADS --sam2bam -i -" echo $cmd timer alignment "$cmd" align.ok}
复制代码
使用 sentieon bwa turbo 进行序列比对（基于 BWA 算法），添加 Read Group 信息（用于后续变异分析），并通过 sentieon util sort 对结果排序，输出 sorted.bam 或 sorted.cram。
基于测序平台选择不同的机器学习模型（$ML_MODEL）优化比对。
成功后创建 align.ok 标记文件。
4. 生成质量评估指标（Metrics）
metrics(){ cmd="sentieon driver --temp_dir $TMPDIR -r $FASTA -t $THREADS -i $SAMPLEID.sorted.$SUFFIX --algo WgsMetricsAlgo $SAMPLEID.WGS_METRICS.txt --algo MeanQualityByCycle $SAMPLEID.mq_metrics.txt --algo QualDistribution $SAMPLEID.qd_metrics.txt --algo GCBias --summary $SAMPLEID.gc_summary.txt $SAMPLEID.gc_metrics.txt --algo AlignmentStat $SAMPLEID.aln_metrics.txt --algo BaseDistributionByCycle $SAMPLEID.bd_metrics.txt --algo QualityYield $SAMPLEID.qy_metrics.txt --algo InsertSizeMetricAlgo $SAMPLEID.is_metrics.txt" timer metrics "$cmd" metrics.ok}
复制代码
使用 sentieon driver 计算多种测序质量指标，包括：
结果输出到多个质量指标的.txt 文件，成功后创建 metrics.ok。
5. 标记重复序列（Dedup）
dedup(){ cmd="sentieon driver -r $FASTA --temp_dir $TMPDIR -t $THREADS -i $SAMPLEID.sorted.$SUFFIX --algo LocusCollector --fun score_info $SAMPLEID.score.txt&&sentieon driver -r $FASTA -t $THREADS -i $SAMPLEID.sorted.$SUFFIX --algo Dedup --score_info $SAMPLEID.score.txt --metrics $SAMPLEID.dedup_metrics.txt $SAMPLEID.deduped.$SUFFIX&&rm $SAMPLEID.sorted.$SUFFIX* $SAMPLEID.score.txt*" timer markdup "$cmd" markdup.ok}
复制代码
标记重复序列分两步：
完成后删除中间文件（排序后的比对文件），创建 markdup.ok。
6. 碱基质量重校正(BQSR 可选)
bqsr(){ sentieon driver -t $THREADS -r $FASTA -i $SAMPLEID.deduped.$SUFFIX --algo QualCal -k $KNOWN_SITES $SAMPLEID.RECAL_DATA.TABLE sentieon driver -t $THREADS -r $FASTA -i $SAMPLEID.deduped.$SUFFIX -q $SAMPLEID.RECAL_DATA.TABLE --algo QualCal -k $KNOWN_SITES $SAMPLEID.RECAL_DATA.TABLE.POST --algo ReadWriter $SAMPLEID.deduped.RECALIBRATED.$SUFFIX timer markdup "$cmd" markdup.ok}
复制代码
此过程分为两个主要部分：首先创建一个校正模型（生成 .TABLE 文件），然后应用该模型生成报告和碱基质量重校正后的 bam。
--algo QualCal: 指定运行的算法是 QualCal（质量校正算法）。
KNOWN_SITES：已知变异数据库的 VCF 文件路径。
RECAL_DATA.TABLE：重校准表的位置和文件名。
RECAL_DATA.TABLE.POST: 临时性的后重校准表的位置和文件名。
--algo ReadWriter：这一步为可选的。Sentieon® 变异检测可以在运行时使用校正前的 BAM 加上重校准表来即时执行重校正。便可以不输出巨大的 BAM 文件，节省磁盘的空间。
若选择执行这一步，第七部分变异检测中请将 $SAMPLEID.deduped.RECALIBRATED.$SUFFIX 设为输入文件。
7. 变异检测（DNAseq）
dnaseq(){ cmd="sentieon driver --temp_dir $TMPDIR -r $FASTA -t $THREADS -i $SAMPLEID.deduped.$SUFFIX --algo Haplotyper --emit_conf=30 --call_conf=30 --emit_mode gvcf --ploidy $PLOIDY $SAMPLEID.hc.gvcf.gz && md5sum $SAMPLEID.hc.gvcf.gz > $SAMPLEID.hc.gvcf.gz.md5" timer Haplotyper "$cmd" hc.ok}
复制代码
使用 Haplotyper 算法进行单倍型分析，生成 GVCF 文件。
计算 GVCF 文件的 MD5 ，确保文件完整性，成功后创建 hc.ok。
8. Joint Calling
8.1 参数检查与使用说明
#!/bin/bash[ $# -eq 0 ]&&echo Usage: $(basename $0) \$FASTA \$GVCF_LIST \$NUM \$DATADIR&&exitstart_time=$(date +%s)
复制代码
如果脚本没有传入任何参数，则打印使用说明并退出。
使用方式：脚本名 FASTA 文件 GVCF 列表文件 NUM 数据目录
记录脚本开始执行的时间
8.2 设置错误处理
-e: 任何命令失败则退出脚本。
-u: 使用未定义的变量时报错。
-o pipefail: 管道中任何一个命令失败则整个管道失败。
8.3 参数赋值
DATADIR=$4FASTA=$1 #"$FASTA_DIR/genomeEN_split.fa"GVCF_LIST=$2NUM=$3
复制代码
$1: 参考基因组 FASTA 文件路径。
$2: 包含所有样本 GVCF 文件路径的列表文件。
$3: 用于命名输出文件的数字标识（如批次号）。
$4: 数据目录，用于存放输出文件。
8.4 设置线程数、工作目录和日志文件
NT=$(nproc) #number of threads to use in computation, set to number of cores in the serverWORKDIR="$DATADIR/JointCall-${NUM}"[ -e $WORKDIR ]||mkdir -p $WORKDIR#[ -e $file ]&&exit 0cd $WORKDIRLOGFILE=$WORKDIR/joint-call${NUM}_run.logexec >$LOGFILE 2>&1
复制代码
NT: 获取当前系统的 CPU 核心数。
WORKDIR: 根据 NUM 创建唯一的工作目录。
如果目录不存在则创建，并进入该目录。
将所有标准输出和标准错误重定向到日志文件。
8.5 执行联合变异检测
root=/APP/u22/x86_com/sentieon/202503/sentieon-genomics-202503cat $GVCF_LIST|$root/bin/sentieon driver -r $FASTA --algo GVCFtyper \ $WORKDIR/output${NUM}-jc.vcf.gz - || { echo echo "GVCFtyper failed"; exit 1; }
复制代码
cat $GVCF_LIST: 读取 GVCF 文件列表。
sentieon driver: 调用 Sentieon 驱动程序。
-r $FASTA: 指定参考基因组。
--algo GVCFtyper: 使用 GVCFtyper 算法进行联合变异检测。
output${NUM}-jc.vcf.gz: 输出的压缩 VCF 文件。
如果命令失败，输出错误信息并退出。
8.6 计算并输出运行时间（可选）
end_time=$(date +%s)cost_time=$[ $end_time-$start_time ]echo "joint calling time is $(($cost_time/60))min $(($cost_time%60))s"echo "the Joint-calling done at `date +%H:%M:%S` !!!"
复制代码
计算脚本运行的总时间，并以“分:秒”格式输出。最后输出完成时间。
9. 总结
该脚本实现了 DNA 测序数据从原始测序数据（FASTQ）到变异检测结果（GVCF）以及 joint calling 的完整分析流程，支持 Illumina/DNBSEQ 平台，可配置输出格式（BAM/CRAM）和中间文件保留策略，通过标记文件和日志确保流程可进行追溯，适用于大规模基因组变异检测场景。特点包括：
参数化：通过命令行参数灵活控制输入、输出和流程选项。
自动化：自动检测测序平台和参考基因组以决定最佳分析策略。
稳健性：使用 set -euxo pipefail 和标记文件实现错误处理和断点续跑。
模块化：将每个分析步骤封装成函数。
高效性：使用商业优化的 sentieon 工具替代金标准 GATK/BWA，速度更快。
可追溯性：详细的日志记录和 MD5 校验确保结果可重现。
要运行此脚本，需要预先安装好 sentieon、fastp 等软件，并准备好对应的模型文件 bundle。
10. 脚本应用示例
使用上述脚本对猪全基因组测序数据分析的测序结果，具体样本信息如下表所示：
参考基因组下载
wget -c -t 0 --wait=5 --progress=bar:force http://ftp.ensembl.org/pub/release-108/fasta/sus_scrofa/dna/Sus_scrofa.Sscrofa11.1.dna.toplevel.fa.gz
复制代码
测试硬件配置
测试结果
使用本文流程对猪的全基因组测序数据进行变异检测分析，下表为不同 CPU 核数下的计算时间和资源调用情况:
本次测试在不同的线程数上进行性能的比较。从数据中可以明显看出，随着线程数的增加，变异检测的整体效率显著提升。从 FastQ 到 gVCF 全流程分析最快用时 14.74 分钟，大幅缩短了猪的全基因组 WGS 分析时间，有效加快动物的分子育种进程。
Sentieon 在不断地优化算法的运行效率，为科研工作者提供更快速、更经济的基因检测方案。若您刚好有需要检测的数据，不妨来申请试用 Sentieon 吧！
更多下载链接等内容，请关注 VX：毅硕科技，了解更多内容。
评论