Sarek: A portable workflow for whole-genome sequencing analysis of germline and somatic variants

Operation: workflow overview and software

Sarek offers a portable workflow for germline and somatic variant detection, annotation and quality control based on WGS, WES or gene panel data, using a range of state-of-the-art software and data resources in the field (Table 1, Figure 1). In the pre-processing step, sequence reads are aligned to the reference genome with BWA-MEM (Li, 2013), followed by deduplication and recalibration with GATK (McKenna et al., 2010). For germline samples, single-nucleotide variants and small insertion/deletions are detected with HaplotypeCaller (McKenna et al., 2010) and Strelka2 (Kim et al., 2018), and structural variations are detected with Manta (Chen et al., 2016) and TIDDIT (Eisfeldt et al., 2017). For somatic samples, somatic single-base mutations (SSM) and small somatic insertion/deletion mutations (SIM) are detected by GATK4 Mutect2 (Cibulskis et al., 2013) and Strelka2 (Kim et al., 2018). Somatic structural variants (including copy-number variation), as well as ploidy and sample purity are detected by Manta (Chen et al., 2016), ASCAT (Van Loo et al., 2010), and Control-FREEC (Boeva et al., 2012). All variants are annotated for potential functional effects with snpEff (Cingolani et al., 2012) and VEP (McLaren et al., 2016). Importantly, Sarek also generates a wide range of quality control metrics using FastQC, QualiMap (Okonechnikov et al., 2016), BCFtools (Li, 2011), Samtools (Li et al., 2009), and VCFtools (Danecek et al., 2011), visualized as an aggregated quality control review across samples with MultiQC (Ewels et al., 2016). All software currently included in Sarek are selected based on the criteria that they should be of high quality, well-maintained, and with robust installation and running performances. Additional alternative or complementing software will be added to Sarek in later updates, based on the input and engagement of the user community.

Figure 1. Schematic overview of the Sarek workflow for analysis of germline and somatic variants.

A schematic overview including some of the main analysis software implemented in the Sarek workflow. A more comprehensive list of the currently implemented software is given in Table 1.

Portability and reproducibility

Sarek is implemented in Nextflow (Di Tommaso et al., 2017), a workflow language designed specifically for bioinformatics applications. Nextflow has a transparent design, making the Sarek code easy to read, adjust and extend. Sarek has well-functioning error reporting to diagnose e.g. software or hardware errors during a run, and incomplete runs are easily restarted from any stage in the workflow process. Compared to the Bpipe workflow language (used in for example DNAp), Nextflow offers superior support for different execution environments, like Slurm, Sun Grid Engine, LSF and Kubernetes, and includes native support for cloud compute environments including Google Cloud and AWS. Support for AWS batch gives the possibility to easily distribute thousands of batch jobs on Amazon Web Services. Sarek is part of a rapidly growing community effort of well documented and community-tested Nextflow pipelines, and adheres to the nf-core portability and documentation guidelines (Ewels et al., 2019). To facilitate easy installation and to ensure reproducibility, all Sarek required tools are installed in Conda, and then pushed to DockerHub (https://hub.docker.com/), making Sarek and all its dependencies directly accessible from a Conda environment, or as Docker or Singularity (Kurtzer et al., 2017) containers. While Docker is a widely appreciated container solution, it is not always allowed at high-performance computing centers because of the involved security risks, making Singularity the preferred choice at these sites (Kurtzer et al., 2017). This is of particular importance for computer environments designed for handling of sensitive personal data, where a high level of data security has to be maintained across multiple projects and users.

Implementation: equipment and resource usage

Sarek can be installed and executed on any POSIX-compatible computer system. To run a full WGS analysis, including both germline and somatic variants from a tumour/normal dataset with 90x/90x read coverage, we recommend a minimum of 16 cores on a node with 128 GB RAM, and at least 4 TB available free storage (in addition to the initial FASTQ files) in the input/output working directory. Of this, about 1.4 TB will be allocated for BAM files, annotated VCF files and CNV files, but excluding GVCF files (Table 2). At the end of the run, 2.3 TB temporary data can be removed, unless the user plans to perform re-runs from intermediate processing states. Many processes are distributed across cores by dividing the genome into smaller chunks, each being handled as a separate core job, with all the results being merged and sorted in a final step. Some of the used software are parallelized by design, while for others Sarek uses a scatter-gather approach to efficiently distribute the processing load across CPU cores and reduce the wall clock runtime.

Installation and testing

Sarek is run from a computer system with a local installation of Nextflow and support for either Conda environments, Docker or Singularity containers. Nextflow can automatically fetch the Sarek source code from GitHub. All software dependencies are encapsulated in Docker or Singularity containers which are downloaded from Docker Hub, or built in a new Conda environment using Bioconda (Grüning et al., 2018). As such, cumbersome software installations by the user are completely avoided. Configuration files allow tailoring to specific user needs. Sarek comes with a small test dataset and a suite of tests to verify the installation. This is also used for Continuous Integration testing with GitHub Actions.

Input data

For a somatic variant analysis, the user should provide the sequencing FASTQ files from both tumour and normal control tissue from the same individual, described in a tab-delimited TSV file (here: samples.tsv). Each line of the TSV file contains information about a sequence data file, including: The identifier of the individual, the gender (XX or XY), the status of the sample (0 for Normal or 1 for Tumour), the identifier of the sample, the sequencing lane (if samples are multiplexed across multiple lanes), and the paths to the FASTQ file of the first and second read in the read-pair. Relapse samples from the same individual are also supported.

Running sarek on WGS data with singularity containers

Running Sarek with Singularity container on a computer system supporting Java 8 requires only installation of Nextflow and Singularity. A full analysis run starting from FASTQ files including mapping, recalibration, variant calling and annotation, as well as generating a full QC report can be invoked by a single Nextflow command:

> nextflow run nf-core/sarek -r 2.5.2 -profile singularity --input samples.tsv --tools Mutect2,Strelka,Manta,TIDDIT,ASCAT,ControlFREEC,snpEff,VEP

Nextflow will recognize the workflow name and will download the specified version (2.5.2) of the pipeline from GitHub, including the corresponding container, as well as fetching the required reference files from AWS-iGenomes. The default reference genome is human GRCh38, but Sarek also supports GRCh37 and nearly 30 other genomes directly accessible from iGenomes. Alternatively, users can manually supply Sarek with other reference genomes. Non-default parameters and links to local reference files are handled in accordance with nf-core guidelines. User configuration profiles can be stored locally or centrally at https://github.com/nf-core/configs.

Output

A full Sarek run will produce a large number of output files, but the main results consist of (i) a set of annotated variants in VCF files from the various included tools for both germline and somatic variants, (ii) tumour sample purity and ploidy results for somatic samples, and (iii) a broad set of QC metrics. A detailed description of all output files is given at the Sarek documentation pages. While Sarek will report variants from all callers included in the run, it is up to the user to decide how to combine and filter the results from different callers, since the optimal post-processing will depend on the particular samples and research questions at hand.

Source data

European Genome-phenome Archive: A comprehensive assessment of somatic mutation detection in cancer using whole genome sequencing. https://www.ebi.ac.uk/ega/datasets/EGAD00001001859. Read sets EGAR00001387019-24 and EGAR00001387025-32 were analysed. These data are held under restricted access. Readers wishing to apply for access to the data must first apply through the ICGC Data Access Compliance Office (https://icgc.org/daco) and complete the data access form. Access will be granted to those whose projects conform to the goals and policies of ICGC. Help with completing the data access form is available at https://icgc.org/daco/help-guide-section.

Sequence Read Archive: NIST Genome in a Bottle, ^~300X sequencing of HG001 (NA12878). ftp://ftp-trace.ncbi.nlm.nih.gov/giab/ftp/data/NA12878/NIST_NA12878_HG001_HiSeq_300x/, read set folders 131219_D00360_005_BH814YADXX [SRA accession number SRR2052337 - SRR2052339, SRR2052342, SRR2052345, SRR20523428], and 131219_D00360_006_AH81VLADXX [SRA accession number SRX1049774 -SRX1049779]). These data are publicly available for direct download.

The workflow itself comes with a prebuilt profile with a complete configuration for automated testing, including links to a small test dataset.