nf-core/raredisease

Edit

Introduction

This document describes the output produced by the pipeline. Most of the plots are taken from the MultiQC report, which summarises results at the end of the pipeline.

The directories listed below will be created in the results directory after the pipeline has finished. All paths are relative to the top-level results directory.

Pipeline overview

The pipeline is built using Nextflow and processes data using the following steps:

• Alignment
  • Mapping
    • Bwa-mem2
    • Sentieon bwa mem
  • Duplicate marking
    • Picard’s MarkDuplicates
    • Sentieon dedup
• Quality control and reporting
  • Quality control
    • FastQC
    • Mosdepth
    • Picard tools
    • Qualimap
    • Sentieon WgsMetricsAlgo
    • TIDDIT’s cov and UCSC WigToBigWig
  • Reporting
    • MultiQC
• Variant calling - SNV
  • DeepVariant
  • Sentieon DNAscope
• Variant calling - SV
  • Manta
  • TIDDIT sv
  • GATK GermlineCNVCaller - CNV calling
  • SVDB merge
• Variant calling - repeat expansions
  • Expansion Hunter
  • Stranger
• Annotation - SNV
  • bcftools roh
  • vcfanno
  • CADD
  • VEP
  • UPD
  • Chromograph
• Annotation - SV
  • SVDB query
  • VEP
• Mitochondrial analysis
  • Alignment and variant calling
    • MT deletion script
  • Annotation:
    • HaploGrep2
    • vcfanno
    • CADD
    • VEP
    • HmtNote
• Rank variants and filtering
  • GENMOD
• Pipeline information

Alignment

Mapping

Bwa-mem2

Bwa-mem2 used to map the reads to a reference genome. The aligned reads are coordinate sorted with samtools sort. These files are treated as intermediates and are not placed in the output folder by default.

Sentieon bwa mem

Sentieon’s bwa mem is the software accelerated version of the bwa-mem algorithm. It is used to efficiently perform the alignment using BWA. Aligned reads are then coordinate sorted using Sentieon’s sort utility. These files are treated as intermediates and are not placed in the output folder by default. It is not the default aligner, but it can be chosen over bwamem2 by setting --aligner option to sentieon.

Duplicate marking

Picard’s MarkDuplicates

Picard MarkDuplicates is used for marking PCR duplicates that can occur during library amplification. This is essential as the presence of such duplicates results in false inflated coverages, which in turn can lead to overly-confident genotyping calls during variant calling. Only reads aligned by Bwa-mem2 are processed by this tool. By default, alignment files are published in bam format. If you would like to store cram files instead, set --save_mapped_as_cram to true.

Sentieon Dedup

Sentieon Dedup is the algorithm used by Sentieon’s driver to remove duplicate reads. Only reads aligned by Sentieon’s implementation of bwa are processed by this algorithm. By default, alignment files are published in bam format. If you would like to store cram files instead, set --save_mapped_as_cram to true.

Quality control and reporting

Quality control

FastQC

FastQC gives general quality metrics about your sequenced reads. It provides information about the quality score distribution across your reads, per base sequence content (%A/T/G/C), adapter contamination and overrepresented sequences. For further reading and documentation see the FastQC help pages.

Mosdepth

Mosdepth is used to report quality control metrics such as coverage, and GC content from alignment files. The global distribution file, generated by this program is passed to MultiQC for generating the following plots,

• Cumulative coverage distribution
• Coverage distribution
• Average coverage per contig

Picard tools

Picard’s CollectMutipleMetrics, CollectHsMetrics, and CollectWgsMetrics We use Picardtools’ CollectWgsMetrics and CollectHsMetrics utilities to calculate metrics about coverage and performance of WGS & WES experiments. In addition to those metrics, we use CollectMultipleMetrics to gather information about alignment summary, insert size, GC content etc., The metrics generated by these three utilites are passed along to MultiQC to generate several plots as well.

Qualimap

Qualimap also allows you to assess the alignment coverage. Qualimap results are used by MultiQC to generate the following plots.

• Coverage histogram
• Cumulative genome coverage
• Insert size histogram
• GC content distribution

Sention WgsMetricsAlgo

Sentieon’s WgsMetricsAlgo is the Sentieon’s equivalent of Picard’s CollectWgsMetrics.

TIDDIT’s cov and UCSC WigToBigWig

TIDDIT’s cov is used to analyse the read depth of a bam file and generates a coverage report in wig format. This file is later passed to UCSC WigToBigWig to convert the file into a bigwig.

Reporting

MultiQC

MultiQC is a visualization tool that generates a single HTML report summarising all samples in your project. Most of the pipeline QC results are visualised in the report and further statistics are available in the report data directory.

Results generated by MultiQC collate pipeline QC from supported tools e.g. FastQC. The pipeline has special steps which also allow the software versions to be reported in the MultiQC output for future traceability. For more information about how to use MultiQC reports, see http://multiqc.info.

Variant calling - SNV

DeepVariant

DeepVariant is a deep learning-based variant caller that takes aligned reads, produces pileup image tensors from them, classifies each tensor using a convolutional neural network and finally reports the results in a standard VCF or gVCF file. Variant calls generated by DeepVariant are joint genotyped using GLnexus, and then normalized using bcftools norm. Only the normalized vcfs are placed in the output folder by default.

NB: In case you are running the separate mitochondrial analysis, mitochondrial calls are filtered from the normalized vcfs before they are published using GATK SelectVariants.

Sentieon DNAscope

The pipeline performs variant calling using Sentieon DNAscope with a machine learning model. This approach identifies the candidate sites with a higher accuracy, and calculates genotypes for each sample at that site. These files are treated as intermediates and are not placed in the output folder by default. DNAscope is not run by default. To use DNAscope instead of DeepVariant, set --variant_caller to sentieon.

Variant calling - SV

Manta

Manta calls structural variants (SVs) and indels from mapped paired-end sequencing reads. It combines paired and split-read evidence during SV discovery and scoring to improve accuracy, but does not require split-reads or successful breakpoint assemblies to report a variant in cases where there is strong evidence otherwise. Output vcf files are treated as intermediates and are not placed in the output folder.

TIDDIT sv

TIDDIT’s sv is used to identify chromosomal rearrangements using sequencing data. TIDDIT identifies intra and inter-chromosomal translocations, deletions, tandem-duplications and inversions, using supplementary alignments as well as discordant pairs. TIDDIT searches for discordant reads and split reads (supplementary alignments). Output vcf files are treated as intermediates and are not placed in the output folder.

GATK GermlineCNVCaller - CNV calling

GATK GermlineCNVCaller is used to identify copy number variants in germline samples given their read counts and a model describing a sample’s ploidy. Output vcf files are treated as intermediates and are not placed in the output folder.

SVDB merge

SVDB merge is used to merge the variant calls from GATK’s GermlineCNVCaller (only if skip_cnv_calling is set to false), Manta, and TIDDIT. Output files are published in the output folder.

Variant calling - repeat expansions

Expansion Hunter

Expansion Hunter aims to estimate sizes of repeat sequences by performing a targeted search through alignments that span, flank, and are fully contained in each repeat.

Stranger

Stranger annotates output files from Expansion Hunter with the pathologic implications of the repeat sizes.

Annotation - SNV

bcftools roh

bcftools roh is a program for detecting runs of homo/autozygosity.from only bi-allelic sites. The output files are not published in the output folder, and is passed to vcfanno for further annotation.

vcfanno

vcfanno allows you to quickly annotate your VCF with any number of INFO fields from any number of VCFs. It uses a simple configuration file to allow the user to specify the source annotation files and fields and how they will be added to the info of the query VCF. Values are pulled by name from the INFO field with special-cases of ID and FILTER to pull from those VCF columns. The output files are not published in the output folder, and is passed to CADD and/or VEP for further annotation.

We recommend using vcfanno to annotate SNVs with precomputed CADD scores (files can be downloaded from here).

CADD

CADD is a tool for scoring the deleteriousness of single nucleotide variants as well as insertion/deletions variants in the human genome. In nf-core/raredisease, SNVs can be annotated with precomputed CADD scores using vcfanno. However, for small indels they will be calculated on the fly by CADD. The output files are not published in the output folder, and is passed to VEP for further annotation.

VEP

VEP determines the effect of your variants on genes, transcripts, and protein sequence, as well as regulatory regions. We recommend annotating with the following plugins:

• LoFtool
• pLI
• SpliceAI
• MaxEntScan

Based on VEP annotations, custom scripts used by the pipeline further annotate each record with the most severe consequence, and pli scores.

NB: Output files described below include mitochondrial annotations only if —skip_mt_analysis is set to true.

UPD

UPD calls regions of uniparental disomy from germline exome/wgs trios. Output from UPD is passed to chromograph for making plots.

Chromograph

Chromograph is a python package to create PNG images from genetics data such as BED and WIG files.

Annotation - SV

SVDB query

SVDB query allows you to quickly annotate your VCF with data from one or more structural variant databases. The output files are not published in the output folder, and is passed to vep for further annotation.

VEP

VEP determines the effect of your variants on genes, transcripts, and protein sequence, as well as regulatory regions. We recommend annotating with pLI plugin, along with any other custom plugins you may want too use. Based on VEP annotations, custom scripts used by the pipeline further annotate each record with the most severe consequence, and pli scores.

Mitochondrial analysis

Mitochondrial analysis is run by default, to turn it off set --skip_mt_analysis to true.

Alignment and variant calling

Alignment and variant calling - GATK Mitochondrial short variant discovery pipeline The mitochondrial genome poses several challenges to the identification and understanding of somatic variants. The circularity of the mitochondrial genome means that the breakpoint in the reference genome is at an arbitrary position in the non-coding control region, creating a challenge in analyzing variation. Additionally, insertions of mitochondrial DNA into the nuclear genome (NuMTs) complicate the mapping of the mitochondrial genome and the distinction between NuMTs and the mitochondrial contig of interest. Lastly, mitochondrial variants often have very low heteroplasmy. Such low allele fraction (AF) variants can thus be mistaken for inherent sequencer noise.

The pipeline for mitochondrial variant discovery, using Mutect2, uses a high sensitivity to low AF and separate alignments using opposite genome breakpoints to allow for the tracing of lineages of rare mitochondrial variants.

MT deletion script

MT deletion script lists the fraction of mitochondrially aligning read pairs (per 1000) that appear discordant, as defined by an insert size of more than 1.2 kb (and less than 15 kb due to the circular nature of the genome) using samtools.

Annotation:

HaploGrep2

HaploGrep2 allows detecting artificial recombinants and missing variants as well as annotating rare and phantom mutations in mitochondria. Haplogrep generates a text report, which is published by default.

vcfanno

vcfanno allows you to quickly annotate your VCF with any number of INFO fields from any number of VCFs. It uses a simple conf file to allow the user to specify the source annotation files and fields and how they will be added to the info of the query VCF. Values are pulled by name from the INFO field with special-cases of ID and FILTER to pull from those VCF columns. The output files are not published in the output folder, and is passed to vep for further annotation.

We recommend using vcfanno to annotate SNVs with precomputed CADD scores (files can be downloaded from here).

CADD

CADD is a tool for scoring the deleteriousness of single nucleotide variants as well as insertion/deletions variants in the human genome. In nf-core/raredisease, SNVs can be annotated with precomputed CADD scores using vcfanno. However, for small indels they will be calculated on the fly by CADD. The output files are not published in the output folder, and is passed to VEP for further annotation.

Hmtnote

HmtNote annotates vcf containing human mitochondrial variants with HmtVar. It will run offline by default with a database within the container.

VEP

VEP determines the effect of your variants on genes, transcripts, and protein sequence, as well as regulatory regions.

Rank variants and filtering

GENMOD

GENMOD is a simple to use command line tool for annotating and analyzing genomic variations in the VCF file format. GENMOD can annotate genetic patterns of inheritance in vcf

Pipeline information

Nextflow provides excellent functionality for generating various reports relevant to the running and execution of the pipeline. This will allow you to troubleshoot errors with the running of the pipeline, and also provide you with other information such as launch commands, run times and resource usage.