nf-core/sarek

CAW is an open-source Nextflow pipeline designed to resolve somatic variants from Whole Genome Sequencing (WGS) data. It is specifically built to process normal/tumor sample pairs (and can accommodate relapse samples) by following GATK best practices for preprocessing and utilizing multiple variant callers for comprehensive analysis.

Key Capabilities:

• Preprocessing: Aligns, realigns, and recalibrates short-read data for both tumor and normal samples in parallel using GATK best practices.
• Somatic Variant Calling: Uses MuTect1, MuTect2, and Strelka to identify somatic SNVs and small indels.
• Structural Variant Calling: Uses Manta to identify structural variants.
• CNV and Heterogeneity Estimation: Uses ASCAT to estimate sample heterogeneity, ploidy, and Copy Number Variations (CNVs).
• Quality Control: Provides integrated quality control reports via MultiQC.
• Flexibility: Can be started from raw FASTQ files, from the realignment step, or directly with specific subsets of variant callers.
• Output: Merges resulting VCF files for downstream processing while retaining individual caller results.

Sarek is an open-source Nextflow pipeline designed for Whole Genome Sequencing (WGS) analysis to resolve both germline and somatic mutations. It follows GATK best practices for short-read data preprocessing and supports parallel processing for tumor/normal pairs.

Key capabilities:

• Sample Modes: Can process normal samples only, tumor/normal pairs, or complex sets including normal, tumor, and multiple relapse samples.
• Variant Calling:
  • Germline: Uses Strelka and GATK HaplotypeCaller.
  • Somatic: Uses MuTect2 and Strelka.
  • Structural Variants: Uses Manta.
• Downstream Analysis: Uses ASCAT to estimate sample heterogeneity, ploidy, and Copy Number Variations (CNVs). VCF files are annotated at the end of the pipeline.
• Quality Control: Provides QC reports via MultiQC.
• Resilience: Includes checkpoints to allow restarting the workflow from different states.
• Reproducibility: Supports execution via Docker and Singularity containers.

The pipeline generates a structured set of output directories within the specified {outdir}. Most plots and summaries are consolidated in the multiqc directory. All paths provided in the documentation are relative to this top-level results directory.

A typical directory structure includes:

• csv/: CSV files generated during processing.
• multiqc/: MultiQC reports and associated data.
• pipeline_info/: Information regarding the pipeline execution.
• preprocessing/: Contains subdirectories for specific preprocessing steps, organized by sample:
  • markduplicates/<sample>/
  • recal_table/<sample>/
  • recalibrated/<sample>/
• reference/: Reference files used during the run.
• reports/: Tool-specific quality control reports, organized by tool name (e.g., reports/<tool1>/).

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

Before calling variants, raw sequencing data must undergo pre-processing to ensure accuracy:

Alignment

Reads are mapped to a specific location in a reference genome. This step is critical for identifying mismatches, insertions, and deletions (INDELs).

Marking Duplicates

To prevent biased allele support, non-independent measurements (duplicates) must be identified. Duplicates can be:

• Library duplicates: Caused by PCR during library preparation.
• Optical duplicates: Caused by the sequencing instrument reading different clusters.

An algorithm identifies these using orientation and 5' position, choosing one representative pair and marking the others to be ignored in downstream steps.

Base Quality Score Recalibration (BQSR)

Sequencing instruments often have systematic errors. BQSR corrects the reported base quality scores by comparing the estimated error rate with the actual (empirical) error rate observed against the reference.

The Mathematical Approach:

1. Calculate empirical error rate ($e_{empirical}$) using a Yates correction: $$e_{empirical} = \frac{n_{mismatches} + 1}{n_{bases} + 2}$$
2. Convert to Phred-scale ($Q_{empirical}$): $$Q_{empirical} = -10 \times \log_{10}(e_{empirical})$$
3. Calculate the difference ($\Delta$) between the average reported quality and the empirical quality.
4. The recalibrated score is the reported score minus the sum of all $\Delta$ values for the bins (e.g., lane, machine cycle, sequencing context) the base belongs to.

Sarek supports two main workflows for processing UMIs: UMI-aware deduplication and consensus read generation.

UMI-aware Deduplication

GATK MarkDuplicates automatically uses UMI-aware deduplication if UMIs are present in the RX tag of the BAM/CRAM file. This occurs if you specify --umi_in_read_header true or use --umi_loc and --umi_len.

Note: GATK MarkDuplicates Spark does not support UMIs.

Consensus Read Generation

1. fgbio: If --umi_read_structure is specified, Sarek uses fgbio to generate consensus reads. You must provide the correct structure (e.g., 2M11S+T for forward-only or 2M11S+T 2M11S+T for both reads) based on your kit.
2. Sentieon: You can enable consensus building within the sentieon_dedup step by setting --sentieon_consensus true. This can be used with or without UMIs.

Configuration Options

• --umi_in_read_header <true|false>: Use if UMIs were transferred to the read header (e.g., via bclconvert).
• --umi_loc <string>: Location of UMI in the read (for fastp extraction).
• --umi_len <integer>: Length of the UMI.
• --umi_read_structure <string>: The fgbio read structure.
• --skip_tools <list>: Use to skip duplicate marking and base quality recalibration if performing post-UMI processing.