文档建议反馈控制台
首页
学习
活动
专区
圈层
工具
MCP广场
文章/答案/技术大牛
发布
社区首页 >专栏 >Sentieon | 猪全基组（WGS）分析流程

Sentieon | 猪全基组（WGS）分析流程

原创
作者头像
INSVAST
发布于 2025-10-11 17:25:33
发布于 2025-10-11 17:25:33
1540
举报

引言

今天给大家介绍的是基于 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  变量定义与参数解析

代码语言：javascript
复制
#!/bin/bash
echo $0 \$SAMPLEID \$WORKDIR \$FASTQ_1 \$FASTQ_2 \$FASTA \$SUFFIX \$DATATYPE \$KEEP_CLEAN \$KEEP_BAM \$PLOIDY
set -euxo pipefail

echo：打印脚本名及输入参数占位符，方便调试时确认参数顺序。

set -euxo pipefail：设置脚本执行规则，增强健壮性。

代码语言：javascript
复制
export SAMPLEID=$1
export WORKDIR=$2
export FASTQ_1=$3
export FASTQ_2=$4
export FASTA=$5
BSUFFIX=$6
TYPE=${7:-"raw"}
KEEP_CLEAN=${8:-keep}
KEEP_BAM=${9:-keep}
PLOIDY=${10:-2}
LOGFILE=$SAMPLEID.run.log
export root=/Path/sentieon/202503/sentieon-genomics-202503/bin/

从命令行读取参数，包括样本 ID、工作目录、输入文件路径等。

部分参数设置默认值。

LOGFILE：定义日志文件名（样本 ID+.run.log）。

root：sentieon工具（基因组分析软件）的安装路径。

更详细的参数表格，如下表所示：

1.3  输入文件有效性检查

代码语言：javascript
复制
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  测序平台判断

代码语言：javascript
复制
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"
    fi
fi

若日志文件已存在，从日志中提取测序平台（Platform字段）。

若日志不存在，通过 FASTQ 文件首行的字段数判断平台。

1.5  输出文件格式确定

代码语言：javascript
复制
FAI=$FASTA.fai
if [ "$BSUFFIX" = "bam" ] || [ "$BSUFFIX" = "cram" ];then
export SUFFIX=$(awk -v preset=$BSUFFIX 'BEGIN{max=0}{if($2>max)max=$2}END{if(max>536870911){print "cram"}else{print preset}}' $FAI)
else
export 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），强制使用cram（更适合大基因组压缩）。
  • 否则根据BSUFFIX或默认bam。

补充说明：对于单个染色体长度>536870911（约512MB）的物种，Sentieon软件可以切换至cram，不用因BAM 文件索引 (.bai) 的格式限制切割染色体。

1.6  环境与目录配置

代码语言：javascript
复制
export TMPDIR=$WORKDIR
export THREADS=$(nproc)
[ -e $WORKDIR ]||mkdir -p $WORKDIR
cd $WORKDIR
exec >>$LOGFILE 2>&1
echo "SampleID:" $SAMPLEID
echo "DataType:" $TYPE
……

配置临时目录、线程数，确保工作目录存在，并将所有输出写入日志文件。

1.7  设置计时函数(可选)

代码语言：javascript
复制
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）

  • 双端测序质控脚本：
代码语言：javascript
复制
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
}
  • 单端测序质控脚本：
代码语言：javascript
复制
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）

  • 双端测序比对脚本：
代码语言：javascript
复制
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
}
  • 单端测序比对脚本：
代码语言：javascript
复制
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）

代码语言：javascript
复制
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计算多种测序质量指标，包括：

  • 全基因组测序指标（WgsMetricsAlgo）。
  • 碱基质量分布、GC 偏差、插入片段长度等。

结果输出到多个质量指标的.txt文件，成功后创建metrics.ok。

5.  标记重复序列（Dedup）

代码语言：javascript
复制
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
}

标记重复序列分两步：

  • LocusCollector：收集位点质量分数信息。
  • Dedup：基于分数标记并去除 PCR 重复序列。

完成后删除中间文件（排序后的比对文件），创建markdup.ok。

6.  碱基质量重校正(BQSR可选)

代码语言：javascript
复制
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）

代码语言：javascript
复制
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  参数检查与使用说明

代码语言：javascript
复制
#!/bin/bash
[ $# -eq 0 ]&&echo Usage: $(basename $0) \$FASTA \$GVCF_LIST  \$NUM \$DATADIR&&exit
start_time=$(date +%s)

如果脚本没有传入任何参数，则打印使用说明并退出。

使用方式：脚本名 FASTA文件 GVCF列表文件 NUM 数据目录

记录脚本开始执行的时间

8.2  设置错误处理

代码语言：javascript
复制
set -euo pipefail

-e: 任何命令失败则退出脚本。

-u: 使用未定义的变量时报错。

-o pipefail: 管道中任何一个命令失败则整个管道失败。

8.3  参数赋值

代码语言：javascript
复制
DATADIR=$4
FASTA=$1 #"$FASTA_DIR/genomeEN_split.fa"
GVCF_LIST=$2
NUM=$3

$1: 参考基因组 FASTA 文件路径。

$2: 包含所有样本 GVCF 文件路径的列表文件。

$3: 用于命名输出文件的数字标识（如批次号）。

$4: 数据目录，用于存放输出文件。

8.4  设置线程数、工作目录和日志文件

代码语言：javascript
复制
NT=$(nproc) #number of threads to use in computation, set to number of cores in the server
WORKDIR="$DATADIR/JointCall-${NUM}"
[ -e $WORKDIR ]||mkdir -p $WORKDIR
#[ -e $file ]&&exit 0
cd $WORKDIR
LOGFILE=$WORKDIR/joint-call${NUM}_run.log
exec >$LOGFILE 2>&1

NT: 获取当前系统的 CPU 核心数。

WORKDIR: 根据 NUM 创建唯一的工作目录。

如果目录不存在则创建，并进入该目录。

将所有标准输出和标准错误重定向到日志文件。

8.5  执行联合变异检测

代码语言：javascript
复制
root=/APP/u22/x86_com/sentieon/202503/sentieon-genomics-202503
cat $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  计算并输出运行时间（可选）

代码语言：javascript
复制
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.  脚本应用示例

使用上述脚本对猪全基因组测序数据分析的测序结果，具体样本信息如下表所示：

参考基因组下载

代码语言：javascript
复制
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为单颗AmpereOne A192-32X
  • 内存为512GB DDR5
  • 系统为Ubuntu 24.04/Kernel 6.8

测试结果

使用本文流程对猪的全基因组测序数据进行变异检测分析，下表为不同CPU核数下的计算时间和资源调用情况:

本次测试在不同的线程数上进行性能的比较。从数据中可以明显看出，随着线程数的增加，变异检测的整体效率显著提升。从FastQ到gVCF全流程分析最快用时14.74分钟，大幅缩短了猪的全基因组WGS分析时间，有效加快动物的分子育种进程。

Sentieon在不断地优化算法的运行效率，为科研工作者提供更快速、更经济的基因检测方案。若您刚好有需要检测的数据，不妨来申请试用Sentieon吧！

更多下载链接等内容，请关注VX：毅硕科技，了解更多内容。

原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。

如有侵权，请联系 cloudcommunity@tencent.com 删除。

生物基因

原创声明：本文系作者授权腾讯云开发者社区发表，未经许可，不得转载。

如有侵权，请联系 cloudcommunity@tencent.com 删除。

生物基因
#Sentieon
#全基因组
#WGS
#标准分析流程
#猪
评论
登录后参与评论
0 条评论
热度
最新
登录 后参与评论
推荐阅读
目录
  • 引言
  • 1.  环境设置与参数解析
    • 1.1  Sentieon环境配置
    • 1.2  变量定义与参数解析
    • 1.3  输入文件有效性检查
    • 1.4  测序平台判断
    • 1.5  输出文件格式确定
    • 1.6  环境与目录配置
    • 1.7  设置计时函数(可选)
  • 2.  数据质控（Raw2clean）
  • 3.  序列比对与排序（Alignment）
  • 4.  生成质量评估指标（Metrics）
  • 5.  标记重复序列（Dedup）
  • 6.  碱基质量重校正(BQSR可选)
  • 7.  变异检测（DNAseq）
  • 8.  Joint Calling
    • 8.1  参数检查与使用说明
    • 8.2  设置错误处理
    • 8.3  参数赋值
    • 8.4  设置线程数、工作目录和日志文件
    • 8.5  执行联合变异检测
    • 8.6  计算并输出运行时间（可选）
  • 9.  总结
  • 10.  脚本应用示例
领券

腾讯云开发者

扫码关注腾讯云开发者

扫码关注腾讯云开发者

领取腾讯云代金券

热门产品

热门推荐

更多推荐

Copyright © 2013 - 2025 Tencent Cloud. All Rights Reserved. 腾讯云 版权所有 

深圳市腾讯计算机系统有限公司 ICP备案/许可证号：粤B2-20090059 深公网安备号 44030502008569

腾讯云计算（北京）有限责任公司 京ICP证150476号 |  京ICP备11018762号 | 京公网安备号11010802020287

问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档

Copyright © 2013 - 2025 Tencent Cloud.

All Rights Reserved. 腾讯云 版权所有

登录 后参与评论