Skip to content

Home

Jump to bottom
CormacKinsella edited this page Mar 23, 2026 · 20 revisions

Welcome to the grave wiki!

For the quickstart, see the main README

Configuration of workflow steps

Certain workflow steps are optional, configurable via the --steps parameter, e.g.:

# This command runs grave in the apptainer pixi environment, using test configuration, and executes only the read preprocessing & graph indexing steps

pixi run -e apptainer nextflow main.nf -profile apptainer,test -params-file tests/test_graph.yml --steps 'preprocess,index'

Optional steps are below, though note that some steps depend on others, so users will be prompted to provide the missing requirements:

preprocess,index,align,merge,deduplicate,profile_pmd,graph_genotype,reads_genotype,variant_call

Summary of workflow steps

Initialisation (not optional)

On startup, grave validates the user-provided parameters and samplesheet, which specifies sample identifiers, library identifiers, repeat numbers, sample type (ancient or modern), merge status, and paths to FASTQ files. Each row must represent a unique combination of sample, library, and repeat number. The workflow validates the requested steps and their dependencies (e.g., alignment requires both preprocessing and indexing to be requested), and checks that a reference file has been provided when required.

Reference processing (always runs, but some tasks are optional/automically run if required)

Reference processing tasks depend on the input and other requested steps.

If indexing is requested, for unfiltered graphs, haplotype indexes are built with vg tools using sample-type appropriate k-mer lengths (i.e., more than one '.hapl' index will be made if users run a mix of modern and aDNA samples). For filtered graphs, appropriate '.dist' and '.min' indexes are precomputed instead. FASTA inputs are standardised to a '.fasta' extension and BWA is used for indexing.

With any graph input, graph statistics are computed with vg stats, while with FASTA, samtools faidx is used.

If linear genomes are required for tasks downstream of graph mapping, reference paths are extracted in FASTA format using vg paths.

If genotyping steps are requested, snarls (bubble features of graphs that represent sites of variation) are computed with vg snarls.

Read preprocessing (optional via 'preprocess')

Read quality is assessed before and after filtering with FastQC. Adapter trimming, duplication estimation, length filtering, and (for paired-end ancient DNA) read merging with base correction are performed with fastp. Modern paired-end reads are not merged. Pre-merged reads are processed in single-end mode. Duplicates detected by the exact hash algorithm of FASTP are not removed by default, as there is a subsequent deduplication step after mapping and sample merging.

Alignment (optional via 'align')

The alignment strategy depends on the reference type and sample type.

Reads are mapped to pangenome graphs with vg giraffe. For unfiltered graphs, a sample-specific subgraph is first generated using a k-mer profile generated from the reads (computed with KMC). The subgraph is indexed and used for mapping that sample. For filtered graphs, the pre-built indexes are used directly, and reads are mapped to the complete graph. Ancient DNA samples are mapped with alignment parameters modelled on BWA ALN (--mismatch 3 --gap-open 11 --gap-extend 4), while modern samples use Giraffe's defaults (modelled on BWA MEM).

The resulting GAM alignments are filtered using vg filter, with entirely configurable parameters. By default, unmapped reads are discarded, alignments are filtered on MAPQ (>= 30) and also on minimum alignment identity (0.8), and end-clipping of ambiguous alignments is not done. Samples with zero passing alignments are flagged as failed and excluded from downstream steps.

Filtered GAM files are then projected onto the linear coordinates of the reference sample using vg surject to produce BAM files, with read group tags added via samtools addreplacerg and coordinate sorting via samtools sort. When the graph contains multiple reference assemblies, surjection is performed separately against each reference sample.

With linear reference inputs, merged aDNA reads are aligned with bwa aln and bwa samse, using aDNA-appropriate parameters (-l 16500 -n 0.01 -o 2). Modern samples are aligned with bwa mem using default parameters. In both cases, unmapped reads are excluded with samtools view, output is coordinate-sorted with samtools sort, and mapping statistics are collected with samtools flagstat.

Sample merging (optional via 'merge')

After alignment, library-level BAM files belonging to the same sample are merged into a single sample-level BAM using samtools merge. When only one library is present for a sample, the merge step is bypassed and the file is simply renamed.

Deduplication (optional via 'deduplicate')

Duplicate reads are identified and removed using Picard MarkDuplicates. For ancient samples, by default deduplication uses both the 5' and 3' mapping coordinates of reads (--FLOW_MODE --FLOW_USE_END_IN_UNPAIRED_READS), taking advantage of the known fragment endpoints in merged reads. For modern samples, only the 5' position is considered. Users can control whether duplicates are removed (yes by default) or merely flagged, and can configure the duplicate tagging policy. Deduplicated BAMs are indexed with samtools index. In multi-reference mode, merging and deduplication are performed independently for each surjection target.

aDNA Damage Profiling (optional via 'profile_pmd')

Post-mortem DNA damage patterns are assessed for ancient samples using DamageProfiler. This step is skipped entirely for modern samples.

Genotyping (optional via 'graph_genotype' & 'reads_genotype')

Genotyping is available only for graph input modes. Graph genotyping uses vg deconstruct to identify variants represented as snarls in the graph, relative to the reference paths, and produces a VCF without requiring read data.

Read-based genotyping uses vg pack to compute per-sample read support across the graph, then vg call to genotype snarls. When a sample has multiple library GAM files, these are concatenated prior to genotyping, but not deduplicated.

In both approaches, output VCFs are filtered with vcfbub, normalised with bcftools norm, and sorted with bcftools sort. In multi-reference mode, genotyping is performed independently relative to each reference assembly.

Variant Calling (optional via 'variant_call')

Variant calling operates on deduplicated BAM files and currently supports two callers, which can be run independently or together: FreeBayes and DeepVariant. Note that DeepVariant is off by default, and is only recommended for human data.

FreeBayes can be run in single-threaded or parallelised mode, though the latter is currently not recommended for large samples due to memory issues. Raw VCFs are normalised with bcftools norm and sorted.

For DeepVariant, in graph mode the pangenome-aware version of DeepVariant (run_pangenome_aware_deepvariant) is used, which takes both the reference FASTA and the GBZ graph as inputs. In linear mode, standard DeepVariant (run_deepvariant) is used. The model type is configurable (default: WGS). DeepVariant output VCFs are normalised with bcftools norm and sorted.

In multi-reference mode, both callers are run independently against each surjected reference assembly. An interactive HTML variant statistics report is produced by DeepVariant for each sample.

Version Reporting

Software versions for all tools used across the run are collected automatically via Nextflow's topic channel mechanism and written to a consolidated report file, ensuring full reproducibility.

Details of individual processes

NOTE: This wiki section is partly out of date as of grave version 2, but much of the detail has not changed. If a specific grave feature or behaviour is not clear, please make a GitHub Issue and it will be prioritised for further/revised documentation.

FAQ

Common questions will be added here!

Clone this wiki locally