DeepVariant全基因组测序案例研究：从数据准备到结果评估-优快云博客

本文链接：https://blog.youkuaiyun.com/gitblog_00564/article/details/148577863

DeepVariant全基因组测序案例研究：从数据准备到结果评估

deepvariant DeepVariant is an analysis pipeline that uses a deep neural network to call genetic variants from next-generation DNA sequencing data. 项目地址: https://gitcode.com/gh_mirrors/de/deepvariant

前言

DeepVariant是Google开发的一款基于深度学习的变异检测工具，它能够从二代测序数据中准确识别单核苷酸变异(SNPs)和小片段插入缺失(INDELs)。本文将详细介绍如何在实际的全基因组测序(WGS)数据上应用DeepVariant，并通过基准测试评估其变异检测的准确性。

实验设计

本案例研究使用以下数据：

参考基因组：GRCh38无alt分析集
测序数据：HG003样本的35x覆盖度NovaSeq PCR-free数据
基准数据集：Genome in a Bottle (GIAB) v4.2.1基准集
分析范围：仅限20号染色体(chr20)

环境准备

软件依赖

运行DeepVariant需要以下工具：

Docker：用于容器化运行DeepVariant和hap.py
curl：用于下载参考基因组和测试数据

参考基因组下载

GRCh38参考基因组可从NCBI获取，我们使用无alt分析集版本：

mkdir -p reference
FTPDIR=ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/001/405/GCA_000001405.15_GRCh38/seqs_for_alignment_pipelines.ucsc_ids

curl ${FTPDIR}/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.gz | gunzip > reference/GRCh38_no_alt_analysis_set.fasta
curl ${FTPDIR}/GCA_000001405.15_GRCh38_no_alt_analysis_set.fna.fai > reference/GRCh38_no_alt_analysis_set.fasta.fai

基准数据集准备

GIAB基准集用于评估DeepVariant的变异检测准确性：

mkdir -p benchmark
FTPDIR=ftp://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/release/AshkenazimTrio/HG003_NA24149_father/NISTv4.2.1/GRCh38

curl ${FTPDIR}/HG003_GRCh38_1_22_v4.2.1_benchmark_noinconsistent.bed > benchmark/HG003_GRCh38_1_22_v4.2.1_benchmark_noinconsistent.bed
curl ${FTPDIR}/HG003_GRCh38_1_22_v4.2.1_benchmark.vcf.gz > benchmark/HG003_GRCh38_1_22_v4.2.1_benchmark.vcf.gz
curl ${FTPDIR}/HG003_GRCh38_1_22_v4.2.1_benchmark.vcf.gz.tbi > benchmark/HG003_GRCh38_1_22_v4.2.1_benchmark.vcf.gz.tbi

测序数据获取

HG003样本的20号染色体BAM文件：

mkdir -p input
HTTPDIR=https://storage.googleapis.com/deepvariant/case-study-testdata

curl ${HTTPDIR}/HG003.novaseq.pcr-free.35x.dedup.grch38_no_alt.chr20.bam > input/HG003.novaseq.pcr-free.35x.dedup.grch38_no_alt.chr20.bam
curl ${HTTPDIR}/HG003.novaseq.pcr-free.35x.dedup.grch38_no_alt.chr20.bam.bai > input/HG003.novaseq.pcr-free.35x.dedup.grch38_no_alt.chr20.bam.bai

DeepVariant运行流程

DeepVariant分析流程包含三个主要步骤：

make_examples：将BAM文件转换为TensorFlow示例
call_variants：使用深度学习模型预测变异
postprocess_variants：后处理生成VCF文件

单命令运行（CPU版本）

mkdir -p output
mkdir -p output/intermediate_results_dir

BIN_VERSION="1.6.1"

sudo docker run \
  -v "${PWD}/input":"/input" \
  -v "${PWD}/output":"/output" \
  -v "${PWD}/reference":"/reference" \
  google/deepvariant:"${BIN_VERSION}" \
  /opt/deepvariant/bin/run_deepvariant \
  --model_type WGS \
  --ref /reference/GRCh38_no_alt_analysis_set.fasta \
  --reads /input/HG003.novaseq.pcr-free.35x.dedup.grch38_no_alt.chr20.bam \
  --output_vcf /output/HG003.output.vcf.gz \
  --output_gvcf /output/HG003.output.g.vcf.gz \
  --num_shards $(nproc) \
  --regions chr20 \
  --intermediate_results_dir /output/intermediate_results_dir

关键参数说明：

--model_type WGS：指定使用全基因组测序模型
--num_shards $(nproc)：使用所有可用的CPU核心并行处理
--regions chr20：限制分析范围到20号染色体
--intermediate_results_dir：保存中间结果（可选）

结果评估

使用hap.py工具将DeepVariant结果与GIAB基准集进行比较：

mkdir -p happy

sudo docker pull jmcdani20/hap.py:v0.3.12

sudo docker run \
  -v "${PWD}/benchmark":"/benchmark" \
  -v "${PWD}/input":"/input" \
  -v "${PWD}/output":"/output" \
  -v "${PWD}/reference":"/reference" \
  -v "${PWD}/happy:/happy" \
  jmcdani20/hap.py:v0.3.12 /opt/hap.py/bin/hap.py \
  /benchmark/HG003_GRCh38_1_22_v4.2.1_benchmark.vcf.gz \
  /output/HG003.output.vcf.gz \
  -f /benchmark/HG003_GRCh38_1_22_v4.2.1_benchmark_noinconsistent.bed \
  -r /reference/GRCh38_no_alt_analysis_set.fasta \
  -o /happy/happy.output \
  --engine=vcfeval \
  --pass-only \
  -l chr20