giraffe (vg giraffe + GATK)

Generate BAM or CRAM output given one or a pair of FASTQ files using the pangenome aligner VG Giraffe ^{1 2}.

Refer to the giraffe Reference section for a detailed listing of all available options.

What is giraffe?

VG Giraffe is a short-read mapping tool developed by Dr. Benedict Paten’s lab at the University of California, Santa Cruz (UCSC). This innovative tool can align reads to a graph representation of multiple reference genomes, enhancing the quality of downstream analyses. By accurately mapping reads to thousands of genomes simultaneously, VG Giraffe offers a substantial improvement over traditional single-reference aligners.

Why giraffe?

By utilizing a graph-based approach, VG Giraffe can more effectively handle genetic diversity and structural variations across populations. Here are three key benefits of using VG Giraffe:

1. Improved accuracy: VG Giraffe achieves higher precision and recall in read mapping compared to linear genome aligners, especially when dealing with complex genomic regions or populations with significant genetic diversity.
2. Reduced reference bias (or mapping bias): By incorporating multiple haplotypes and known variants into its graph structure, VG Giraffe minimizes the reference bias inherent in traditional linear genome aligners. This leads to more comprehensive and unbiased characterization of genetic variation, especially for samples that diverge significantly from the standard reference genome.
3. Faster performance: Despite working with more complex graph structures, VG Giraffe is significantly faster than its predecessor VG Map and comparable in speed to popular linear genome mappers. It can map sequencing reads to thousands of human genomes at a speed similar to methods that map to a single reference genome.

How Should I Use giraffe in Parabricks?

VG Giraffe can be used within Parabricks, a software suite designed for accelerated secondary analysis in genomics. Our wrapper (pbrun giraffe) will run our GPU-accelerated VG Giraffe and sort the output BAM by coordinate.

While users can build custom reference graphs for VG Giraffe using the VG Autoindex tool, pre-built pangenome graphs are also available. Dr. Paten’s lab and the Human Pangenome Consortium have made these resources publicly accessible, allowing researchers to leverage high-quality, ready-to-use pangenome graphs for their analyses (HPRC data).

Generating Required Index Files

The index files .gbz, .dist, .min, and .zipcodes are required to run Giraffe. A reference paths file is also needed to define the set of paths used for BAM output.

Supported GBZ version

Parabricks Giraffe reads GBZ version 1. Starting with vg 1.72.0, vg autoindex and vg gbwt write GBZ version 2 by default, which stores the node sequences zstd-compressed. That format is not supported: loading such a graph fails right away with an error naming the conversion command below.

An existing version 2 graph can be converted with the same vg v1.75.1 image used for index generation:

docker run --rm --volume $(pwd):/workdir \
    --workdir /workdir \
    --user $(id -u):$(id -g) \
    quay.io/vgteam/vg:v1.75.1 \
    vg gbwt -Z graph.gbz -g graph-v1.gbz --gbz-v1

--gbz-v1 is part of the GBWTGraph construction step of vg gbwt: -Z reads the version 2 graph and -g writes the version 1 copy. The option was added in vg 1.74.0, so the conversion requires 1.74.0 or later even if the graph itself is older. Refer to GBZ versions in the upstream vg wiki for more details.

Conversion only changes how the node sequences are stored, so existing .dist, .min, and .zipcodes files built from the same graph stay valid.

Building the index files

The index files can be generated from a GBZ version 1 graph using vg autoindex. The following example uses the HPRC v1.1 Minigraph-Cactus pangenome graph aligned to GRCh38. That download is already version 1, and vg autoindex reads the graph without rewriting it, so the graph stays version 1 even when indexed with vg v1.75.1:

# Download GBZ
# https://s3-us-west-2.amazonaws.com/human-pangenomics/pangenomes/freeze/freeze1/minigraph-cactus/hprc-v1.1-mc-grch38/hprc-v1.1-mc-grch38.d9.gbz
aws s3 cp \
    s3://human-pangenomics/pangenomes/freeze/freeze1/minigraph-cactus/hprc-v1.1-mc-grch38/hprc-v1.1-mc-grch38.d9.gbz \
    . \
    --no-sign-request
# Extract index files from GBZ
docker run --rm --volume $(pwd):/workdir \
    --workdir /workdir \
    --user $(id -u):$(id -g) \
    quay.io/vgteam/vg:v1.75.1 \
    vg autoindex \
        -p hprc-v1.1-mc-grch38.d9.autoindex.1.75.1 \
        -G hprc-v1.1-mc-grch38.d9.gbz \
        -w giraffe
# Extract paths from GBZ
docker run --rm \
    --user $(id -u):$(id -g) \
    --volume $(pwd):/workdir \
    --workdir /workdir \
    quay.io/vgteam/vg:v1.75.1 \
    vg paths -x hprc-v1.1-mc-grch38.d9.gbz \
    -L --paths-by GRCh38 > hprc-v1.1-mc-grch38.d9.paths
# As per best practices, remove decoys, unplaced/unlocalized contigs,
# and other non-primary paths unnecessary for pangenome graph alignment.
grep -v _decoy hprc-v1.1-mc-grch38.d9.paths \
    | grep -v _random \
    | grep -v chrUn_ \
    | grep -v chrEBV \
    | grep -v chrM \
    | grep -v chain_ > hprc-v1.1-mc-grch38.d9.paths.sub

This produces the following files:

• hprc-v1.1-mc-grch38.d9.gbz — GBZ graph
• hprc-v1.1-mc-grch38.d9.autoindex.1.75.1.dist — distance index
• hprc-v1.1-mc-grch38.d9.autoindex.1.75.1.shortread.withzip.min — minimizer index
• hprc-v1.1-mc-grch38.d9.autoindex.1.75.1.shortread.zipcodes — zipcodes
• hprc-v1.1-mc-grch38.d9.paths.sub — filtered reference paths

Quick Start

Before running giraffe, ensure you have generated the required index files. Refer to the index generation section above for instructions.

# This command assumes all the inputs are in the current working directory and all the outputs go to the same place.
docker run --rm --gpus all --volume $(pwd):/workdir --volume $(pwd):/outputdir \
    --workdir /workdir \
    nvcr.io/nvidia/clara/clara-parabricks:4.8.0-1 \
    pbrun giraffe --read-group "sample_rg1" \
    --sample "sample-name" --read-group-library "library" \
    --read-group-platform "platform" --read-group-pu "pu" \
    --gbz-name /workdir/hprc-v1.1-mc-grch38.d9.gbz \
    --dist-name /workdir/hprc-v1.1-mc-grch38.d9.autoindex.1.75.1.dist \
    --minimizer-name /workdir/hprc-v1.1-mc-grch38.d9.autoindex.1.75.1.shortread.withzip.min \
    --zipcodes-name /workdir/hprc-v1.1-mc-grch38.d9.autoindex.1.75.1.shortread.zipcodes \
    --ref-paths /workdir/hprc-v1.1-mc-grch38.d9.paths.sub \
    --in-fq /workdir/${INPUT_FASTQ_1} /workdir/${INPUT_FASTQ_2} \
    --out-bam /outputdir/${OUTPUT_BAM}

CRAM Output

Giraffe writes CRAM instead of BAM when --out-bam is given a path ending in .cram; there is no separate output option. CRAM encoding needs a linear reference, so --ref is required in this case. The reference must match the graph passed to --gbz-name and the paths listed in --ref-paths; refer to Extracting the Reference FASTA for how to produce it from the graph.

• --run-ref-verification runs a pre-flight check of the reference (.fai plus the sequences themselves) against the GBZ graph. It requires the FASTA index, which can be created with samtools faidx. Without this option, Giraffe warns that an inconsistent reference may produce an invalid CRAM file, so enabling it is recommended whenever a reference is used for the first time.

• --align-only cannot be combined with CRAM output, because CRAM is written by the coordinate-sorting stage that --align-only skips.

• --gpuwrite uses one GPU to accelerate CRAM compression. A .crai index is written alongside the coordinate-sorted CRAM output.

# This command assumes all the inputs are in the current working directory and all the outputs go to the same place.
docker run --rm --gpus all --volume $(pwd):/workdir --volume $(pwd):/outputdir \
    --workdir /workdir \
    nvcr.io/nvidia/clara/clara-parabricks:4.8.0-1 \
    pbrun giraffe --read-group "sample_rg1" \
    --sample "sample-name" --read-group-library "library" \
    --read-group-platform "platform" --read-group-pu "pu" \
    --gbz-name /workdir/hprc-v1.1-mc-grch38.d9.gbz \
    --dist-name /workdir/hprc-v1.1-mc-grch38.d9.autoindex.1.75.1.dist \
    --minimizer-name /workdir/hprc-v1.1-mc-grch38.d9.autoindex.1.75.1.shortread.withzip.min \
    --zipcodes-name /workdir/hprc-v1.1-mc-grch38.d9.autoindex.1.75.1.shortread.zipcodes \
    --ref-paths /workdir/hprc-v1.1-mc-grch38.d9.paths.sub \
    --ref /workdir/hprc-v1.1-mc-grch38.d9.fa \
    --run-ref-verification \
    --in-fq /workdir/${INPUT_FASTQ_1} /workdir/${INPUT_FASTQ_2} \
    --out-bam /outputdir/${OUTPUT_CRAM}

System Requirements and Useful Options for Performance

To ensure optimal performance with VG Giraffe, consider the following system requirements based on your GPU configuration. Giraffe uses 16 CPU threads per GPU by default (--num-cpu-threads-per-gpu):

• A 2 GPU system should have at least 100GB CPU RAM and at least 32 CPU threads.

• A 4 GPU system should have at least 200GB CPU RAM and at least 64 CPU threads.

• An 8 GPU system should have at least 400GB CPU RAM and at least 128 CPU threads.

• For GPUs with less than 22 GB of device memory, use --low-memory.

Auto Mode

By default, --nstreams is set to auto, which enables auto mode. In this mode, Giraffe automatically configures the number of CUDA streams, batch size, and GPU acceleration options based on GPU device memory. Auto mode is designed to provide sensible defaults, but may still require further optimization for each specific system. The following table summarizes the auto mode configuration based on the total device memory of the smallest GPU, in MiB, as reported by nvidia-smi --query-gpu=memory.total. Each range includes its lower bound and excludes its upper bound:

Auto mode also takes total host memory into account:

• With less than 80 GiB of host memory, the batch size is reduced to 5000 and the work queue capacity to 10.
• --minimizers-gpu is disabled if host memory is insufficient for the total number of streams across all GPUs.

Note that --minimizers-gpu is only enabled for single-end (SE) reads. For paired-end (PE) reads, the number of streams and batch size are configured as shown above, but --minimizers-gpu is always disabled.

For best performance, auto mode can be overridden by explicitly setting --nstreams and other options. The following configurations are, for example, recommended for L4, H100 and RTX PRO 6000 Blackwell Server Edition GPUs:

• L4 (16 GB): --batch-size 5000 --nstreams 2

• H100 (80 GB): --nstreams 5 --num-cpu-threads-per-gpu 24 --minimizers-gpu

• RTX PRO 6000 Blackwell Server Edition (96 GB): --nstreams 6 --num-cpu-threads-per-gpu 24 --minimizers-gpu

Note: While a fixed base memory allocation exists per device, the number of streams and batch size are the primary factors affecting total device memory consumption.

Using Giraffe in Variant Calling Workflows

To use Giraffe-aligned BAM files for variant calling, you need to extract the appropriate reference file from the Giraffe index files.

Extracting the Reference FASTA

Run the following commands from the directory containing the Giraffe index files:

# Extract the sequences corresponding to the list of paths to a FASTA file
docker run --rm \
    --user $(id -u):$(id -g) \
    --volume $(pwd):/workdir \
    --workdir /workdir \
    quay.io/vgteam/vg:v1.75.1 \
    vg paths -x hprc-v1.1-mc-grch38.d9.gbz \
        -p hprc-v1.1-mc-grch38.d9.paths.sub \
        -F > hprc-v1.1-mc-grch38.d9.fa
# Index the FASTA file
docker run --rm \
    --user $(id -u):$(id -g) \
    --volume $(pwd):/workdir \
    --workdir /workdir \
    quay.io/biocontainers/samtools:1.17--hd87286a_2 \
    samtools faidx hprc-v1.1-mc-grch38.d9.fa

These commands generate a FASTA file (hprc-v1.1-mc-grch38.d9.fa) and the corresponding index (hprc-v1.1-mc-grch38.d9.fa.fai).

The FASTA can be used as the reference for variant calling, for CRAM output, and for BQSR (bqsr). We can now run Giraffe to obtain the aligned BAM as follows:

# This command assumes all the inputs are in the current working directory and all the outputs go to the same place.
docker run --rm --gpus all --volume $(pwd):/workdir --volume $(pwd):/outputdir \
    --workdir /workdir \
    nvcr.io/nvidia/clara/clara-parabricks:4.8.0-1 \
    pbrun giraffe --read-group "sample_rg1" \
    --sample "sample-name" --read-group-library "library" \
    --read-group-platform "platform" --read-group-pu "pu" \
    --gbz-name /workdir/hprc-v1.1-mc-grch38.d9.gbz \
    --dist-name /workdir/hprc-v1.1-mc-grch38.d9.autoindex.1.75.1.dist \
    --minimizer-name /workdir/hprc-v1.1-mc-grch38.d9.autoindex.1.75.1.shortread.withzip.min \
    --zipcodes-name /workdir/hprc-v1.1-mc-grch38.d9.autoindex.1.75.1.shortread.zipcodes \
    --ref-paths /workdir/hprc-v1.1-mc-grch38.d9.paths.sub \
    --in-fq /workdir/${INPUT_FASTQ_1} /workdir/${INPUT_FASTQ_2} \
    --out-bam /outputdir/${OUTPUT_BAM}

Once you have the Giraffe-aligned BAM file and the extracted reference FASTA, you can proceed with variant calling using HaplotypeCaller, Deepvariant, Pangenome_aware_deepvariant, or the end-to-end pangenome_germline pipeline.

# Haplotype Caller
# This command assumes all the inputs are in the current working directory and all the outputs go to the same place.
docker run --rm --gpus all --volume $(pwd):/workdir --volume $(pwd):/outputdir \
    --workdir /workdir \
    nvcr.io/nvidia/clara/clara-parabricks:4.8.0-1 \
    pbrun haplotypecaller \
    --ref /workdir/hprc-v1.1-mc-grch38.d9.fa \
    --in-bam /workdir/${INPUT_BAM} \
    --out-variants /outputdir/${OUTPUT_VCF}
# Deepvariant
# This command assumes all the inputs are in the current working directory and all the outputs go to the same place.
docker run --rm --gpus all --volume $(pwd):/workdir --volume $(pwd):/outputdir \
    --workdir /workdir \
    nvcr.io/nvidia/clara/clara-parabricks:4.8.0-1 \
    pbrun deepvariant \
    --ref /workdir/hprc-v1.1-mc-grch38.d9.fa \
    --in-bam /workdir/${INPUT_BAM} \
    --out-variants /outputdir/${OUTPUT_VCF}
# Pangenome_aware_deepvariant
# This command assumes all the inputs are in the current working directory and all the outputs go to the same place.
docker run --rm --gpus all --volume $(pwd):/workdir --volume $(pwd):/outputdir \
    --workdir /workdir \
    nvcr.io/nvidia/clara/clara-parabricks:4.8.0-1 \
    pbrun pangenome_aware_deepvariant \
    --ref /workdir/hprc-v1.1-mc-grch38.d9.fa \
    --pangenome /workdir/hprc-v1.1-mc-grch38.d9.gbz \
    --in-bam /workdir/${INPUT_BAM} \
    --out-variants /outputdir/${OUTPUT_VCF}

For more detailed instructions on variant calling, refer to the tool-specific documentation (haplotypecaller, deepvariant, pangenome_aware_deepvariant, pangenome_germline).

Using Giraffe with Haplotype Sampling

Giraffe’s haplotype sampling functionality, activated using arguments --haplotype-name and --kff-name, was introduced to significantly enhance alignment accuracy by tailoring the reference graph to the specific genetic profile of a sample. This process begins by analyzing sequencing reads with a kmer counter to identify patterns of kmer presence and frequency. Using this information, Giraffe sub-samples the GBWT (using the original .hapl and .gbz files) to select haplotypes that best represent the sample, creating a customized graph. From this tailored graph, Giraffe also generates new index files (.dist, .min, and .zipcodes) that are optimized for the sample to be analyzed.

These steps can be performed using the baseline VG container for graph customization and index generation, followed by Parabricks’ accelerated Giraffe for high-performance alignment, as demonstrated below.

The required .gbz file can be downloaded as follows:

aws s3 cp \
    s3://human-pangenomics/pangenomes/freeze/freeze1/minigraph-cactus/hprc-v1.1-mc-grch38/hprc-v1.1-mc-grch38.gbz \
    . \
    --no-sign-request

The .hapl file published alongside this graph predates the current haplotype information format and is rejected with a version error, so it has to be regenerated locally. vg haplotypes builds it from a distance index and an r-index, which it locates by file name, so both must sit next to the .gbz file and share its base name:

# Build the distance index
docker run --rm \
    --user $(id -u):$(id -g) \
    --volume $(pwd):/workdir \
    --workdir /workdir \
    quay.io/vgteam/vg:v1.75.1 \
    vg index -j hprc-v1.1-mc-grch38.dist \
        --no-nested-distance hprc-v1.1-mc-grch38.gbz
# Build the r-index
docker run --rm \
    --user $(id -u):$(id -g) \
    --volume $(pwd):/workdir \
    --workdir /workdir \
    quay.io/vgteam/vg:v1.75.1 \
    vg gbwt -p --num-threads 16 \
        -r hprc-v1.1-mc-grch38.ri \
        -Z hprc-v1.1-mc-grch38.gbz
# Build the haplotype information
docker run --rm \
    --user $(id -u):$(id -g) \
    --volume $(pwd):/workdir \
    --workdir /workdir \
    quay.io/vgteam/vg:v1.75.1 \
    vg haplotypes -v 2 -t 16 \
        -H hprc-v1.1-mc-grch38.hapl \
        hprc-v1.1-mc-grch38.gbz

These three commands write .dist, .ri, and .hapl files rather than a graph, so they do not affect the GBZ version of the input.

You will also need to extract and filter reference paths for this graph:

docker run --rm \
    --user $(id -u):$(id -g) \
    --volume $(pwd):/workdir \
    --workdir /workdir \
    quay.io/vgteam/vg:v1.75.1 \
    vg paths -x hprc-v1.1-mc-grch38.gbz \
    -L --paths-by GRCh38 > hprc-v1.1-mc-grch38.paths
grep -v _decoy hprc-v1.1-mc-grch38.paths \
    | grep -v _random \
    | grep -v chrUn_ \
    | grep -v chrEBV \
    | grep -v chrM \
    | grep -v chain_ > hprc-v1.1-mc-grch38.paths.sub

# Run KMC on the input reads to obtain the .kff file
mkdir kmc_tmp_dir
cat > input.fq.paths <<- EOM
${INPUT_FASTQ_1}
${INPUT_FASTQ_2}
EOM
docker run --rm --volume $(pwd):/workdir \
    --workdir /workdir \
    quay.io/biocontainers/kmc:3.2.4--haf24da9_3 \
    kmc \
        -k29 \
        -m128 \
        -okff \
        -t64 \
        @input.fq.paths \
        input.fq.distr kmc_tmp_dir
# Compute the sampled .gbz file using the baseline container
docker run --rm --volume $(pwd):/workdir \
    --workdir /workdir \
    quay.io/vgteam/vg:v1.75.1 \
    vg haplotypes \
        -v 2 -t 64 \
        --include-reference \
        --diploid-sampling \
        -i hprc-v1.1-mc-grch38.hapl \
        -k input.fq.distr.kff \
        -g hprc-v1.1-mc-grch38.sampled.gbz \
        hprc-v1.1-mc-grch38.gbz

vg haplotypes always writes the sampled graph as GBZ version 2, and has no option to write version 1, so the sampled graph has to be converted before Parabricks Giraffe can read it. Refer to Supported GBZ version for background:

# Convert the sampled graph to GBZ version 1
docker run --rm --volume $(pwd):/workdir \
    --workdir /workdir \
    --user $(id -u):$(id -g) \
    quay.io/vgteam/vg:v1.75.1 \
    vg gbwt -Z hprc-v1.1-mc-grch38.sampled.gbz \
        -g hprc-v1.1-mc-grch38.sampled.v1.gbz --gbz-v1

The index files are then generated from the converted graph, and the converted graph is also the one passed to --gbz-name:

# Generate index files from the sampled graph using autoindex
docker run --rm --volume $(pwd):/workdir \
    --workdir /workdir \
    --user $(id -u):$(id -g) \
    quay.io/vgteam/vg:v1.75.1 \
    vg autoindex \
        -p hprc-v1.1-mc-grch38.sampled.autoindex.1.75.1 \
        -G hprc-v1.1-mc-grch38.sampled.v1.gbz \
        -w giraffe
# Align the reads to the sampled graph using Parabricks Giraffe
# This command assumes all the inputs are in the current working directory and all the outputs go to the same place.
docker run --rm --gpus all --volume $(pwd):/workdir --volume $(pwd):/outputdir \
    --workdir /workdir \
    nvcr.io/nvidia/clara/clara-parabricks:4.8.0-1 \
    pbrun giraffe --read-group "sample_rg1" \
    --sample "sample-name" --read-group-library "library" \
    --read-group-platform "platform" --read-group-pu "pu" \
    --gbz-name hprc-v1.1-mc-grch38.sampled.v1.gbz \
    --dist-name hprc-v1.1-mc-grch38.sampled.autoindex.1.75.1.dist \
    --minimizer-name hprc-v1.1-mc-grch38.sampled.autoindex.1.75.1.shortread.withzip.min \
    --zipcodes-name hprc-v1.1-mc-grch38.sampled.autoindex.1.75.1.shortread.zipcodes \
    --ref-paths hprc-v1.1-mc-grch38.paths.sub \
    --in-fq ${INPUT_FASTQ_1} ${INPUT_FASTQ_2} \
    --out-bam /outputdir/${OUTPUT_BAM}

Aligning Ultima CRAM Input

Giraffe can align or re-align reads directly from a single-end CRAM file with --in-se-cram. Reads are always taken in sequencing orientation, so input that is already aligned is re-aligned from scratch. Paired-end CRAM input is not supported.

• --in-cram-ref supplies the reference FASTA used to decode the input. It is required only when the input contains aligned reads.

• --cram-read-threads sets the number of threads used to decode the input (default: 4).

Ultima Genomics data carries per-base tp and t0 auxiliary tags. Passing --ultima-tags carries these tags from the --in-se-cram input through to the output. This is supported for single-end input only and requires --in-se-cram.

Because alignment can reverse-complement a read, the order of the per-base values matters. Giraffe determines the order of the incoming values from the @RG lines of the input header:

• tp:reference — the values follow the SEQ field, so for reverse-strand records they run opposite to the sequencing direction.

• tp:synthesis — the values are in sequencing direction regardless of strand.

If the input contains aligned reads and no @RG line declares either tp:reference or tp:synthesis, Giraffe exits with an error, because the order of the values cannot be determined. Add the declaration to the read group before re-aligning such a file.

The order written to the output is controlled by --tp-output-order:

• reference (default) — the values match the output SEQ, and are reversed for reverse-strand alignments.

• synthesis — the values are kept in sequencing direction regardless of strand.

# This command assumes all the inputs are in the current working directory and all the outputs go to the same place.
docker run --rm --gpus all --volume $(pwd):/workdir --volume $(pwd):/outputdir \
    --workdir /workdir \
    nvcr.io/nvidia/clara/clara-parabricks:4.8.0-1 \
    pbrun giraffe --read-group "sample_rg1" \
    --sample "sample-name" --read-group-library "library" \
    --read-group-platform "platform" --read-group-pu "pu" \
    --gbz-name /workdir/hprc-v1.1-mc-grch38.d9.gbz \
    --dist-name /workdir/hprc-v1.1-mc-grch38.d9.autoindex.1.75.1.dist \
    --minimizer-name /workdir/hprc-v1.1-mc-grch38.d9.autoindex.1.75.1.shortread.withzip.min \
    --zipcodes-name /workdir/hprc-v1.1-mc-grch38.d9.autoindex.1.75.1.shortread.zipcodes \
    --ref-paths /workdir/hprc-v1.1-mc-grch38.d9.paths.sub \
    --in-se-cram /workdir/${INPUT_CRAM} \
    --in-cram-ref /workdir/${INPUT_CRAM_REFERENCE} \
    --ultima-tags \
    --tp-output-order reference \
    --out-bam /outputdir/${OUTPUT_BAM}

Note that --ultima-tags and --tp-output-order are available on pbrun giraffe only. The pangenome_germline pipeline accepts CRAM input but passes the alignments to a variant caller that has no use for the per-base tags.

Compatible CPU-based vg giraffe and GATK4 Commands

The commands below are the vg-1.75.1 and GATK4 counterpart of the Parabricks command above. The output from these commands will be identical to the output from the above command. Refer to the Output Comparison page for comparing the results. Note that the baseline needs a few upstream follow-up commits on top of the v1.75.1 tag; refer to Source of Mismatches for the build recipe.

The index files used below are generated in the index generation section.

# Run giraffe and pipe the output to create a sorted BAM.
$ vg giraffe \
    -t 16 \
    -Z /workdir/hprc-v1.1-mc-grch38.d9.gbz \
    -d /workdir/hprc-v1.1-mc-grch38.d9.autoindex.1.75.1.dist \
    -m /workdir/hprc-v1.1-mc-grch38.d9.autoindex.1.75.1.shortread.withzip.min \
    -z /workdir/hprc-v1.1-mc-grch38.d9.autoindex.1.75.1.shortread.zipcodes \
    --ref-paths /workdir/hprc-v1.1-mc-grch38.d9.paths.sub \
    -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
# Mark duplicates.
$ gatk MarkDuplicates \
    -I cpu.bam \
    -O cpu.markdup.bam \
    -M metrics.txt

Source of Mismatches

When comparing output with the CPU counterpart the following can be sources of small differences.

• Baseline VG Container

  Alignment output is validated against a CPU vg binary built from the upstream v1.75.1 tag with a few upstream follow-up commits cherry-picked on top. All of them were made after the release, so the unmodified v1.75.1 container produces different output: it cannot report comparable supplementary alignments, and its Dozeu implementation predates the two fixes that Parabricks matches. The commits below are the upstream fixes for both.

# Clone the repo (full history needed for cherry-pick)
git clone https://github.com/vgteam/vg.git
cd vg
# Checkout v1.75.1 tag and create a patch branch
git checkout v1.75.1
git checkout -b v1.75.1-suppl-fix-dozeu-fix
# The elfutils submodule URL in the tag no longer resolves
git config submodule.deps/elfutils.url https://github.com/vgteam/elfutils.git
git submodule update --init --recursive
# Supplementary-alignment fixes
git cherry-pick 4f873b0ccb37789aa53281e7b918bb7fcc5e4e57   # gate the per-read log message
git cherry-pick a34eb7d2e494166a7b23c1aa7f2d5949e37b109e   # first round of fixes
git cherry-pick d7227c4454a13c78b11d26499583dbc09b0468db   # follow-up fix
git cherry-pick e8389c3955505b11dfff7dce73e7e654ca2c331d   # std::move fix
# Dozeu fixes
git cherry-pick 05d80785dd69561b232c3b099391d9ca5bb16f0a
git cherry-pick ea9a9150fa4903e43040f5c3ff298a3d967e856e
git submodule update --init --recursive

The patched checkout is then built with upstream vg’s own Dockerfile, which carries the full dependency list and build steps:

# Record the version string compiled into the binary
make version
# Build the container from the patched checkout
docker build --build-arg THREADS=16 -t vg:v1.75.1-suppl-fix-dozeu-fix .
# Verify the resulting binary
docker run --rm vg:v1.75.1-suppl-fix-dozeu-fix vg version

THREADS controls build parallelism and is clamped to the number of available cores. Skipping make version is possible, but then the version must be supplied with --build-arg VG_GIT_VERSION=<value> or the binary reports its version as unknown.

• Unmapped reads

  • Parabricks giraffe sorts unmapped reads slightly differently than baseline GATK SortSam. Unmapped reads can be filtered with samtools by running samtools view -F 4.

Options