nf-core/sarek

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:

• Directory Structure
• Preprocessing
  • Preparation of input files (FastQ or (u)BAM)
    • Clip and filter read length
    • Trim adapters
    • Split FastQ files
    • UMI consensus
  • Map to Reference
    • BWA
    • BWA-mem2
    • DragMap
    • Sentieon BWA mem
  • Mark Duplicates
    • GATK MarkDuplicates (Spark)
  • Sentieon LocusCollector and Dedup
  • Base Quality Score Recalibration
    • GATK BaseRecalibrator (Spark)
    • GATK ApplyBQSR (Spark)
  • CSV files
• Variant Calling
  • SNVs and small indels
    • bcftools
    • DeepVariant
    • FreeBayes
    • GATK HaplotypeCaller
      • GATK Germline Single Sample Variant Calling
      • GATK Joint Germline Variant Calling
    • GATK Mutect2
    • Sentieon DNAscope
      • Sentieon DNAscope joint germline variant calling
    • Sentieon Haplotyper
      • Sentieon Haplotyper joint germline variant calling
    • Strelka
    • Lofreq
  • Structural Variants
    • Indexcov
    • Manta
    • TIDDIT
  • Sample heterogeneity, ploidy and CNVs
    • ASCAT
    • CNVKit
    • Control-FREEC
  • Microsatellite instability (MSI)
    • MSIsensorPro
  • Concatenation
• Variant annotation
  • snpEff
  • VEP
  • BCFtools annotate
• Quality control and reporting
  • Quality control
    • FastQC
    • FastP
    • Mosdepth
    • NGSCheckMate
    • GATK MarkDuplicates reports
    • Sentieon Dedup reports
    • samtools stats
    • bcftools stats
    • VCFtools
    • snpEff reports
    • VEP reports
  • Reporting
    • MultiQC
  • Pipeline information
• Reference files

Directory Structure

The default directory structure is as follows

{outdir}
├── csv
├── multiqc
├── pipeline_info
├── preprocessing
│   ├── markduplicates
│       └── <sample>
│   ├── recal_table
│       └── <sample>
│   └── recalibrated
│       └── <sample>
├── reference
└── reports
    ├── <tool1>
    └── <tool2>
work/
.nextflow.log

Preprocessing

Sarek pre-processes raw FastQ files or unmapped BAM files, based on GATK best practices.

Preparation of input files (FastQ or (u)BAM)

FastP is a tool designed to provide all-in-one preprocessing for FastQ files and as such is used for trimming and splitting. By default, these files are not published. However, if publishing is enabled, please be aware that these files are only published once, meaning if trimming and splitting is enabled, then the resulting files will be sharded FastQ files with trimmed reads. If only one of them is enabled then the files contain either trimmed or split reads, respectively.

Clip and filter read length

FastP enables efficient clipping of reads from either the 5’ end (--clip_r1, --clip_r2) or the 3’ end (--three_prime_clip_r1, --three_prime_clip_r2). Additionally, FastP allows the filtering of reads based on insert size by specifying a minimum required length with the --length_required parameter (default: 15bp). It is recommended to optimize these parameters according to the specific characteristics of your data.

Trim adapters

FastP supports global trimming, which means it trims all reads in the front or the tail. This function is useful since sometimes you want to drop some cycles of a sequencing run. In the current implementation in Sarek --detect_adapter_for_pe is set by default which enables auto-detection of adapter sequences. For more information on how to fine-tune adapter trimming, take a look into the parameter docs.

The resulting files are intermediate and by default not kept in the final files delivered to users. Set --save_trimmed to enable publishing of the files in:

Split FastQ files

FastP supports splitting of one FastQ file into multiple files allowing parallel alignment of sharded FastQ file. To enable splitting, the number of reads per output can be specified. For more information, take a look into the parameter --split_fastqin the parameter docs.

These files are intermediate and by default not placed in the output-folder kept in the final files delivered to users. Set --save_split to enable publishing of these files to:

UMI consensus

Sarek can process UMI-reads, using fgbio tools.

These files are intermediate and by default not placed in the output-folder kept in the final files delivered to users. Set --save_split to enable publishing of these files to:

Map to Reference

BWA

BWA is a software package for mapping low-divergent sequences against a large reference genome. The aligned reads are then coordinate-sorted (or name-sorted if GATK MarkDuplicatesSpark is used for duplicate marking) with samtools.

BWA-mem2

BWA-mem2 is a software package for mapping low-divergent sequences against a large reference genome.The aligned reads are then coordinate-sorted (or name-sorted if GATK MarkDuplicatesSpark is used for duplicate marking) with samtools.

DragMap

DragMap is an open-source software implementation of the DRAGEN mapper, which the Illumina team created so that we would have an open-source way to produce the same results as their proprietary DRAGEN hardware. The aligned reads are then coordinate-sorted (or name-sorted if GATK MarkDuplicatesSpark is used for duplicate marking) with samtools.

These files are intermediate and by default not placed in the output-folder kept in the final files delivered to users. Set --save_mapped to enable publishing, furthermore add the flag save_output_as_bam for publishing in BAM format.

Sentieon BWA mem

Sentieon bwa mem is a subroutine for mapping low-divergent sequences against a large reference genome. It is part of the proprietary software package DNAseq from Sentieon.

The aligned reads are coordinate-sorted with Sentieon.

Mark Duplicates

During duplicate marking, read pairs that are likely to have originated from duplicates of the same original DNA fragments through some artificial processes are identified. These are considered to be non-independent observations, so all but a single read pair within each set of duplicates are marked, causing the marked pairs to be ignored by default during the variant discovery process.

For further reading and documentation see the data pre-processing for variant discovery from the GATK best practices.

GATK MarkDuplicates (Spark)

By default, Sarek will use GATK MarkDuplicates.

To use the corresponding spark implementation GATK MarkDuplicatesSpark, please specify --use_gatk_spark markduplicates. The resulting files are converted to CRAM with either samtools, when GATK MarkDuplicates is used, or, implicitly, by GATK MarkDuplicatesSpark.

The resulting CRAM files are delivered to the users.

Sentieon LocusCollector and Dedup

The subroutines LocusCollector and Dedup are part of Sentieon DNAseq packages with speedup versions of the standard GATK tools, and together those two subroutines correspond to GATK’s MarkDuplicates.

The subroutine LocusCollector collects read information that will be used for removing or tagging duplicate reads; its output is the score file indicating which reads are likely duplicates.

The subroutine Dedup marks or removes duplicate reads based on the score file supplied by LocusCollector, and produces a BAM or CRAM file.

Base Quality Score Recalibration

During Base Quality Score Recalibration, systematic errors in the base quality scores are corrected by applying machine learning to detect and correct for them. This is important for evaluating the correct call of a variant during the variant discovery process. However, this is not needed for all combinations of tools in Sarek. Notably, this should be turned off when having UMI tagged reads or using DragMap (see here) as mapper.

For further reading and documentation see the technical documentation by GATK.

GATK BaseRecalibrator (Spark)

GATK BaseRecalibrator generates a recalibration table based on various co-variates.

To use the corresponding spark implementation GATK BaseRecalibratorSpark, please specify --use_gatk_spark baserecalibrator.

GATK ApplyBQSR (Spark)

GATK ApplyBQSR recalibrates the base qualities of the input reads based on the recalibration table produced by the GATK BaseRecalibrator tool.

Specify --use_gatk_spark baserecalibrator to use GATK ApplyBQSRSpark instead, the respective spark implementation.

The resulting recalibrated CRAM files are delivered to the user. Recalibrated CRAM files are usually 2-3 times larger than the duplicate-marked CRAM files.

CSV files

The CSV files are auto-generated and can be used by Sarek for further processing and/or variant calling.

See the input section in the usage documentation for further reading and documentation on how to make the most of them.

Variant Calling

The results regarding variant calling are collected in {outdir}/variantcalling/. If some results from a variant caller do not appear here, please check out the --tools section in the parameter documentation.

(Recalibrated) CRAM files can used as an input to start the variant calling.

SNVs and small indels

For single nucleotide variants (SNVs) and small indels, multiple tools are available for normal (germline), tumor-only, and tumor-normal (somatic) paired data. For a list of the appropriate tool(s) for the data and sequencing type at hand, please check here.

bcftools

bcftools mpileup generates pileup of a CRAM file, followed by bcftools call and filtered with -i 'count(GT==\"RR\")==0. For further reading and documentation see the bcftools manual.

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. For further documentation take a look here.

FreeBayes

FreeBayes is a Bayesian genetic variant detector designed to find small polymorphisms, specifically SNPs, indels, MNPs, and complex events smaller than the length of a short-read sequencing alignment. For further reading and documentation see the FreeBayes manual.

GATK HaplotypeCaller

GATK HaplotypeCaller calls germline SNPs and indels via local re-assembly of haplotypes.

GATK Germline Single Sample Variant Calling

GATK Single Sample Variant Calling uses HaplotypeCaller in its default single-sample mode to call variants. The VCF that HaplotypeCaller emits errors on the side of sensitivity, therefore they are filtered by first running the CNNScoreVariants tool. This tool annotates each variant with a score indicating the model’s prediction of the quality of each variant. To apply filters based on those scores run the FilterVariantTranches tool with SNP and INDEL sensitivity tranches appropriate for your task.

If the haplotype-called VCF files are not filtered, then Sarek should be run with at least one of the options --dbsnp or --known_indels.

GATK Joint Germline Variant Calling

GATK Joint germline Variant Calling uses Haplotypecaller per sample in gvcf mode. Next, the gVCFs are consolidated from multiple samples into a GenomicsDB datastore. After joint genotyping, VQSR is applied for filtering to produce the final multisample callset with the desired balance of precision and sensitivity.

GATK Mutect2

GATK Mutect2 calls somatic SNVs and indels via local assembly of haplotypes. When --joint_mutect2 is used, Mutect2 subworkflow outputs will be saved in a subfolder named with the patient ID and {patient}.mutect2.vcf.gz file will contain variant calls from all of the normal and tumor samples of the patient. For further reading and documentation see the Mutect2 manual. It is not required, but recommended to have a panel of normals (PON) using at least 40 normal samples to get filtered somatic calls. When using --genome GATK.GRCh38, a panel-of-normals file is available. However, it is highly recommended to create one matching your tumor samples. Creating your own panel-of-normals is currently not natively supported by the pipeline. See here for how to create one manually.

Sentieon DNAscope

Sentieon DNAscope is a variant-caller which aims at outperforming GATK’s Haplotypecaller in terms of both speed and accuracy. DNAscope allows you to use a machine learning model to perform variant calling with higher accuracy by improving the candidate detection and filtering.

The output from Sentieon’s DNAscope can be controlled through the option --sentieon_dnascope_emit_mode for Sarek, see Basic usage of Sentieon functions.

Unless dnascope_filter is listed under --skip_tools in the nextflow command, Sentieon’s DNAModelApply is applied to the unfiltered VCF-files in order to obtain filtered VCF-files.

Sentieon DNAscope joint germline variant calling

In Sentieon’s package DNAscope, joint germline variant calling is done by first running Sentieon’s Dnacope in emit-mode gvcf for each sample and then running Sentieon’s GVCFtyper on the set of gVCF-files. See Basic usage of Sentieon functions for information on how joint germline variant calling can be done in Sarek using Sentieon’s DNAscope.

Sentieon Haplotyper

Sentieon Haplotyper is Sention’s speedup version of GATK’s Haplotypecaller (see above).

The output from Sentieon’s Haplotyper can be controlled through the option --sentieon_haplotyper_emit_mode for Sarek, see Basic usage of Sentieon functions.

Unless haplotyper_filter is listed under --skip_tools in the nextflow command, GATK’s CNNScoreVariants and FilterVariantTranches (see above) is applied to the unfiltered VCF-files in order to obtain filtered VCF-files.

Sentieon Haplotyper joint germline variant calling

In Sentieon’s package DNAseq, joint germline variant calling is done by first running Sentieon’s Haplotyper in emit-mode gvcf for each sample and then running Sentieon’s GVCFtyper on the set of gVCF-files. See Basic usage of Sentieon functions for information on how joint germline variant calling can be done in Sarek using Sentieon’s DNAseq. After joint genotyping, Sentieon’s version of VQSR (VarCal and ApplyVarCal) is applied for filtering to produce the final multisample callset with the desired balance of precision and sensitivity.

Strelka

Strelka is a fast and accurate small variant caller optimized for analysis of germline variation in small cohorts and somatic variation in tumor/normal sample pairs. For further reading and documentation see the Strelka user guide. If Strelka is used for somatic variant calling and Manta is also specified in tools, the output candidate indels from Manta are used according to Strelka Best Practices. For further downstream analysis, take a look here.

Lofreq

Lofreq is a fast and sensitive variant-caller for inferring SNVs and indels from next-generation sequencing data. It makes full use of base-call qualities and other sources of errors inherent in sequencing, which are usually ignored by other methods or only used for filtering. For further reading and documentation see the Lofreq user guide.

Structural Variants

indexcov

indexcov quickly estimate coverage from a whole-genome bam or cram index. A bam index has 16KB resolution and it is used as a coverage estimate . The output is scaled to around 1. So a long stretch with values of 1.5 would be a heterozygous duplication. This is useful as a quick QC to get coverage values across the genome.

Output directory: {outdir}/variantcalling/indexcov/

In addition to the interactive HTML files, indexcov outputs a number of text files:

• <sample>-indexcov.ped: a .ped/.fam file with the inferred sex in the appropriate column if the sex chromosomes were found. the CNX and CNY columns indicating the floating-point estimate of copy-number for those chromosomes. bins.out: how many bins had a coverage value outside of (0.85, 1.15). high values can indicate high-bias samples. bins.lo: number of bins with value < 0.15. high values indicate missing data. bins.hi: number of bins with value > 1.15. bins.in: number of bins with value inside of (0.85, 1.15) p.out: bins.out/bins.in PC1...PC5: PCA projections calculated with depth of autosomes.

• <sample>-indexcov.roc: tab-delimited columns of chrom, scaled coverage cutoff, and $n_samples columns where each indicates the proportion of 16KB blocks at or above that scaled coverage value.

• <sample>-indexcov.bed.gz: a bed file with columns of chrom, start, end, and a column per sample where the values indicate there scaled coverage for that sample in that 16KB chunk.

Manta

Manta calls structural variants (SVs) and indels from mapped paired-end sequencing reads. It is optimized for analysis of germline variation in small sets of individuals and somatic variation in tumor/normal sample pairs. Manta provides a candidate list for small indels that can be fed to Strelka following Strelka Best Practices. For further reading and documentation see the Manta user guide.

TIDDIT

TIDDIT identifies intra and inter-chromosomal translocations, deletions, tandem-duplications and inversions. For further reading and documentation see the TIDDIT manual.

Sample heterogeneity, ploidy and CNVs

ASCAT

ASCAT is a software for performing allele-specific copy number analysis of tumor samples and for estimating tumor ploidy and purity (normal contamination). It infers tumor purity and ploidy and calculates whole-genome allele-specific copy number profiles. The ASCAT process gives several images as output, described in detail in this book chapter. Running ASCAT on NGS data requires that the BAM files are converted into BAF and LogR values. This is done internally using the software AlleleCount. For further reading and documentation see the ASCAT manual.

CNVKit

CNVKit is a toolkit to infer and visualize copy number from high-throughput DNA sequencing data. It is designed for use with hybrid capture, including both whole-exome and custom target panels, and short-read sequencing platforms such as Illumina. For further reading and documentation, see the CNVKit Documentation

Control-FREEC

Control-FREEC is a tool for detection of copy-number changes and allelic imbalances (including loss of heterozygoity (LOH)) using deep-sequencing data. Control-FREEC automatically computes, normalizes, segments copy number and beta allele frequency profiles, then calls copy number alterations and LOH. It also detects subclonal gains and losses and evaluates the most likely average ploidy of the sample. For further reading and documentation see the Control-FREEC Documentation.

Microsatellite instability (MSI)

Microsatellite instability is a genetic condition associated with deficiencies in the mismatch repair (MMR) system which causes a tendency to accumulate a high number of mutations (SNVs and indels). An altered distribution of microsatellite length is associated with a missed replication slippage which would be corrected under normal MMR conditions.

MSIsensorPro

MSIsensorPro is a tool to detect the MSI status of a tumor scanning the length of the microsatellite regions. It requires a normal sample for each tumour to differentiate the somatic and germline cases. For further reading see the MSIsensor paper.

Concatenation

Germline VCFs from DeepVariant, FreeBayes, HaplotypeCaller, Haplotyper, Manta, bcftools mpileup, Strelka, or Tiddit are concatenated with bcftools concat. The field SOURCE is added to the VCF header to report the variant caller.

Variant annotation

This directory contains results from the final annotation steps: two tools are used for annotation, snpEff and VEP. Both results can also be combined by setting --tools merge. All variants present in the called VCF files are annotated. For some variant callers this can mean that the variants are already filtered by PASS, for some this needs to be done during post-processing.

snpEff

snpeff is a genetic variant annotation and effect prediction toolbox. It annotates and predicts the effects of variants on genes (such as amino acid changes) using multiple databases for annotations. The generated VCF header contains the software version and the used command line. For further reading and documentation see the snpEff manual.

VEP

VEP (Variant Effect Predictor), based on Ensembl, is a tool to determine the effects of all sorts of variants, including SNPs, indels, structural variants, CNVs. The generated VCF header contains the software version, also the version numbers for additional databases like Clinvar or dbSNP used in the VEP line. The format of the consequence annotations is also in the VCF header describing the INFO field. For further reading and documentation see the VEP manual.

Currently, it contains:

• Consequence: impact of the variation, if there is any
• Codons: the codon change, i.e. cGt/cAt
• Amino_acids: change in amino acids, i.e. R/H if there is any
• Gene: ENSEMBL gene name
• SYMBOL: gene symbol
• Feature: actual transcript name
• EXON: affected exon
• PolyPhen: prediction based on PolyPhen
• SIFT: prediction by SIFT
• Protein_position: Relative position of amino acid in protein
• BIOTYPE: Biotype of transcript or regulatory feature

plus any additional filed selected via the plugins: dbNSFP, LOFTEE, SpliceAI, SpliceRegion.

BCFtools annotate

BCFtools annotate is used to add annotations to VCF files. The annotations are added to the INFO column of the VCF file. The annotations are added to the VCF header and the VCF header is updated with the new annotations. For further reading and documentation see the BCFtools annotate manual.

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.

The plots display:

• Sequence counts for each sample.
• Sequence Quality Histograms: The mean quality value across each base position in the read.
• Per Sequence Quality Scores: The number of reads with average quality scores. Shows if a subset of reads has poor quality.
• Per Base Sequence Content: The proportion of each base position for which each of the four normal DNA bases has been called.
• Per Sequence GC Content: The average GC content of reads. Normal random library typically have a roughly normal distribution of GC content.
• Per Base N Content: The percentage of base calls at each position for which an N was called.
• Sequence Length Distribution.
• Sequence Duplication Levels: The relative level of duplication found for each sequence.
• Overrepresented sequences: The total amount of overrepresented sequences found in each library.
• Adapter Content: The cumulative percentage count of the proportion of your library which has seen each of the adapter sequences at each position.

FastP

FastP is a tool designed to provide all-in-one preprocessing for FastQ files and is used for trimming and splitting. The tool then determines QC metrics for the processed reads.

Mosdepth

Mosdepth reports information for the evaluation of the quality of the provided alignment data. In short, the basic statistics of the alignment (number of reads, coverage, GC-content, etc.) are summarized and a number of useful graphs are produced. For further reading and documentation see the Mosdepth documentation.

Plots will show:

• cumulative coverage distribution
• absolute coverage distribution
• average coverage per contig/chromosome

NGSCheckMate

NGSCheckMate is a tool for determining whether samples come from the same genetic individual, using a set of commonly heterozygous SNPs. This enables for the detecting of sample mislabelling events. The output includes a text file indicating whether samples have matched or not according to the algorithm, as well as a dendrogram visualising these results.

GATK MarkDuplicates reports

More information in the GATK MarkDuplicates section

Duplicates can arise during sample preparation e.g. library construction using PCR. Duplicate reads can also result from a single amplification cluster, incorrectly detected as multiple clusters by the optical sensor of the sequencing instrument. These duplication artifacts are referred to as optical duplicates. If GATK MarkDuplicates is used, the metrics file generated by the tool is used, if GATK MarkDuplicatesSpark is used the report is generated by GATK4 EstimateLibraryComplexity on the mapped BAM files. For further reading and documentation see the MarkDuplicates manual.

The plot will show:

• duplication statistics

Sentieon Dedup reports

Sentieon’s DNAseq subroutine Dedup produces a metrics report much like the one produced by GATK’s MarkDuplicates. The Dedup metrics are imported into MultiQC as custom content and displayed in a table.

samtools stats

samtools stats collects statistics from CRAM files and outputs in a text format. For further reading and documentation see the samtools manual.

The plots will show:

• Alignment metrics.

bcftools stats

bcftools stats produces a statistics text file which is suitable for machine processing and can be plotted using plot-vcfstats. For further reading and documentation see the bcftools stats manual.

Plots will show:

• Stats by non-reference allele frequency, depth distribution, stats by quality and per-sample counts, singleton stats, etc.
• Note: When using Strelka, there will be no depth distribution plot, as Strelka does not report the INFO/DP field

VCFtools

VCFtools is a program package designed for working with VCF files. For further reading and documentation see the VCFtools manual.

Plots will show:

• the summary counts of each type of transition to transversion ratio for each FILTER category.
• the transition to transversion ratio as a function of alternative allele count (using only bi-allelic SNPs).
• the transition to transversion ratio as a function of SNP quality threshold (using only bi-allelic SNPs).

snpEff reports

snpeff is a genetic variant annotation and effect prediction toolbox. It annotates and predicts the effects of variants on genes (such as amino acid changes) using multiple databases for annotations. For further reading and documentation see the snpEff manual.

The plots will show:

• locations of detected variants in the genome and the number of variants for each location.
• the putative impact of detected variants and the number of variants for each impact.
• the effect of variants at protein level and the number of variants for each effect type.
• the quantity as function of the variant quality score.

VEP reports

VEP (Variant Effect Predictor), based on Ensembl, is a tool to determine the effects of all sorts of variants, including SNPs, indels, structural variants, CNVs. For further reading and documentation see the VEP manual

Reporting

MultiQC

MultiQC is a visualization tool that generates a single HTML report summarizing 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 collect 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.

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.

Reference files

Contains reference folders generated by the pipeline. These files are only published, if --save_reference is set.