giraffe (vg giraffe + GATK)
请注意,Parabricks GPU 加速的 Giraffe 工具目前处于 beta 测试阶段。
使用泛基因组比对器 VG Giraffe [1] [2] 生成一个或一对 FASTQ 文件的 BAM 输出。
有关所有可用选项的详细列表,请参阅 giraffe 参考部分。
VG Giraffe 是由加州大学圣克鲁兹分校 (UCSC) Benedict Paten 博士实验室开发的短读比对工具。这款创新工具可以将 reads 比对到多个参考基因组的图表示,从而提高下游分析的质量。通过将 reads 同时准确地比对到数千个基因组,VG Giraffe 相对于传统的单参考基因组比对器有了显著的改进。
通过利用基于图的方法,VG Giraffe 可以更有效地处理跨人群的遗传多样性和结构变异。以下是使用 VG Giraffe 的三个主要优势
提高准确性:与线性基因组比对器相比,VG Giraffe 在 reads 比对中实现了更高的精度和召回率,尤其是在处理复杂的基因组区域或具有显著遗传多样性的人群时。
减少参考偏倚(或比对偏倚):通过将多个单倍型和已知变异纳入其图结构,VG Giraffe 最大程度地减少了传统线性基因组比对器中固有的参考偏倚。这使得能够更全面和公正地表征遗传变异,特别是对于与标准参考基因组差异显著的样本。
更快的性能:尽管 VG Giraffe 使用更复杂的图结构,但它比其前身 VG Map 快得多,并且速度与流行的线性基因组比对器相当。它可以将测序 reads 比对到数千个人类基因组,速度与比对到单个参考基因组的方法相似。
VG Giraffe 可以在 Parabricks 中使用,Parabricks 是一套为基因组学加速二级分析而设计的软件套件。我们的封装器 (pbrun giraffe
) 将运行我们的 GPU 加速 VG Giraffe,并按坐标对输出 BAM 进行排序。
虽然用户可以使用 VG Autoindex 工具为 VG Giraffe 构建自定义参考图,但预构建的泛基因组图也可用。Paten 博士的实验室和人类泛基因组联盟已将这些资源公开,使研究人员能够利用高质量、即用型的泛基因组图进行分析 (HPRC 数据)。
下面测试中使用的索引文件集可以使用 aws cli
如下方式下载
aws s3 cp \
s3://human-pangenomics/pangenomes/freeze/freeze1/minigraph-cactus/filtered . \
--no-sign-request --recursive --exclude "*" --include "hprc-v1.0-mc-grch38-minaf.0.1.*"
然后,要从索引文件所在的目录生成所需的 .gbz
文件
docker run --rm \
-v $(pwd):/workdir \
--workdir /workdir \
quay.io/vgteam/vg:v1.59.0 \
vg gbwt --gbz-format \
-g hprc-v1.0-mc-grch38-minaf.0.1.gbz \
-I hprc-v1.0-mc-grch38-minaf.0.1.gg hprc-v1.0-mc-grch38-minaf.0.1.gbwt
在运行时,索引数据被加载到 GPU 设备内存中。由于索引的大小,这可能会限制其他操作的可用内存。以下选项会影响设备内存使用量和性能
对于 16GB 设备(例如 T4):使用
--low-memory
选项对于 16GB-40GB 设备(例如 L4、A10):通过调整以下参数优化性能
--nstreams
: 控制每个 GPU 的 CUDA 流数量--batch-size
: 调整批处理中处理的 reads 数量对于 L4,使用
--nstreams 2 --batch-size 10000
获得了最佳性能
对于 >40GB 设备:默认参数已足够;但是,通过调整上述参数,仍有进一步优化的潜力。
注意:虽然每个设备都存在固定的基本内存分配,但流的数量和批处理大小是影响总设备内存消耗的主要因素。
# This command assumes all the inputs are in INPUT_DIR and all the outputs go to OUTPUT_DIR.
docker run --rm --gpus all --volume INPUT_DIR:/workdir --volume OUTPUT_DIR:/outputdir \
--workdir /workdir \
nvcr.io/nvidia/clara/clara-parabricks:4.4.0-1 \
pbrun giraffe --read-group "sample_rg1" \
--sample "sample-name" --read-group-library "library" \
--read-group-platform "platform" --read-group-pu "pu" \
--dist-name /workdir/hprc-v1.0-mc-grch38-minaf.0.1.dist \
--minimizer-name /workdir/hprc-v1.0-mc-grch38-minaf.0.1.min \
--gbz-name /workdir/hprc-v1.0-mc-grch38-minaf.0.1.giraffe.gbz \
--xg-name /workdir/hprc-v1.0-mc-grch38-minaf.0.1.xg \
--graph-name /workdir/hprc-v1.0-mc-grch38-minaf.0.1.gg \
--gbwt-name /workdir/hprc-v1.0-mc-grch38-minaf.0.1.gbwt \
--in-fq /workdir/${INPUT_FASTQ_1} /workdir/${INPUT_FASTQ_2} \
--out-bam /outputdir/${OUTPUT_BAM}
要将 Giraffe 比对的 BAM 文件用于变异检测,您需要从 Giraffe 索引文件中提取相应的参考文件。从包含 Giraffe 索引文件的目录运行以下命令
# Extract the list of paths corresponding to GRCh38
docker run --rm \
-v $(pwd):/workdir \
--workdir /workdir \
quay.io/vgteam/vg:v1.59.0 \
vg paths -x hprc-v1.0-mc-grch38-minaf.0.1.xg -L > hprc-v1.0-mc-grch38-minaf.0.1.paths
# Filter paths list
grep -v _decoy hprc-v1.0-mc-grch38-minaf.0.1.paths \
| grep -v _random \
| grep -v chrUn_ \
| grep -v chrEBV \
| grep -v chrM \
| grep -v chain_ > hprc-v1.0-mc-grch38-minaf.0.1.paths.sub
# Extract the corresponding sequences to a FASTA file
docker run --rm \
-v $(pwd):/workdir \
--workdir /workdir \
quay.io/vgteam/vg:v1.59.0 \
vg paths -x hprc-v1.0-mc-grch38-minaf.0.1.xg -p hprc-v1.0-mc-grch38-minaf.0.1.paths.sub -F > hprc-v1.0-mc-grch38-minaf.0.1.fa
# Index the fasta file
samtools faidx hprc-v1.0-mc-grch38-minaf.0.1.fa
这些命令将生成一个 FASTA 文件 (hprc-v1.0-mc-grch38-minaf.0.1.fa
) 和相应的索引 (hprc-v1.0-mc-grch38-minaf.0.1.fa.fai
),它们可以用作变异检测的参考。请注意,这些文件也可用于 BQSR (bqsr)。
获得 Giraffe 比对的 BAM 文件和提取的参考 FASTA 后,您可以使用 DeepVariant 或 HaplotypeCaller 继续进行变异检测。
# Deepvariant
# This command assumes all the inputs are in INPUT_DIR and all the outputs go to OUTPUT_DIR.
docker run --rm --gpus all --volume INPUT_DIR:/workdir --volume OUTPUT_DIR:/outputdir \
--workdir /workdir \
nvcr.io/nvidia/clara/clara-parabricks:4.4.0-1 \
pbrun deepvariant \
--ref /workdir/hprc-v1.0-mc-grch38-minaf.0.1.fa \
--in-bam /workdir/${INPUT_BAM} \
--out-variants /outputdir/${OUTPUT_VCF}
# Haplotype Caller
# This command assumes all the inputs are in INPUT_DIR and all the outputs go to OUTPUT_DIR.
docker run --rm --gpus all --volume INPUT_DIR:/workdir --volume OUTPUT_DIR:/outputdir \
--workdir /workdir \
nvcr.io/nvidia/clara/clara-parabricks:4.4.0-1 \
pbrun haplotypecaller \
--ref /workdir/hprc-v1.0-mc-grch38-minaf.0.1.fa \
--in-bam /workdir/${INPUT_BAM} \
--in-recal-file /workdir/${INPUT_RECAL_FILE} \
--out-variants /outputdir/${OUTPUT_VCF}
有关变异检测的更详细说明,请参阅特定工具的文档 (deepvariant, haplotypecaller)。
以下命令是上述 Parabricks 命令的 vg-1.59.0 和 GATK4 对等命令。这些命令的输出将与上述命令的输出相同。有关结果比较,请参阅输出比较页面。
# Run giraffe and pipe the output to create a sorted BAM.
$ vg giraffe \
-t 16 \
-d /workdir/hprc-v1.0-mc-grch38-minaf.0.1.dist \
-m /workdir/hprc-v1.0-mc-grch38-minaf.0.1.min \
-x /workdir/hprc-v1.0-mc-grch38-minaf.0.1.xg \
-g /workdir/hprc-v1.0-mc-grch38-minaf.0.1.gg \
-H /workdir/hprc-v1.0-mc-grch38-minaf.0.1.gbwt \
-f /workdir/${INPUT_FASTQ_1} \
-f /workdir/${INPUT_FASTQ_2} \
--output-format bam | \
gatk SortSam \
--java-options -Xmx30g \
--MAX_RECORDS_IN_RAM 5000000 \
-I /dev/stdin \
-O cpu.bam \
--SORT_ORDER coordinate
将输出与 CPU 对等项进行比较时,以下可能是造成细微差异的来源。
基线 VG 容器
在比较基线 giraffe 和 Parabricks 加速版本的输出时,如果您打算使用基线 vg 容器 (quay.io/vgteam/vg:v1.59.0
),则需要使用 Ubuntu 22.04 基础重建容器。这是因为底层操作系统默认 gcc 版本的 C++ 标准库发生了变化 (#4391)。修改其 Dockerfile 的第 6 行,将 20.04 替换为引用 mirror.gcr.io/library/ubuntu:22.04
。使用以下命令重建容器。
git clone https://github.com/vgteam/vg.git
cd vg
git checkout v1.59.0
git submodule update --init --recursive
make version
docker build --no-cache -f Dockerfile --build-arg THREADS=16 --tag \
<YOUR_CONTIANER_NAME> --network host ./
未比对的 reads
Parabricks
giraffe
对未比对的 reads 的排序方式与基线 GATK SortSam 略有不同。可以使用 samtools 过滤未比对的 reads,方法是运行samtools view -F 4
。
将 reads 比对到泛基因组图。
输入/输出文件选项
- --in-fq [IN_FQ ...]
-
配对末端 FASTQ 文件的路径。文件必须为 fastq 或 fastq.gz 格式。示例 1:--in-fq sampleX_1_1.fastq.gz sampleX_1_2.fastq.gz。(默认值:None)
- --in-se-fq [IN_SE_FQ ...]
-
单端 FASTQ 文件的路径。文件必须为 fastq 或 fastq.gz 格式。(默认值:None)
- -d DIST_NAME, --dist-name DIST_NAME
-
使用此距离索引进行聚类。(默认值:None)
选项为必填项。
- -m MINIMIZER_NAME, --minimizer-name MINIMIZER_NAME
-
使用此 minimizer 索引。(默认值:None)
选项为必填项。
- -Z GBZ_NAME, --gbz-name GBZ_NAME
-
映射到此 GBZ 图。(默认值:None)
选项为必填项。
- -x XG_NAME, --xg-name XG_NAME
-
用于 BAM 输出的 XG 图。(默认值:None)
选项为必填项。
- -g GRAPH_NAME, --graph-name GRAPH_NAME
-
用于映射的 GBWTGraph。(默认值:None)
选项为必填项。
- -H GBWT_NAME, --gbwt-name GBWT_NAME
-
用于映射的 GBWT 索引。(默认值:None)
选项为必填项。
- --out-bam OUT_BAM
-
输出 BAM 文件的路径。(默认值:None)
选项为必填项。
工具选项
- --read-group READ_GROUP
-
此运行的读取组 ID。(默认值:None)
选项为必填项。
- --sample SAMPLE
-
此运行中读取组的样本 (SM) 标签。(默认值:None)
选项为必填项。
- --read-group-library READ_GROUP_LIBRARY
-
此运行中读取组的文库 (LB) 标签。(默认值:None)
- --read-group-platform READ_GROUP_PLATFORM
-
此运行中读取组的平台 (PL) 标签;指用于生成 reads 的平台/技术。(默认值:None)
- --read-group-pu READ_GROUP_PU
-
此运行中读取组的平台单元 (PU) 标签。(默认值:None)
- --prune-low-cplx
-
在线性格式重新比对期间修剪短的和低复杂度的锚点。(默认值:None)
- --max-fragment-length MAX_FRAGMENT_LENGTH
-
在估计片段长度分布时,假设片段长度应小于 INT。(默认值:None)
- --fragment-mean FRAGMENT_MEAN
-
强制片段长度分布具有此平均值。(默认值:None)
- --fragment-stdev FRAGMENT_STDEV
-
强制片段长度分布具有此标准偏差。(默认值:None)
- --align-only
-
在 vg-giraffe 比对后生成输出 BAM。输出将不会按坐标排序。(默认值:None)
- --monitor-usage
-
监控执行期间(在排序和后排序期间可用)的大致 CPU 利用率和主机内存使用量。(默认值:None)
- --max-read-length MAX_READ_LENGTH
-
用于 giraffe 和过滤 FASTQ 输入的最大 read 长度/大小(即序列长度)(默认值:480)
- --min-read-length MIN_READ_LENGTH
-
用于 giraffe 和过滤 FASTQ 输入的最小 read 长度/大小(即序列长度)(默认值:1)
性能选项
- --nstreams NSTREAMS
-
每个 GPU 要使用的流数量;注意:更多流会增加设备内存使用量。(默认值:3)
- --num-primary-cpus-per-gpu NUM_PRIMARY_CPUS_PER_GPU
-
每个 GPU 驱动其关联线程池的主 CPU 线程数。(默认值:16)
- --cpu-thread-pool CPU_THREAD_POOL
-
每个主 CPU 线程的处理线程数。(默认值:1)
- --batch-size BATCH_SIZE
-
用于处理比对的批处理大小。(默认值:10000)
- --write-threads WRITE_THREADS
-
用于写入和预排序输出的线程数。(默认值:4)
- --gpuwrite
-
使用一个 GPU 加速写入最终 BAM/CRAM。(默认值:None)
- --gpuwrite-deflate-algo GPUWRITE_DEFLATE_ALGO
-
选择与 --gpuwrite 一起使用的 nvCOMP DEFLATE 算法。请注意,这些选项与 CPU DEFLATE 选项不对应。有效选项为 1、2 和 4。选项 1 最快,而选项 2 和 4 的吞吐量逐渐降低,但压缩率更高。当用户未提供输入时(即 None),默认值为 1 (默认值:None)
- --gpusort
-
使用 GPU 加速排序和标记。(默认值:None)
- --use-gds
-
使用 GPUDirect Storage (GDS) 启用 GPU 内存和存储之间直接内存访问 (DMA) 传输的直接数据路径。必须与 --gpuwrite 同时使用。有关如何设置和使用 GPUDirect Storage 的信息,请参阅 Parabricks 文档 > 最佳性能。(默认值:None)
- --memory-limit MEMORY_LIMIT
-
排序和后排序期间的系统内存限制,以 GB 为单位。默认情况下,限制为系统总内存的一半。(默认值:62)
- --low-memory
-
使用低内存模式;将降低每个 GPU 的流数量并减小批处理大小。(默认值:None)
常用选项
- --logfile LOGFILE
-
日志文件的路径。如果未指定,消息将仅写入标准错误输出。(默认值:None)
- --tmp-dir TMP_DIR
-
临时文件将存储在其中的目录的完整路径。
- --with-petagene-dir WITH_PETAGENE_DIR
-
PetaGene 安装目录的完整路径。默认情况下,这应已安装在 /opt/petagene。使用此选项还需要通过设置 LD_PRELOAD 环境变量来预加载 PetaLink 库。可选地设置用于数据和凭据的 PETASUITE_REFPATH 和 PGCLOUD_CREDPATH 环境变量 (默认值:None)
- --keep-tmp
-
完成后不要删除存储临时文件的目录。
- --no-seccomp-override
-
不要覆盖 docker 的 seccomp 选项 (默认值:None)。
- --version
-
查看兼容的软件版本。
GPU 选项
- --num-gpus NUM_GPUS
-
运行要使用的 GPU 数量。将使用 GPU 0..(NUM_GPUS-1)。
Jouni Sirén 等人,《Pangenomics enables genotyping of known structural variants in 5202 diverse genomes》。Science 374, abg 8871 (2021)。DOI: 10.1126/science.abg8871
基线 VG Giraffe: https://github.com/vgteam/vg