Introduction

nf-core/genomicrelatedness is a bioinformatics pipeline for estimating genomic relatedness from low-coverage whole-genome sequencing (lcWGS) data. It performs read mapping, optional base quality score recalibration, variant calling with GATK and BCFtools, and downstream relatedness estimation. For many non-model organisms, no high-confidence variant set is available. The pipeline provides an automated multi-round bootstrapping workflow to generate one. The resulting standardized outputs include genotype likelihood-based variant calls, filtered VCF files, and relatedness estimates from several independent algorithms, enabling robust inference even from very sparse sequencing data.

Quick Start

1. Install nextflow

2. Install a container environment for full pipeline reproducibility, e.g. Docker, Singularity, Podman, Shifter or Charliecloud (please only use Conda as a last resort; see docs)

3. Download the pipeline and test it on a minimal dataset with a single command:

   nextflow run nf-core/genomicrelatedness -profile test,<docker/singularity/podman/shifter/charliecloud/conda/institute>
   

   Many Institutes have custom config files to run nextflow pipelines available. You can check nf-core/configs to see if one exists for your Institute. Then you can simply use -profile <institute> in your command. This will enable either docker or singularity and set the appropriate execution settings for your local compute environment, e.g. HPC cluster.

4. Start analysing your own data!

   nextflow run nf-core/genomicrelatedness -profile <docker/singularity/podman/conda/institute> --input 'samplesheet.csv' --fasta '<REFGENOME>.fasta' --outdir <OUTDIR>
   

5. Once your run has completed successfully, you'll find an overview of the run in the MultiQC report located at ./outdir/MultiQC/multiqc_report.html

See usage docs for all of the available options when running the pipeline and output docs for performed analyses and results files.

Modifications to the default pipeline are easily made using various options as described in the documentation.

Pipeline summary

The pipeline consists of the following four main sections that perform the major processing steps:

1. Preprocessing section: Prepares the input files for downstream analyses.

   • Reference genome preparation (--prepare_genome): This step prepares the input reference genome and generates the relevant files for subsequent processing steps. Key processes include:

     • Indexing of the reference genome with BWA-mem2 -index to generate alignment index files.

     • Indexing of the reference genome with samtools -faidx to generate FASTA index (.fai)

     • Creating sequence dictionary with the GATK4 -createsequencedictionary to generate sequence dictionary file (.dict)

   • Interval preparation (--prepare_intervals): This step splits the indexed reference genome into intervals for computational more efficient downstrean analyses.

     • Build and split intervals of the reference genome with [gawk -build_intervals] and [split_intervals] to generate interval file (.bed)

   • Preprocessing of raw sequencing reads (--preprocess): This step performs all essential steps to provide aligned sequence reads and quality metrics.

     • Input parsing & metadata setup: Reads a CSV samplesheet describing the input FASTQ or SPRING files with raw sequencing reads and their read-group information.

     • Raw read quality control: performs quality control, trimming, filtering and merging of paired reads with fastp, generating trimmed and filtered paired reads (.merged.fastq.gz)

     • Read mapping: Aligns the paired reads to the reference genome with BWA-mem2 -mem, sorts reads with samtools -sort, adds read group information with GATK4 -addorreplacereadgroups.

     • Merge files stemming from the same sample with samtools -merge.

     • Remove duplicates with GATK4 -markduplicates.

     • Indexing of the mapped reads with samtools -index, providing sorted and indexed compressed alignment files with proper read group annotations (.cram) and corresponding index files (.crai).

     • Quality metrics of the sequencing data are compiled with MultiQC and include overviews of library complexity analysed with preseq -c_curve and preseq -lcextrap, read mapping statics over the various preprocessing steps with samtools -stats, and genome coverage using mosdepth.

2. Bootstrapping section: corrects systematic errors introduced during sequencing by adjusting base quality scores.

   • Variant bootstrapping (--bootstrap_variant_set): (OPTIONAL) If no known variant set is available, this step automatically generates one via bootsrapping. It iteratively refines and stabilises the set of high-confidence variants for downstream use.

     • Variant calling performs joint variant discovery for all samples with GATK4 -HaplotypeCaller, GATK4 -GenomicsDBImport, GATK4 -GenotypeGVCFs, GATK4 -mergevcfs.

     • Hard filtering of variants to create a high confidence reference variant set (.vcf) with GATK4 -VariantFiltration (Qual >=100, QD < 2.0; MQ < 35.0; FS >60; HaplotypeScore > 13.0; MQRankSum < -12.5; ReadPosRankSum < -8.0) and GATK4 -SelectVariants.

   • BQSR (--base-quality-score-recalibration): uses the produced or a provided variant set to adjust quality scores in alignment files:

     • BQSR with GATK4 -BaseRecalibrator, GATK4 -GatherBQSRReport, GATK4 -ApplyBQSR, samtools -merge, and samtools -index to produce recalibrated alignment files (.cram).

     • BQSR diagnostics with GATK4 -AnalyzeCovariates to evaluate the effectiveness of recalibration.

     • Iteration (DEFAULT: 2 rounds) of recalibration.

3. Variant calling section: calls variants with different algorithms to produce genotype likelihoods for all samples. The default is that both subworkflows run in parallel, but can be adjusted to one or the other.

   • Variant calling with GATK (--call_variants-gatk) performs joint variant discovery for all samples with GATK4 -HaplotypeCaller , GATK4 -GenomicsDBImport, GATK4 -Genoty