SomaticBunny

SomaticBunny is an ensemble pipeline for somatic variant calling that integrates multiple tools to provide comprehensive variant detection in tumor-normal paired samples.

Table of Contents

• Overview
• Prerequisites
• Installation
• Configuration
• Usage
  • Basic Usage
  • Choosing Where Jobs Run
  • Custom Reference Genome
  • Advanced Options
• Output
• Tool Versions
• Benchmark

Overview

The SomaticBunny pipeline integrates multiple variant calling tools to improve detection accuracy for:

• Single Nucleotide Variants (SNVs)
• Insertions and Deletions (INDELs)
• Structural Variants (SVs)
• Copy Number Variants (CNVs)

This pipeline is designed to work with paired tumor-normal whole-genome or whole-exome sequencing data.

Prerequisites

• Conda/Mamba for environment management
• A SLURM cluster is recommended for WGS data (see Choosing Where Jobs Run)

Installation

1. Clone this repository

2. Install Nextflow as a conda environment

3. Make sure you have the following required files in your working directory:

   • main.nf
   • nextflow.config
   • conf/base.config
   • sample.csv

4. Prepare your sample.csv file according to your analysis requirements:

   For starting from fastq (click to expand)

   patient,sample,status,fastq_1,fastq_2
   UCSD101,A,normal,/PATH/TO/ERR5285401_1.fastq.gz,/PATH/TO/ERR5285401_2.fastq.gz
   UCSD101,A,tumor,/PATH/TO/ERR5285402_1.fastq.gz,/PATH/TO/ERR5285402_2.fastq.gz
   UCSD101,B,tumor,/PATH/TO/ERR5285403_1.fastq.gz,/PATH/TO/ERR5285403_2.fastq.gz

   For starting from bams (click to expand)

   patient,sample,status,bam,bai,sex
   UCSD101,A,normal,/PATH/TO/UCSD101_A_normal.bam,/PATH/TO/UCSD101_A_normal.bam.bai,XY
   UCSD101,A,tumor,/PATH/TO/UCSD101_A_tumor.bam,/PATH/TO/UCSD101_A_tumor.bam.bai,XY
   UCSD101,B,tumor,/PATH/TO/UCSD101_B_tumor.bam,/PATH/TO/UCSD101_B_tumor.bam.bai,XY

   For ASCAT analysis (click to expand)

   patient,sample,status,fastq_1,fastq_2,sex
   UCSD101,A,normal,/PATH/TO/ERR5285401_1.fastq.gz,/PATH/TO/ERR5285401_2.fastq.gz,XY
   UCSD101,A,tumor,/PATH/TO/ERR5285402_1.fastq.gz,/PATH/TO/ERR5285402_2.fastq.gz,XY
   UCSD101,B,tumor,/PATH/TO/ERR5285403_1.fastq.gz,/PATH/TO/ERR5285403_2.fastq.gz,XY

Configuration

For memory management, modify the following parameters if needed:

• Set $params.mkdup_temp_dir in main.nf to your ideal location (Default: $projectDir/mkdup_tmp)
• Set $workDir in nextflow.config to your ideal location (Default: ./work)

// In main.nf
params.mkdup_temp_dir = "/path/to/ideal/mkdup_tmp/folder/"

// In nextflow.config
workDir = "/path/to/ideal/intermediate/work/folder/"

Usage

Basic Usage

1. Request an interactive compute node (TSCC example; on another cluster, use your own partition and account):

srun -N 1 -n 1 -c 16 --mem 250G -t 80:00:00 -p platinum -q hcp-ddp302 -A ddp302 --pty bash

2. Activate your Nextflow environment:

conda activate env_nf

3. Set TSCC temporary directory and Nextflow tmp directory, and memory:

export TMPDIR=/path/to/restricted/folder/

export NXF_OPTS="-Djava.io.tmpdir=${TMPDIR} -Xms4g -Xmx16g"

4. Run the pipeline, choosing where jobs run with -profile (see Choosing Where Jobs Run). This example runs on TSCC:

nextflow run main.nf -profile tscc_slurm,conda --type [genome|exome] --genome [GRCh38|GRCh37|mm39|RN7] --first_step [mapping|markdup|recalibration|variant_calling] --tool [ascat,manta,...] --database_path [/path/to/SomaticBunny_database] --ref [/path/to/ref/] --ref_fai [/path/to/ref_fai] --ref_dict [/path/to/ref_dict] --bed [/path/to/bed] --bed_tbi [/path/to/bed_tbi]

Pipeline Steps and Tools

The SomaticBunny pipeline workflow is divided into multiple steps that can be run individually or in sequence using the --first_step parameter.

Available first_step

First_step	Description	Required Input
mapping	This will start from performing alignment of FASTQ files to the reference genome	fastq_1 and fastq_2
markdup	This will start from mark duplicate reads using the aligned BAM files	raw bam and bai
recalibration	This will start from performing base quality score recalibration	mark duplicated bam and bai
variant_calling	This will start from executing variant calling tools	recalibrated bam and bai

Available reference genomes

Genomes
GRCh38
GRCh37
mm39
RN7

Available Tools

Use the --tool parameter to specify which variant callers to run. You can select multiple tools by separating them with commas.

Category	Available Tools	Description
SV Callers	manta	Structural variant and indel caller
	delly	Integrated structural variant detection
CNV Callers	ascat	Allele-specific copy number analysis (requires sex information)
	cnvkit	Copy number variation detection from targeted DNA sequencing

Example of running multiple tools:

# Run Manta for SVs and ASCAT for CNVs
nextflow run main.nf -profile tscc_slurm,conda --genome GRCh38 --type genome --first_step mapping --tool manta,ascat

Choosing Where Jobs Run

Use -profile (one dash, it is a Nextflow option) to choose the scheduler, combined with conda for the software. Profiles are separated by a comma, with no space.

Profile	Runs on
tscc_slurm,conda	UCSD TSCC: each job is submitted to SLURM with the lab's partition, QoS and account
slurm,conda	Any SLURM cluster: each job is submitted to SLURM with the cluster's defaults, or with your own options (see below)
conda	The current machine only (no scheduler)

Warning

Without tscc_slurm or slurm, all jobs run on the machine where Nextflow is started, which is not practical for WGS data.

Running on other SLURM clusters

If your cluster needs a partition, account or QoS, put them in a small config file of your own and pass it with -c. You don't need to edit SomaticBunny's nextflow.config, so your settings are kept when you update the pipeline.

1. Create my_cluster.config (any name works):

// my_cluster.config
process.clusterOptions = '-p <partition> -A <account>'

2. Run with the slurm profile and your config:

nextflow run main.nf -profile slurm,conda -c my_cluster.config --type genome ...

process.clusterOptions is passed to every job as extra sbatch options, so anything sbatch accepts works, e.g. '--partition=compute --account=mylab --qos=normal'.

If your nodes are smaller than the default requests (up to 16 CPUs and 250 GB for some steps), cap what any job can request:

nextflow run main.nf -profile slurm,conda -c my_cluster.config --max_cpus 32 --max_memory 192.GB --max_time 72.h ...

Software

SomaticBunny installs its tools with Conda (-profile conda, or -profile mamba for faster environment creation). Container engines such as Docker and Singularity are not supported.

Custom Reference Genome

When starting the pipeline from BAM files (--first_step markdup, recalibration, or variant_calling), you must provide your own reference genome files. This is because BAM files are already aligned to a specific reference, and the pipeline needs the matching reference for downstream analysis.

Required parameters

Parameter	Description	Example
--ref	Reference FASTA file	/path/to/hg38_genome.fa
--ref_fai	FASTA index file	/path/to/hg38_genome.fa.fai
--ref_dict	Sequence dictionary	/path/to/hg38_genome.dict
--bed	Callable regions BED file (bgzipped)	/path/to/callable_regions.bed.gz
--bed_tbi	Tabix index for the BED file	/path/to/callable_regions.bed.gz.tbi

All five files must be provided together. The pipeline will exit with an error if any are missing.

Example

nextflow run main.nf \
  -profile tscc_slurm,conda \
  --type genome \
  --genome GRCh38 \
  --first_step variant_calling \
  --sample sample.csv \
  --database_path /path/to/SomaticBunny_database \
  --ref /path/to/hg38_genome.fa \
  --ref_fai /path/to/hg38_genome.fa.fai \
  --ref_dict /path/to/hg38_genome.dict \
  --bed /path/to/hg38_genome_callable_regions.bed.gz \
  --bed_tbi /path/to/hg38_genome_callable_regions.bed.gz.tbi

When a custom reference is provided, the pipeline will automatically generate the scattered interval lists required by GATK tools (Mutect2, BaseRecalibrator, etc.) from your reference. When starting from mapping, these parameters are optional and will fall back to the pipeline's built-in defaults for the selected --genome.

Important

Your custom reference genome must use the chr prefix for contig names (e.g., chr1, chr2, ..., chrX, chrY). The pipeline's built-in database files (dbSNP, gnomAD, Panel of Normals, known indels, etc.) use chr-prefixed contig names. A reference without the chr prefix (e.g., 1, 2, ..., X, Y) will cause tools such as Mutect2, MuSE, Strelka, SAGE, and BQSR to fail due to contig name mismatches.

Different GRCh38 builds that include additional or fewer contigs (e.g., alt contigs, decoys, HLA) are supported as long as the standard chromosomes use the chr prefix.

Advanced Options

Publishing Intermediate Files

By default, intermediate files are not saved to reduce disk usage. Use the --publish parameter to save specific file types:

nextflow run main.nf -profile tscc_slurm,conda --type genome --genome GRCh38 --first_step mapping --publish raw_bam,markdup_bam

Available options:

• raw_bam: Publish raw BAM files from BWA-MEM
• markdup_bam: Publish mark duplicates BAM files
• recal_bam: Publish recalibrated BAM files

Resuming Failed Runs

If your pipeline terminates with an error or the interactive node is killed, resume with:

nextflow run main.nf -profile tscc_slurm,conda --type genome --genome GRCh38 --first_step mapping --tool ascat,manta -resume

Email Notifications

Receive completion notification:

nextflow run main.nf -profile tscc_slurm,conda --type genome --genome GRCh38 --first_step mapping --tool ascat,manta -N your_email@example.com

Output

All results and reports are stored in the RESULTS folder with the following structure:

RESULTS/
├── ASCAT (optional)
├── BAM (optional)
├── CNVkit (optional)
├── Conpair
├── Delly (optional)
├── FASTQC
├── MANTA (optional)
├── MKDUP
├── mosdepth
├── MuSE
├── Mutect2
├── RECALIBRATE (optional)
├── REPORT
├── SAGE
└── STRELKA

Tool Versions

Category	Tool	Version
QC	FastQC	v0.12.1
Data Processing	Picard	v3.4.0
	samtools	v1.24
	bwa-mem2	v2.3
	Conpair	v0.2
	GATK4	v4.6.2.0
	mosdepth	v0.3.14
SNV/INDEL Callers	Strelka2	v2.9.10
	Mutect2	v4.6.2.0
	SAGE	v4.2
	MuSE2	v2.1.2
Structural Variants	Delly	v2.6.0
	Manta	v1.6.0
Copy Number Variants	CNVkit	v0.9.13
	ASCAT	v3.2.0

Benchmark

Code to reproduce the benchmark against nf-core/sarek is in benchmark/.