我有一个 hg19 对齐的 BAM,我希望为其生成一个 DeepVariant VCF。我使用 samtools 提取标题并确保 hg19 参考 FASTA 索引包含相同的重叠群和位置。我最初的目标是在这个 WGS BAM 上只运行一个外显子组模型,使用以下模型和区域:
MODEL=gs://deepvariant/models/DeepVariant/0.7.2/DeepVariant-inception_v3-0.7.2+data-wes_standard
--regions gs://deepvariant/exome-case-study-testdata/refseq.coding_exons.b37.extended50.bed
不幸的是,脚本抗议说 BED 和 BAM / FASTA 参考之间的匹配为 0。我决定运行相同的外显子组模型,但没有指定区域。这是我的脚本:
#!/bin/bash
set -euo pipefail
# Set common settings.
PROJECT_ID=<MY PROJECT>
OUTPUT_BUCKET=gs://<MY BUCKET>
STAGING_FOLDER_NAME=staging
OUTPUT_FILE_NAME=output.vcf
# Model for calling whole genome sequencing data.
MODEL=gs://deepvariant/models/DeepVariant/0.7.2/DeepVariant-inception_v3-0.7.2+data-wes_standard
IMAGE_VERSION=0.7.2
DOCKER_IMAGE=gcr.io/deepvariant-docker/deepvariant:"${IMAGE_VERSION}"
COMMAND="/opt/deepvariant_runner/bin/gcp_deepvariant_runner \
--project ${PROJECT_ID} \
--zones us-west1-* \
--docker_image ${DOCKER_IMAGE} \
--outfile ${OUTPUT_BUCKET}/${OUTPUT_FILE_NAME} \
--staging ${OUTPUT_BUCKET}/${STAGING_FOLDER_NAME} \
--model ${MODEL} \
--bam gs://my-bucket/wgs_data.bam \
--ref gs://my-bucket/human_g1k_v37.fa \
--shards 512 \
--make_examples_workers 32 \
--make_examples_cores_per_worker 16 \
--make_examples_ram_per_worker_gb 60 \
--make_examples_disk_per_worker_gb 200 \
--call_variants_workers 32 \
--call_variants_cores_per_worker 32 \
--call_variants_ram_per_worker_gb 60 \
--call_variants_disk_per_worker_gb 50 \
--gcsfuse"
# Run the pipeline.
gcloud alpha genomics pipelines run \
--project "${PROJECT_ID}" \
--service-account-scopes="https://www.googleapis.com/auth/cloud-platform" \
--logging "${OUTPUT_BUCKET}/${STAGING_FOLDER_NAME}/runner_logs_$(date +%Y%m%d_%H%M%S).log" \
--regions us-west1 \
--docker-image gcr.io/cloud-genomics-pipelines/gcp-deepvariant-runner \
--command-line "${COMMAND}"
BAM 有对应的 BAI,FA 有 FAI 文件。DeepVariant QuickStart 表明这些设置将在 1-2 小时内生成 VCF,但我的管道现在已经运行了 7 多个小时。暂存文件夹现在有一个 call_variants,其中包含 32 个 GZ 文件中的 31 个。Genomics 管道视图显示了 11 个运行 call_variant 的管道,因此我怀疑它正在处理最后一个文件,准备将所有文件合并到一个 VCF 中。
我只是不明白为什么要花这么长时间。我故意排除了抢占式实例,并且文档说外显子组管道应该只需要 20 分钟(WGS 需要 1-2 小时)。为什么会这么慢?