-
Notifications
You must be signed in to change notification settings - Fork 0
Home
For the quickstart, see the main README
Certain workflow steps are optional, configurable via the --steps parameter, which takes a comma separated string, 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'
The optional steps are listed as a string below, though note that some steps depend on others, so users will be prompted to provide any missing requirements:
preprocess,index,align,merge,deduplicate,profile_pmd,graph_genotype,reads_genotype,variant_call
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 samplesheet or reference file has been provided when required.
Reference processing tasks depend on the input and other requested steps.
If indexing is requested, for filtered graphs, sample-type appropriate '.dist' and '.min' indexes are computed with vg tools (i.e., more than one '.min' index will be made if users run a mix of modern and aDNA samples). For unfiltered graphs, haplotype indexes.
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 more about snarls here.
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 dedicated deduplication step after mapping and sample merging.
The alignment strategy depends on the reference type and sample type.
Reads are mapped to pangenome graphs with vg giraffe. 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 defaults (modelled on bwa mem). 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.
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.
Alignment statistics are reported with vg stats both before and after alignment filtering.
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.
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.
Duplicate reads are identified and removed using Picard MarkDuplicates. For ancient samples, deduplication considers 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. This behaviour mimics that of the popular aDNA tool DeDup, read more here. 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.
Post-mortem DNA damage patterns are assessed for ancient samples using DamageProfiler. This step is skipped for modern samples.
Genotyping is only available for graph input modes. Graph genotyping uses vg deconstruct to identify variants represented as snarls (i.e., bubbles) 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 note that they are 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 operates on deduplicated BAM files and currently supports two callers, which can be run independently or together: freebayes and DeepVariant. Note that DeepVariant is not run by default, and is only recommended for human data.
Freebayes can be run in single-threaded (default) or parallelised mode, though the latter is currently experimental. 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). 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.
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.
Users can optionally produce details reports of resource usage, run time, and other information, by providing the parameter --tracing on the command line (or tracing: true if supplied in a params.yml)
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.
- INDEX_FILTERED_GRAPH
- INDEX_UNFILTERED_GRAPH
- PROCESS_GRAPH
- COMPUTE_SNARLS
- VG_DECONSTRUCT
- FASTP
- FASTQC
- PANGENOME_MAP
- VG_GENOTYPE
- GAM_TO_TAGGED_SORTED_BAM
- MERGE_BAMS
- DEDUPLICATE
- PROFILE_PMD
- FREEBAYES
- DEEPVARIANT
- PROCESS_DEEPVARIANT
Common questions will be added here!