diff --git a/.gitattributes b/.gitattributes index 8f61a8e..59c24d9 100644 --- a/.gitattributes +++ b/.gitattributes @@ -1,2 +1,19 @@ # SCM syntax highlighting pixi.lock linguist-language=YAML linguist-generated=true + +# Force windows not to reformat line-endings +*.sh text eol=lf +*.bash text eol=lf +*.py text eol=lf +*.pl text eol=lf +*.R text eol=lf +*.pbs text eol=lf +*.slurm text eol=lf +*.config text eol=lf +*.yaml text eol=lf +*.nf text eol=lf +*.md text eol=lf +*.Rmd text eol=lf +*.groovy text eol=lf +*.json text eol=lf +*.ssi text eol=lf diff --git a/.github/workflows/nextflow-test-conda.yml b/.github/workflows/nextflow-test-conda.yml index 81bb8ec..81164fd 100644 --- a/.github/workflows/nextflow-test-conda.yml +++ b/.github/workflows/nextflow-test-conda.yml @@ -11,6 +11,9 @@ jobs: test: runs-on: ubuntu-latest + env: + NXF_SYNTAX_PARSER: v1 + steps: - name: Checkout repository uses: actions/checkout@v4 diff --git a/.github/workflows/nextflow-test-container.yml b/.github/workflows/nextflow-test-container.yml index 45b9eed..8d332dc 100644 --- a/.github/workflows/nextflow-test-container.yml +++ b/.github/workflows/nextflow-test-container.yml @@ -10,7 +10,10 @@ on: jobs: test: runs-on: ubuntu-latest - + + env: + NXF_SYNTAX_PARSER: v1 + steps: - name: Checkout repository uses: actions/checkout@v4 @@ -28,7 +31,7 @@ jobs: key: docker-cache-${{ runner.os }}-${{ hashFiles('**/nextflow.config', '**/modules.config') }} restore-keys: | docker-cache-${{ runner.os }}- - + - name: Install Nextflow run: | curl -s https://get.nextflow.io | bash diff --git a/.gitignore b/.gitignore index 1c7e94f..3f436c8 100644 --- a/.gitignore +++ b/.gitignore @@ -14,7 +14,6 @@ dist/ downloads/ eggs/ .eggs/ -lib/ lib64/ parts/ sdist/ @@ -185,3 +184,4 @@ vcftable.tsv pixi.toml pixi.lock .DS_Store +.vscode diff --git a/README.md b/README.md index aa772aa..bf0955f 100644 --- a/README.md +++ b/README.md @@ -43,19 +43,24 @@ Configurations can be supplied a configuration file, see ```config.json```, and nextflow run main.nf \ (-with-docker/-with-apptainer/-with-conda) \ --pooltable \ - --decodetable \ + (--decodetable \) --reference_genome \ --bedfile \ --ploidy ``` -The ```pooltable.tsv``` should connect (user assigned) pool id's to input FASTQ files; one entry for each pool. +The ```pooltable.tsv``` should connect (user assigned) pool id's and their row/column arrangement to input FASTQ files; one tab-separated line for each pool. ```Bash -pool_row_1 path/to/sample1_R1.fq.gz path/to/sample1_R2.fq.gz -pool_column_1 path/to/sample2_R1.fq.gz path/to/sample2_R2.fq.gz +pool_1 row path/to/sample1_R1.fq.gz path/to/sample1_R2.fq.gz +pool_2 row path/to/sample2_R1.fq.gz path/to/sample2_R2.fq.gz +pool_3 column path/to/sample3_R1.fq.gz path/to/sample3_R2.fq.gz +pool_4 column path/to/sample4_R1.fq.gz path/to/sample4_R2.fq.gz ``` -The ```decodetable.tsv``` should map (user assigned) individual id's in the matrix to the corresponding row and column id's of each pool; one entry for each element in the matrix. +The optional ```decodetable.tsv``` should map (user assigned) individual id's in the matrix to the corresponding row and column id's of each pool; one entry for each element in the matrix. ```Bash -individual1 pool_row_1 pool_column_1 +individual1 pool_1 pool_3 +individual2 pool_2 pool_3 +individual3 pool_1 pool_4 +individual4 pool_2 pool_4 ``` ### 6. Pipeline output The workflow will output a results folder containing multiple config dependent output files: @@ -63,9 +68,11 @@ The workflow will output a results folder containing multiple config dependent o results ├── pinpointables.vcf # Merged VCF file containing all assigned variants ├── cram/ # CRAM files for each pool +├── context/ # TSV and JSON files with matrix information and pool-individual linkage. ├── logs/ # Log files for each process ├── variants/ # VCF files for each pool ├── variant_tables/ # TSV files converted from pool VCFs +├── variant_compilation/ # TSV files with aggregated variants, annotations and rescue probabilities └── pinpoint_variants/ ├── all_pins/ # All pinpointables for each sample in individual vcfs (*note) ├── unique_pins/ # All unique pinpointables for each sample in individual vcfs (*note) @@ -75,6 +82,34 @@ results ``` A central files is the ```pinpointables.vcf```. This file contains all individually assigned variants. Since each variant contains information from two pools, these a presented as the sample columns: ROW and COLUMN. +## Variant annotation +Two annotations workflows are currently available. SnpEff for VCF output files and VEP for tabular compiled output. They can be applied by adding the following configurations and absolute paths to the JSON params file. +```Bash +config.json +{ + # Annotate output VCFs with SnpEff and clinvar + annotate: true + snpeff_db: "GRCh38.99" + snpeff_config: "snpEff.config" + snpeff_cache: "cache/" + clinvar_db: "clinvar_20230903.vcf.gz" # (*.tbi in same folder) + + # Annotate tabular output with VEP + annotate_vep: true + vep_cache: "cache/" # (Current version 111.x) + # Optional VEP input: + danmac_db: "danmac.vcf.gz" # (*.tbi in same folder) + blacklist_bed: "hg38-blacklist.v2.sorted.bed.gz" # (*.tbi in same folder) + repeatmasker_bed: "repeatmasker.sorted.bed.gz" # (*.tbi in same folder) + gnomad_vcf: "gnomad.vcf.bgz" # (*.tbi in same folder) + utr_file: "uORF_5UTR_GRCh38_PUBLIC.txt" + alphamissense_tsv: "AlphaMissense_hg38.tsv.gz" + loftee_gerp_bw: "gerp_conservation_scores.homo_sapiens.GRCh38.bw" + loftee_human_ancestor: "human_ancestor.fa.gz" + loftee_sqlite: "loftee.sql" +} +``` + # Workflow repository contents: ```Bash @@ -95,6 +130,7 @@ DoBSeqWF │ │ ├── pooltable.tsv │ │ ├── snvlist.tsv │ │ └── target_calling.bed +│ ├── filter/ # Filter model modules and parameters │ └── helper_scripts │ └── simulator.py # Script for simulating minimal pipeline data ├── bin # Executable pipeline scripts @@ -106,11 +142,12 @@ DoBSeqWF ├── envs │ └── / │ └── environment.yaml # Conda environment definitions +├── lib/ # Pipeline groovy utility functions ├── main.nf # Main workflow ├── modules/ │ └── .nf # Module scripts ├── subworkflows/ -│ └── .nf # Module scripts +│ └── .nf # Subworkflow scripts ├── next.pbs # Helper script for running on NGC-HPC └── nextflow.config # Workflow parameters ``` @@ -180,8 +217,8 @@ predisposed ├── / │ ├── DoBSeqWF/ # Clone repository here │ │ ├── config.json # Configuration file - │ │ ├── pooltable.tsv # Pool table (create with helper script) - │ │ └── decodetable.tsv # Decode table (we need a convention for this) + │ │ ├── pooltable.tsv # Pool table + │ │ └── decodetable.tsv # Decode table │ └── results │ ├── cram/ # CRAM files for each pool │ ├── logs/ # Log files for each process diff --git a/assets/data/test_data/pooltable.tsv b/assets/data/test_data/pooltable.tsv index 77ba74d..1b733c3 100644 --- a/assets/data/test_data/pooltable.tsv +++ b/assets/data/test_data/pooltable.tsv @@ -1,4 +1,4 @@ -B0_H0 assets/data/test_data/pools/B0_H0_1.fq.gz assets/data/test_data/pools/B0_H0_2.fq.gz -B0_H1 assets/data/test_data/pools/B0_H1_1.fq.gz assets/data/test_data/pools/B0_H1_2.fq.gz -B0_V0 assets/data/test_data/pools/B0_V0_1.fq.gz assets/data/test_data/pools/B0_V0_2.fq.gz -B0_V1 assets/data/test_data/pools/B0_V1_1.fq.gz assets/data/test_data/pools/B0_V1_2.fq.gz +B0_H0 row assets/data/test_data/pools/B0_H0_1.fq.gz assets/data/test_data/pools/B0_H0_2.fq.gz +B0_H1 row assets/data/test_data/pools/B0_H1_1.fq.gz assets/data/test_data/pools/B0_H1_2.fq.gz +B0_V0 column assets/data/test_data/pools/B0_V0_1.fq.gz assets/data/test_data/pools/B0_V0_2.fq.gz +B0_V1 column assets/data/test_data/pools/B0_V1_1.fq.gz assets/data/test_data/pools/B0_V1_2.fq.gz diff --git a/assets/data/test_data/target_calling.bed b/assets/data/test_data/target_calling.bed index d1ac8d8..64865d4 100644 --- a/assets/data/test_data/target_calling.bed +++ b/assets/data/test_data/target_calling.bed @@ -1,2 +1,2 @@ -small_ref 1 710 -small_ref 720 1300 +small_ref 1 710 gene1 +small_ref 720 1300 gene2 diff --git a/bin/build_matrix.py b/bin/build_matrix.py new file mode 100755 index 0000000..77c8e8f --- /dev/null +++ b/bin/build_matrix.py @@ -0,0 +1,219 @@ +#!/usr/bin/env python3 +""" +Build a matrix context from pools (ground truth) and an optional decode table. + +Inputs +------ +1) Pools TSV (required): columns: pool_iddimension + - dimension must be "row" or "column" + - pool_id must be unique + +2) Decode TSV (optional): columns: sample_idrow_pool_idcol_pool_id + - each (row_pool_id, col_pool_id) should appear at most once + +Outputs +------- +A) JSON context (for machines): + - schema_version, layout_hash + - pools (rows/columns with index + labels) + - matrix metadata and full cell list (row/col pool, indices, labels, alias, sample_id) + +B) TSV (for humans / spreadsheets): + - one row per intersection with: row/col pool IDs, indices, labels, cell_label (A01), + alias (same as cell_label by default), sample_id (if any) + +Usage +----- +python build_matrix.py \ + --pools pools.tsv \ + --decode decode.tsv \ + --out-json matrix_context.json \ + --out-tsv matrix_cells.tsv \ + --pad-width auto +""" +import argparse, csv, hashlib, json, sys, re +from pathlib import Path +from typing import Dict, List, Tuple, Optional + +def read_pools(pools_path: Path) -> Tuple[List[str], List[str]]: + rows, cols = [], [] + with pools_path.open() as fh: + reader = csv.reader(fh, delimiter="\t") + first = next(reader, None) + if first is None: + raise ValueError("Pools TSV is empty.") + records = [first] + for r in reader: + records.append(r) + for rec in records: + if len(rec) < 2: + raise ValueError(f"Bad pools row (need 2 columns): {rec}") + pool_id, dim = rec[0].strip(), rec[1].strip().lower() + if dim not in {"row", "column"}: + raise ValueError(f"Invalid dimension '{dim}' for pool '{pool_id}'") + (rows if dim == "row" else cols).append(pool_id) + # Check uniqueness + all_ids = rows + cols + if len(all_ids) != len(set(all_ids)): + dup = [p for p in set(all_ids) if all_ids.count(p) > 1] + raise ValueError(f"Duplicate pool_id(s) in pools TSV: {dup}") + if not rows or not cols: + raise ValueError("Need at least one row pool and one column pool.") + return rows, cols + +def read_decode(decode_path: Optional[Path]) -> Dict[Tuple[str,str], str]: + mapping = {} + if not decode_path: return mapping + with decode_path.open() as fh: + reader = csv.reader(fh, delimiter="\t") + + first = next(reader, None) + if first is None: + return mapping + + records = [first] + for r in reader: + records.append(r) + for rec in records: + if len(rec) < 3: + raise ValueError(f"Bad decode row (need 3 columns): {rec}") + sample_id, rpid, cpid = rec[0].strip(), rec[1].strip(), rec[2].strip() + key = (rpid, cpid) + if key in mapping and mapping[key] != sample_id: + raise ValueError(f"Decode conflict for {key}: '{mapping[key]}' vs '{sample_id}'") + mapping[key] = sample_id + return mapping + +def column_label_from_index(idx: int) -> str: + """Create excel style column label 1->A, 26->Z, 27->AA""" + label = [] + n = idx + while n > 0: + n -= 1 + label.append(chr(65 + (n % 26))) + n //= 26 + return "".join(reversed(label)) + +def layout_hash(rows_ordered: List[str], cols_ordered: List[str]) -> str: + m = hashlib.sha256() + m.update(b"rows\0"); [m.update((r+"\0").encode()) for r in rows_ordered] + m.update(b"cols\0"); [m.update((c+"\0").encode()) for c in cols_ordered] + return m.hexdigest()[:16] + +def build_context(rows: List[str], + cols: List[str], + decode: Dict[Tuple[str,str], str], + pad_width: Optional[int]) -> dict: + # Force deterministic ordering + rows_ord = sorted(rows, key=str.lower) + cols_ord = sorted(cols, key=str.lower) + + # Indices and labels + row_index = {pid: i for i, pid in enumerate(rows_ord)} + col_index = {pid: i for i, pid in enumerate(cols_ord)} + col_label = {pid: column_label_from_index(col_index[pid] +1 ) for pid in cols_ord} + n_rows = len(rows_ord) + width = pad_width if pad_width and pad_width > 0 else len(str(n_rows)) + row_label = {pid: str(row_index[pid] + 1).zfill(width) for pid in rows_ord} + + # Pools section + pools_rows = [{"pool_id": pid, "index": row_index[pid], "label": row_label[pid]} for pid in rows_ord] + pools_cols = [{"pool_id": pid, "index": col_index[pid], "label": col_label[pid]} for pid in cols_ord] + + # Cells (Cartesian product, ie. all combinations) + cells = [] + for rpid in rows_ord: + for cpid in cols_ord: + ri, ci = row_index[rpid], col_index[cpid] + rl, cl = row_label[rpid], col_label[cpid] + cell = { + "row_pool_id": rpid, + "col_pool_id": cpid, + "row_index": ri, + "col_index": ci, + "row_label": rl, + "col_label": cl, + "cell_label": f"{cl}{rl}", + "alias": f"{cl}{rl}", + "sample_id": decode.get((rpid, cpid), None), + } + cells.append(cell) + + ctx = { + "schema_version": "1.0", + "layout_hash": layout_hash(rows_ord, cols_ord), + "pools": { + "rows": pools_rows, + "columns": pools_cols, + }, + "matrix": { + "n_rows": len(rows_ord), + "n_cols": len(cols_ord), + "row_label_width": width, + "cells": cells, + }, + } + return ctx + +def write_json(ctx: dict, path: Path): + with path.open("w") as fh: + json.dump(ctx, fh, indent=2, sort_keys=False) + +def write_tsv(ctx: dict, path: Path): + fields = [ + "row_pool_id","col_pool_id", + "row_index","col_index", + "row_label","col_label", + "cell_label","alias","sample_id" + ] + with path.open("w", newline="") as fh: + w = csv.writer(fh, delimiter="\t") + w.writerow(fields) + for cell in ctx["matrix"]["cells"]: + w.writerow([ + cell["row_pool_id"], cell["col_pool_id"], + cell["row_index"], cell["col_index"], + cell["row_label"], cell["col_label"], + cell["cell_label"], cell["alias"], cell.get("sample_id") + ]) + +def write_decode(ctx: dict, path: Path): + with path.open("w", newline="") as fh: + w = csv.writer(fh, delimiter="\t") + for cell in ctx["matrix"]["cells"]: + w.writerow([ + cell["alias"], cell["row_pool_id"], cell["col_pool_id"] + ]) + +def main(argv=None): + p = argparse.ArgumentParser(description="Build matrix context from pools and optional decode.") + p.add_argument("--pools", required=True, type=Path, help="Pools TSV (pool_id, dimension)") + p.add_argument("--decode", required=False, type=Path, help="Decode TSV (sample_id, row_pool_id, col_pool_id)") + p.add_argument("--out-json", required=True, type=Path, help="Output JSON path") + p.add_argument("--out-tsv", required=True, type=Path, help="Output TSV path") + p.add_argument("--out-decode", required=True, type=Path, help="Output decode TSV path") + p.add_argument("--pad-width", default="auto", help='Column label zero-padding width (int or "auto")') + args = p.parse_args(argv) + + rows, cols = read_pools(args.pools) + decode = read_decode(args.decode) if args.decode else {} + pad_width = None if str(args.pad_width).lower() == "auto" else int(args.pad_width) + + # Validate decode keys are known pools + for (rpid, cpid) in decode.keys(): + if rpid not in rows: + raise ValueError(f"Decode references unknown row_pool_id '{rpid}'") + if cpid not in cols: + raise ValueError(f"Decode references unknown col_pool_id '{cpid}'") + + ctx = build_context(rows, cols, decode, pad_width) + args.out_json.parent.mkdir(parents=True, exist_ok=True) + args.out_tsv.parent.mkdir(parents=True, exist_ok=True) + args.out_decode.parent.mkdir(parents=True, exist_ok=True) + write_json(ctx, args.out_json) + write_tsv(ctx, args.out_tsv) + write_decode(ctx, args.out_decode) + +if __name__ == "__main__": + main() + \ No newline at end of file diff --git a/bin/pin_basic.py b/bin/pin_basic.py new file mode 100755 index 0000000..2538f39 --- /dev/null +++ b/bin/pin_basic.py @@ -0,0 +1,425 @@ +#!/usr/bin/env python +import json, sys, argparse +from pathlib import Path +from dataclasses import dataclass, fields +from pysam import VariantFile +from typing import Dict, Set, List, Tuple, Optional, Any, Union + +## This is pin_basic.py - script for processing pipeline output for easy downstream analysis. +# mads - 2025-09-05 + +@dataclass +class VariantEntry: + varid: str + CHROM: str + POS: int + REF: str + ALT: str + QUAL: float = 0.0 + AC: int = 0 + FS: float = 0.0 + SOR: float = 0.0 + GQ: int = 0 + DP: int = 0 + REF_AD: int = 0 + ALT_AD: int = 0 + +def get_other_variants(pools, idx): + return set().union(*(pools[i] for i in range(len(pools)) if i != idx)) + +def pin(vcf_folder: Path, caller: str, matrix_context: str) -> Dict[str, Dict[str, Set[str]]]: + """ + Performs the actual pinpointing logic and returns a dictionary of variants per sample. + + Parameters: + vcf_folder (Path): Path to the folder containing pool vcfs. + caller (str): Variant caller name. + + Returns: + A dictionary mapping sample IDs to their corresponding unique variants and all pinpointable variants. + """ + sample_variants: Dict[str, Dict[str, Set[str]]] = {} + + row_vcfs = [vcf_folder / f"{pool}.{caller}.vcf.gz" for pool in matrix_context.row_pools] + column_vcfs = [vcf_folder / f"{pool}.{caller}.vcf.gz" for pool in matrix_context.col_pools] + + row_pools = [set(collect_variant_ids(vcf)) for vcf in row_vcfs] + column_pools = [set(collect_variant_ids(vcf)) for vcf in column_vcfs] + + # Pinning logic: + for cell in matrix_context._ctx["matrix"]["cells"]: + sample_id = cell["alias"] + row_idx = cell["row_index"] + col_idx = cell["col_index"] + # Extract all variants from all other pools in each dimension + other_row_variants = get_other_variants(row_pools, row_idx) + other_col_variants = get_other_variants(column_pools, col_idx) + + # Extract unique and pool-specific variants + unique_row_variants = row_pools[row_idx].difference(other_row_variants) + unique_col_variants = column_pools[col_idx].difference(other_col_variants) + unique_pins = unique_row_variants.intersection(unique_col_variants) + + # Extract all pinnable variants. Only required to be unique in one dimension + unique_one_dimension_row = unique_row_variants.intersection(column_pools[col_idx]) + unique_one_dimension_col = unique_col_variants.intersection(row_pools[row_idx]) + all_pins = unique_one_dimension_row.union(unique_one_dimension_col) + + sample_variants[sample_id] = { + "unique_pins": unique_pins, + "all_pins": all_pins + } + + return sample_variants + + +def collect_variant_ids(vcf_path: str, sample_id: Optional[str] = None) -> Set[str]: + """ + Read VCF and return set of variant IDs for a single sample in format 'chrom:pos:ref:alt'. + """ + with VariantFile(vcf_path) as vcf: + if sample_id is not None: + if sample_id not in vcf.header.samples: + raise ValueError(f"Sample '{sample_id}' not found in VCF header.") + chosen_sample = sample_id + else: + if len(vcf.header.samples) == 0: + raise ValueError("No samples found in VCF header.") + elif len(vcf.header.samples) > 1: + raise ValueError( + f"Multi-sample VCF ({len(vcf.header.samples)} samples) " + f"requires specifying sample_id. Available: {list(vcf.header.samples)}" + ) + chosen_sample = next(iter(vcf.header.samples)) + + variant_ids: Set[str] = set() + + for rec in vcf.fetch(): + if not rec.alts or len(rec.alts) != 1: + raise ValueError( + f"Non-biallelic record at {rec.chrom}:{rec.pos} " + f"(REF={rec.ref}, ALTS={rec.alts}). Expected exactly one ALT." + ) + + # Only include if sample has a non-missing call + sample = rec.samples[chosen_sample] + if sample.alleles and None not in sample.alleles: + variant_id = f"{rec.chrom}:{rec.pos}:{rec.ref}:{rec.alts[0]}" + variant_ids.add(variant_id) + + return variant_ids + +def _first_scalar(x, default=0): + if x is None: + return default + if isinstance(x, (list, tuple)): + return x[0] if x else default + return x + +def collect_variant_entries(vcf_path: str, sample_id: Optional[str] = None) -> List[VariantEntry]: + """ + Read VCF and return VariantEntry per record. + + If sample_id is given, extract that sample's fields; otherwise use the first sample if present. + For site-only VCFs, sample-level fields remain 0. + """ + vcf = VariantFile(vcf_path) + + # Choose sample once (or None if site-only) + if sample_id is not None: + if sample_id not in vcf.header.samples: + raise ValueError(f"Requested sample_id '{sample_id}' not found in VCF header.") + chosen_sample = sample_id + else: + chosen_sample = next(iter(vcf.header.samples)) if len(vcf.header.samples) > 0 else None + + out: List[VariantEntry] = [] + + for rec in vcf.fetch(): + # Only accept one alt allele - multi-allelic sites is split before input. + if not rec.alts or len(rec.alts) != 1: + raise ValueError( + f"Non-biallelic record at {rec.chrom}:{rec.pos} " + f"(REF={rec.ref}, ALTS={rec.alts}). Expected exactly one ALT." + ) + alt = rec.alts[0] + + # Site-level fields + entry = VariantEntry( + varid=f"{rec.chrom}:{rec.pos}:{rec.ref}:{alt}", + CHROM=rec.chrom, + POS=rec.pos, + REF=rec.ref, + ALT=alt, + QUAL=rec.qual or 0.0, + AC=int(_first_scalar(rec.info.get("AC") if "AC" in vcf.header.info else 0, 0)), + FS=float(_first_scalar(rec.info.get("FS") if "FS" in vcf.header.info else 0.0, 0.0)), + SOR=float(_first_scalar(rec.info.get("SOR") if "SOR" in vcf.header.info else 0.0, 0.0)), + ) + + # Sample-level fields (if any) + if chosen_sample is not None: + s = rec.samples[chosen_sample] + entry.GQ = int(s.get("GQ", 0) or 0) + entry.DP = int(s.get("DP", 0) or 0) + ad = s.get("AD", []) + + entry.REF_AD = int(ad[0]) if isinstance(ad, (list, tuple)) and len(ad) > 0 else 0 + entry.ALT_AD = int(ad[1]) if isinstance(ad, (list, tuple)) and len(ad) > 1 else 0 + out.append(entry) + return out + +def get_variant_entry(vcf_path: str, variant_id: str, sample_id: Optional[str] = None) -> VariantEntry: + """ + Locate a specific variant, extract and return a VariantEntry. + """ + try: + chrom, pos_str, ref, alt = variant_id.split(":") + pos = int(pos_str) + except ValueError: + raise ValueError(f"Invalid variant_id format '{variant_id}'. Expected 'chrom:pos:ref:alt'") + + with VariantFile(vcf_path) as vcf: + if sample_id is not None: + if sample_id not in vcf.header.samples: + raise ValueError(f"Requested sample_id '{sample_id}' not found in VCF header.") + chosen_sample = sample_id + else: + chosen_sample = next(iter(vcf.header.samples)) if len(vcf.header.samples) > 0 else None + + # Search for the specific variant + for rec in vcf.fetch(chrom, pos - 1, pos): + if (rec.pos == pos and + rec.ref == ref and + rec.alts == (alt,)): + + # Only accept one alt allele - multi-allelic sites should be split before input + if not rec.alts or len(rec.alts) != 1: + raise ValueError( + f"Non-biallelic record at {rec.chrom}:{rec.pos} " + f"(REF={rec.ref}, ALTS={rec.alts}). Expected exactly one ALT." + ) + + # Site-level fields + entry = VariantEntry( + varid=variant_id, + CHROM=rec.chrom, + POS=rec.pos, + REF=rec.ref, + ALT=rec.alts[0], + QUAL=rec.qual or 0.0, + AC=int(_first_scalar(rec.info.get("AC") if "AC" in vcf.header.info else 0, 0)), + FS=float(_first_scalar(rec.info.get("FS") if "FS" in vcf.header.info else 0.0, 0.0)), + SOR=float(_first_scalar(rec.info.get("SOR") if "SOR" in vcf.header.info else 0.0, 0.0)), + ) + + # Sample-level fields (if any) + if chosen_sample is not None: + s = rec.samples[chosen_sample] + entry.GQ = int(s.get("GQ", 0) or 0) + entry.DP = int(s.get("DP", 0) or 0) + ad = s.get("AD", []) + + entry.REF_AD = int(ad[0]) if isinstance(ad, (list, tuple)) and len(ad) > 0 else 0 + entry.ALT_AD = int(ad[1]) if isinstance(ad, (list, tuple)) and len(ad) > 1 else 0 + return entry + raise ValueError(f"Variant '{variant_id}' not found in VCF '{vcf_path}'") + +class MatrixContext: + """ + A helper class for accessing matrix context data generated by build_matrix.py. + + Provides convenience functions to navigate DoBSeq matrices safely. + """ + + def __init__(self, ctx: dict): + self._ctx = ctx + + # Sort pools by index for consistent ordering + self._rows = sorted(ctx["pools"]["rows"], key=lambda r: r["index"]) + self._cols = sorted(ctx["pools"]["columns"], key=lambda c: c["index"]) + + # Pool ID lists (ordered) + self.row_pools: List[str] = [r["pool_id"] for r in self._rows] + self.col_pools: List[str] = [c["pool_id"] for c in self._cols] + + # Index lookups + self.row_index_by_pool: Dict[str, int] = {r["pool_id"]: r["index"] for r in self._rows} + self.col_index_by_pool: Dict[str, int] = {c["pool_id"]: c["index"] for c in self._cols} + + # Label lookups + self.row_label_by_pool: Dict[str, str] = {r["pool_id"]: r["label"] for r in self._rows} + self.col_label_by_pool: Dict[str, str] = {c["pool_id"]: c["label"] for c in self._cols} + + # Combined dimension lookup + self._dim_by_pool: Dict[str, str] = { + **{p: "row" for p in self.row_pools}, + **{p: "column" for p in self.col_pools} + } + + # Label lists (ordered) + self.row_labels: List[str] = [r["label"] for r in self._rows] + self.col_labels: List[str] = [c["label"] for c in self._cols] + + # Cell lookups - single source of truth + cells = ctx["matrix"]["cells"] + self.cells: Dict[Tuple[str, str], Dict[str, Any]] = { + (c["row_pool_id"], c["col_pool_id"]): c for c in cells + } + + # Alternative access methods + self.by_alias: Dict[str, Dict[str, Any]] = { + c["alias"]: c for c in cells + } + self.by_sample: Dict[str, Dict[str, Any]] = { + c["sample_id"]: c for c in cells if c.get("sample_id") + } + + @classmethod + def load(cls, path: Union[Path, str]) -> "MatrixContext": + """Load matrix context from JSON file.""" + with open(path, "r") as fh: + data = json.load(fh) + return cls(data) + + @property + def shape(self) -> Tuple[int, int]: + """Return matrix shape as (n_rows, n_cols).""" + return (self._ctx["matrix"]["n_rows"], self._ctx["matrix"]["n_cols"]) + + @property + def n_rows(self) -> int: + """Number of rows in the matrix.""" + return self._ctx["matrix"]["n_rows"] + + @property + def n_cols(self) -> int: + """Number of columns in the matrix.""" + return self._ctx["matrix"]["n_cols"] + + @property + def layout_hash(self) -> str: + """Get the layout hash for this matrix configuration.""" + return self._ctx["layout_hash"] + + def pool_dim(self, pool_id: str) -> str: + """Get dimension ('row' or 'column') for a pool ID.""" + if pool_id not in self._dim_by_pool: + raise ValueError(f"Unknown pool_id: {pool_id}") + return self._dim_by_pool[pool_id] + + def pool_index(self, pool_id: str) -> int: + """Get the 1-based index for a pool ID.""" + if pool_id in self.row_index_by_pool: + return self.row_index_by_pool[pool_id] + elif pool_id in self.col_index_by_pool: + return self.col_index_by_pool[pool_id] + else: + raise ValueError(f"Unknown pool_id: {pool_id}") + + def pool_label(self, pool_id: str) -> str: + """Get the formatted label for a pool ID.""" + if pool_id in self.row_label_by_pool: + return self.row_label_by_pool[pool_id] + elif pool_id in self.col_label_by_pool: + return self.col_label_by_pool[pool_id] + else: + raise ValueError(f"Unknown pool_id: {pool_id}") + + def cell(self, row_pool: str, col_pool: str) -> Optional[Dict[str, Any]]: + """Get cell data for the intersection of row and column pools.""" + return self.cells.get((row_pool, col_pool)) + + def cell_from_alias(self, alias: str) -> Optional[Dict[str, Any]]: + """Get cell data from cell alias (e.g., 'A01').""" + return self.by_alias.get(alias) + + def cell_from_sample(self, sample_id: str) -> Optional[Dict[str, Any]]: + """Get cell data from sample ID.""" + return self.by_sample.get(sample_id) + + def cell_label(self, row_pool_id: str, col_pool_id: str) -> Optional[str]: + """Get the cell label (e.g., 'A01') for pool intersection.""" + cell = self.cell(row_pool_id, col_pool_id) + return cell["cell_label"] if cell else None + + def get_samples(self) -> List[str]: + """Get all sample IDs in the matrix.""" + return list(self.by_sample.keys()) + + def __repr__(self) -> str: + return (f"MatrixContext({self.n_rows}×{self.n_cols}, " + f"{len(self.get_samples())} samples)") + +def main(): + p = argparse.ArgumentParser(description="Generate pool and pinpoint variant tables from matrix context + pool VCFs") + p.add_argument("--context", "-c", type=Path, required=True, help="Path to matrix_context.json") + p.add_argument("--vcf-folder", "-v", type=Path, required=True, help="Folder containing per-pool VCFs") + p.add_argument("--output", "-o", type=Path, required=True, help="Output folder for result TSVs") + p.add_argument("--caller", "-C", type=str, default="GATK", help="Variant caller name (used in file naming, default: GATK)") + p.add_argument("--pin-type", "-p", type=str, choices=["unique_pins", "all_pins"], default="unique_pins", help="Which pinpoint set to use in pinpoint_variants.tsv (default: unique_pins)") + args = p.parse_args() + + output_path = args.output + matrix_context = args.context + pool_vcfs_path = args.vcf_folder + pin_type = args.pin_type + caller = args.caller + + output_path.mkdir(parents=True, exist_ok=True) + mc = MatrixContext.load(matrix_context) + + pools = (mc.row_pools,mc.col_pools) + + ###### + ## Create table with all variants: + + variant_entry_fields = [f.name for f in fields(VariantEntry)] + header = ["dim", "pool_id", "pool_index", "pool_label", *variant_entry_fields] + + with open(output_path / 'all_pool_variants.tsv', 'w') as fout: + print(*header, sep="\t", file=fout) + for n, dim in enumerate(['row','column']): + for pool_id in pools[n]: + filepath = pool_vcfs_path / (f"{pool_id}.{caller}.vcf.gz") + variants = collect_variant_entries(filepath) + for v in variants: + values = [getattr(v, col) for col in variant_entry_fields] + print(dim, pool_id, mc.pool_index(pool_id), mc.pool_label(pool_id), *values, sep="\t", file=fout) + + ###### + # Create table with pinpointables only: + pins = pin(vcf_folder=pool_vcfs_path, matrix_context=mc, caller=caller) + + common_fields = ['varid','CHROM','POS','REF','ALT'] + variant_entry_header_fields = [f"{dim}.{f.name}" for dim in ["row","column"] for f in fields(VariantEntry) if f.name not in common_fields] + header = ["uvarid", "sample_alias", "sample_id", "row_id", "row_index", "row_label", "column_id", "column_index", "column_label", *common_fields, *variant_entry_header_fields] + + with open(output_path / 'pinpoint_variants.tsv', 'w') as fout: + print(*header, sep="\t", file=fout) + for cell in mc._ctx["matrix"]["cells"]: + sample_alias = cell["alias"] + sample_id = cell["sample_id"] + row_label = cell["row_label"] + col_label = cell["col_label"] + row_idx = cell["row_index"] + col_idx = cell["col_index"] + row_id = cell["row_pool_id"] + col_id = cell["col_pool_id"] + for pinpoint in pins[sample_alias][pin_type]: + row_info = get_variant_entry(vcf_path=pool_vcfs_path / (f"{row_id}.{caller}.vcf.gz"), variant_id=pinpoint) + col_info = get_variant_entry(vcf_path=pool_vcfs_path / (f"{col_id}.{caller}.vcf.gz"), variant_id=pinpoint) + + assert row_info.varid == col_info.varid, "Row and column variant IDs do not match" + + uvarid = f"{sample_alias}:{row_info.varid}" + print( + uvarid, sample_alias, sample_id, row_id, row_idx, row_label, col_id, col_idx, col_label, + *[getattr(row_info, col) for col in variant_entry_fields], + *[getattr(col_info, col) for col in variant_entry_fields if col not in common_fields], + sep='\t', + file=fout + ) + +if __name__ == '__main__': + main() diff --git a/bin/report.Rmd b/bin/report.Rmd new file mode 100644 index 0000000..4bf6141 --- /dev/null +++ b/bin/report.Rmd @@ -0,0 +1,269 @@ +--- +title: "DoBSeq Report" +output: html_document +date: "2025-09-17" +params: + matrix_context: null + pinpoints: null + all_variants: null + annotations: null + rescue: null + output: null +--- + +```{r setup, include=FALSE} +knitr::opts_chunk$set(echo = F, message = FALSE, warning = FALSE) +knitr::opts_knit$set(root.dir = params$project_root) +``` + +```{r} +library('ggthemes') +library('jsonlite') +library("dplyr") +library("stringr") +library("ggplot2") +library("ggbeeswarm") +library("tidyr") +library("readr") +library("tibble") +library("forcats") +library("scales") +``` + + +```{r} +have <- function(p) !is.null(p) && length(p) == 1 && fs::file_exists(p) +``` + + +```{r eval=have(params$pinpoints)} +matrix_context_path <- params$matrix_context +pinpoints_path <- params$pinpoints +all_variants_path <- params$all_variants +annotations_path <- params$annotations +rescue_path <- params$rescue +output_path <- params$output +``` + +```{r eval=!have(params$pinpoints)} +here::i_am("bin/report.Rmd") + +matrix_context_path <- here::here('results','context','matrix_context.json') +pinpoints_path <- here::here('results','variant_compilation','pinpoint_variants.tsv') +all_variants_path <- here::here('results','variant_compilation','all_pool_variants.tsv') +annotations_path <- here::here('results','vep_annotate','annotations_w_varid.tsv') +rescue_path <- here::here('results','rescue_probabilities','predictions.tsv') +output_path <- here::here('reports') +``` + +```{r} +pinvars <- + read_tsv(pinpoints_path, + col_types = cols(.default = col_character())) %>% + type_convert(guess_integer = TRUE) %>% + mutate(row_label = as.integer(row_label), + is_pool_pin = 1) + +poolvars <- + read_tsv(all_variants_path, + col_types = cols(.default = col_character())) %>% + type_convert(guess_integer = TRUE) +``` + +```{r eval=have(params$rescue)} +rescue <- + read_tsv(rescue_path) %>% + mutate(row_label = as.integer(row_label)) %>% + select(-sample_id) + +pinvars <- + pinvars %>% + full_join(rescue, by=c('uvarid','row_id','column_id','row_label','column_label','varid','sample_alias', 'is_pool_pin')) %>% + mutate(is_pool_pin = if_else(is.na(is_pool_pin),0,is_pool_pin)) +``` + +```{r eval=have(params$annotations)} + +annotations <- + read_tsv(annotations_path, + na = '-', col_types = cols(.default = col_character())) %>% + type_convert(guess_integer = TRUE) %>% + mutate(is_lof = if_else(LoF == 'HC',1,0)) %>% + mutate(clinvar_stars = case_when( + str_detect(ClinVar_CLNREVSTAT, "practice_guideline") ~ 4, + str_detect(ClinVar_CLNREVSTAT, "reviewed_by_expert_panel") ~ 3, + str_detect(ClinVar_CLNREVSTAT, "criteria_provided,_multiple_submitters,_no_conflicts") ~ 2, + str_detect(ClinVar_CLNREVSTAT, "criteria_provided,_conflicting_interpretations") ~ 1, + str_detect(ClinVar_CLNREVSTAT, "criteria_provided,_single_submitter") ~ 1, + TRUE ~ NA_real_ + )) %>% + mutate(is_p = if_else(str_detect(ClinVar_CLNSIG, 'Pathogenic') | str_detect(ClinVar_CLNSIG, "Likely_pathogenic"), 1,0)) %>% + mutate(is_lofp = if_else(is_lof == 1 | is_p == 1,1,0)) + +pinvars <- + pinvars %>% + left_join(annotations, by='varid') + +poolvars <- + poolvars %>% + left_join(annotations, by='varid') +``` + +```{r} +matrix_ctx <- fromJSON(matrix_context_path) +matrix_size <- matrix_ctx$matrix$n_rows + +all_sample_ids <- + bind_rows(matrix_ctx$matrix$cells) %>% + filter(!is.na(sample_id) & sample_id != "") %>% + pull(sample_id) %>% + unique() + +all_sample_alias <- + bind_rows(matrix_ctx$matrix$cells) %>% + filter(!is.na(alias) & alias != "") %>% + pull(alias) %>% + unique() + +all_pool_alias <- + bind_rows(matrix_ctx$matrix$cells) %>% + filter(!is.na(sample_id) & sample_id != "") %>% + pull(sample_id) %>% + unique() + +# Create dynamic factor levels +row_levels <- 1:matrix_size +col_levels <- LETTERS[1:matrix_size] + +pinvars <- + pinvars %>% + mutate(row_factor = fct_rev(factor(row_label, levels = row_levels, ordered = TRUE)), + column_factor = factor(column_label, levels = col_levels, ordered = TRUE)) +``` + +#### Assigned variants per individual +```{r} +plot_data <- + pinvars %>% + filter(is_pool_pin == 1) %>% + group_by(sample_id,row_factor, column_factor) %>% + summarise(private_variants = n()) + +plot_data %>% + ggplot(aes(x = column_factor, y = row_factor, fill = private_variants)) + + geom_point(shape = 22, size = 14, color='black') + + scale_x_discrete(drop = FALSE, position = "top") + + scale_y_discrete(drop = FALSE, limits = levels(plot_data$row_factor)) + + geom_text(aes(label=private_variants),size=5, color='white') + + ggthemes::scale_fill_continuous_tableau() + + coord_fixed() + + labs(x = "COLUMN POOL", y = "ROW POOL") + + theme_minimal() + + theme(legend.title = element_blank(), + legend.position = '', + panel.grid.minor = element_line(size = 2), + panel.grid.major = element_line(size = 2), + axis.title = element_text(face = 'bold', color = '#7f8c8d', size = 20), + axis.text = element_text(face = 'bold', color = '#7f8c8d', size = 15)) + + ggtitle('All assigned') +``` + +```{r eval=have(params$annotations)} +plot_data <- + pinvars %>% + filter(is_pool_pin == 1) %>% + filter(is_lofp == 1) %>% + group_by(sample_id,row_factor, column_factor) %>% + summarise(private_variants = n()) + +plot_data %>% + ggplot(aes(x = column_factor, y = row_factor, fill = private_variants)) + + geom_point(shape = 22, size = 14, color='black') + + scale_x_discrete(drop = FALSE, position = "top") + + scale_y_discrete(drop = FALSE, limits = levels(plot_data$row_factor)) + + geom_text(aes(label=private_variants),size=5, color='white') + + ggthemes::scale_fill_continuous_tableau() + + coord_fixed() + + labs(x = "COLUMN POOL", y = "ROW POOL") + + theme_minimal() + + theme(legend.title = element_blank(), + legend.position = '', + panel.grid.minor = element_line(size = 2), + panel.grid.major = element_line(size = 2), + axis.title = element_text(face = 'bold', color = '#7f8c8d', size = 20), + axis.text = element_text(face = 'bold', color = '#7f8c8d', size = 15)) + + ggtitle('pLoF/p only') +``` + + + + + +```{r eval=have(params$annotations) & have(params$rescue)} +plot_data <- + pinvars %>% + filter((probability >= 0.5 & is_single) | is_pool_pin == 1) %>% + filter(is_lofp == 1) %>% + group_by(sample_id,row_factor, column_factor) %>% + summarise(private_variants = n()) + +plot_data %>% + ggplot(aes(x = column_factor, y = row_factor, fill = private_variants)) + + geom_point(shape = 22, size = 14, color='black') + + scale_x_discrete(drop = FALSE, position = "top") + + scale_y_discrete(drop = FALSE, limits = levels(plot_data$row_factor)) + + geom_text(aes(label=private_variants),size=5, color='white') + + ggthemes::scale_fill_continuous_tableau() + + coord_fixed() + + labs(x = "COLUMN POOL", y = "ROW POOL") + + theme_minimal() + + theme(legend.title = element_blank(), + legend.position = '', + panel.grid.minor = element_line(size = 2), + panel.grid.major = element_line(size = 2), + axis.title = element_text(face = 'bold', color = '#7f8c8d', size = 20), + axis.text = element_text(face = 'bold', color = '#7f8c8d', size = 15)) + + ggtitle('pLoF/p only after rescue') +``` + +#### Variant count by pool +```{r} +poolvars %>% + group_by(pool_id, dim) %>% + summarise(n = n()) %>% + ungroup() %>% + mutate(din = if_else(dim == 'column', 1,0)) %>% + ggplot(aes(x=reorder(pool_id,din),y=n,fill=dim)) + + geom_col() + + coord_flip() + + theme_minimal() + + xlab('Pool ID') + + ylab('Variant count') + + ggthemes::scale_fill_tableau() + + theme_bw() + + theme(axis.text=element_text(size=10), + legend.position = 'bottom', + legend.title = element_blank()) +``` + +#### Mean VAF count by individual +```{r} +plot_df <- + pinvars %>% + group_by(sample_alias) %>% + mutate(mean_vaf = (column.ALT_AD+row.ALT_AD)/(column.DP+row.DP)) + +ggplot(plot_df, aes(x = sample_alias, y = mean_vaf)) + + ggbeeswarm::geom_beeswarm(corral='wrap',corral.width = 0.9,alpha=0.6) + + geom_hline(yintercept = 1/(matrix_size*2), linetype = 'dashed', color='darkgrey') + + geom_hline(yintercept = 1/matrix_size, linetype = 'dashed', color='darkgrey') + + coord_flip() + + theme_minimal() + + ggthemes::scale_color_fivethirtyeight() + + theme(legend.title = element_blank()) + + ylab('Mean VAF') + + xlab('Sample ID') +``` + + + diff --git a/conf/container.config b/conf/container.config index 9713ad4..9b162af 100644 --- a/conf/container.config +++ b/conf/container.config @@ -2,29 +2,74 @@ // mads 2024-02-09 // Base container image for shell scripts +// Docker runoptions fixes ownership confusion in some containers process { - container = 'ubuntu:20.04' + container = 'ubuntu:24.04' + containerOptions = '-u $(id -u):$(id -g)' } params { - 'container' { - gatk = "quay.io/biocontainers/gatk4:4.5.0.0--py36hdfd78af_0" - samtools = "quay.io/biocontainers/samtools:1.21--h50ea8bc_0" - fastqc = "quay.io/biocontainers/fastqc:0.12.1--hdfd78af_0" - bwa = "quay.io/biocontainers/mulled-v2-e5d375990341c5aef3c9aff74f96f66f65375ef6:2d15960ccea84e249a150b7f5d4db3a42fc2d6c3-0" - bbmap = "quay.io/biocontainers/bbmap:38.89-h1296035_0" - fgbio = "quay.io/biocontainers/fgbio:1.5.1-0" - multiqc = "quay.io/biocontainers/multiqc:1.26--pyhdfd78af_0" - mosdepth = "quay.io/biocontainers/mosdepth:0.3.10--h4e814b3_1" - lofreq = "quay.io/biocontainers/lofreq:2.1.5--py37h3a01142_1" - freebayes = "quay.io/biocontainers/freebayes:1.2.0--py35h82df9c4_2" - bcftools = "quay.io/biocontainers/bcftools:1.19--h8b25389_1" - pinpy = "madscort/dswf_pinpy:1.0" - snpeff = "quay.io/biocontainers/snpeff:4.3--2" - sambamba = "quay.io/biocontainers/sambamba:1.0.1--h6f6fda4_1" - octopus = "quay.io/biocontainers/octopus:0.7.4--hacbf28e_1" - deepvariant = "quay.io/biocontainers/deepvariant:1.5.0--py36hf3e76ba_0" - r_env = "rocker/tidyverse:latest" - filter_variants = "madscort/dswf_filter:1.0" + container { + docker { + gatk = "quay.io/biocontainers/gatk4:4.5.0.0--py36hdfd78af_0" + samtools = "quay.io/biocontainers/samtools:1.21--h50ea8bc_0" + fastqc = "quay.io/biocontainers/fastqc:0.12.1--hdfd78af_0" + bwa = "quay.io/biocontainers/mulled-v2-e5d375990341c5aef3c9aff74f96f66f65375ef6:2d15960ccea84e249a150b7f5d4db3a42fc2d6c3-0" + fgbio = "quay.io/biocontainers/fgbio:1.5.1-0" + multiqc = "quay.io/biocontainers/multiqc:1.26--pyhdfd78af_0" + mosdepth = "quay.io/biocontainers/mosdepth:0.3.10--h4e814b3_1" + lofreq = "quay.io/biocontainers/lofreq:2.1.5--py37h3a01142_1" + freebayes = "quay.io/biocontainers/freebayes:1.3.6--hbfe0e7f_2" + bcftools = "quay.io/biocontainers/bcftools:1.19--h8b25389_1" + snpeff = "quay.io/biocontainers/snpeff:4.3--2" + snpsift = "quay.io/biocontainers/snpsift:5.1d--hdfd78af_0" + sambamba = "quay.io/biocontainers/sambamba:1.0.1--h6f6fda4_1" + deepvariant = "quay.io/biocontainers/deepvariant:1.5.0--py36hf3e76ba_0" + + // Custom and non-singularity containers + ngsep = "madscort/ngsep:5.1.0" + octopus = "quay.io/biocontainers/octopus:0.7.4--hacbf28e_1" + bbmap = "community.wave.seqera.io/library/bbmap_pigz:07416fe99b090fa9" + picard = "community.wave.seqera.io/library/picard:3.4.0--e9963040df0a9bf6" + r_env = "madscort/dswf_r:1.0" + pinpy = "madscort/dswf_pinpy:1.0" + filter_variants = "madscort/dswf_filter:1.0" + vep = "madscort/vep:111.0" + marbl = "madscort/marbl:0.6.0" + } + singularity { + gatk = "https://depot.galaxyproject.org-singularity-gatk4-4.5.0.0--py36hdfd78af_0" + bwa = "https://depot.galaxyproject.org-singularity-mulled-v2-e5d375990341c5aef3c9aff74f96f66f65375ef6-2d15960ccea84e249a150b7f5d4db3a42fc2d6c3-0" + samtools = "https://depot.galaxyproject.org-singularity-samtools-1.20--h50ea8bc_0" + picard = "https://depot.galaxyproject.org-singularity-picard-3.1.1--hdfd78af_0" + fastqc = "https://depot.galaxyproject.org-singularity-fastqc-0.12.1--hdfd78af_0" + fastp = "https://depot.galaxyproject.org-singularity-fastp-0.23.4--h5f740d0_0" + mosdepth = "https://depot.galaxyproject.org-singularity-mosdepth-0.3.8--hd299d5a_0" + multiqc = "https://depot.galaxyproject.org-singularity-multiqc-1.23--pyhdfd78af_0" + deepvariant = "https://depot.galaxyproject.org-nf-core-deepvariant-1.5.0" + strelka = "https://depot.galaxyproject.org-singularity-strelka-2.9.10--h9ee0642_1" + freebayes = "https://depot.galaxyproject.org-singularity-freebayes-1.3.6--hbfe0e7f_2" + lofreq = "https://depot.galaxyproject.org-singularity-lofreq-2.1.5--py38h588ecb2_4" + manta = "https://depot.galaxyproject.org-singularity-manta-1.6.0--h9ee0642_1" + bcftools = "https://depot.galaxyproject.org-singularity-bcftools-1.20--h8b25389_0" + sambamba = "https://depot.galaxyproject.org/singularity/sambamba:1.0.1--h6f6fda4_0" + snpeff = "https://depot.galaxyproject.org-singularity-snpeff-5.1--hdfd78af_2" + snpsift = "https://depot.galaxyproject.org/singularity/snpsift:5.1d--hdfd78af_0" + bedtools = "https://depot.galaxyproject.org-singularity-bedtools-2.31.1--hf5e1c6e_0" + sentieon = "https://depot.galaxyproject.org-singularity-sentieon-202308.02--h43eeafb_0" + fgbio = "https://depot.galaxyproject.org-singularity-fgbio-2.1.0--hdfd78af_0" + bbmap = "https://depot.galaxyproject.org/singularity/bbmap:38.22--h14c3975_1" + + // Custom and non-singularity containers + ngsep = "madscort/ngsep:5.1.0" + octopus = "quay.io/biocontainers/octopus:0.7.4--hacbf28e_1" + bbmap = "community.wave.seqera.io/library/bbmap_pigz:07416fe99b090fa9" + picard = "community.wave.seqera.io/library/picard:3.4.0--e9963040df0a9bf6" + r_env = "madscort/dswf_r:1.0" + pinpy = "madscort/dswf_pinpy:1.0" + filter_variants = "madscort/dswf_filter:1.0" + vep = "madscort/vep:111.0" + marbl = "madscort/marbl:0.6.0" + } } -} \ No newline at end of file +} diff --git a/conf/ngc.config b/conf/ngc.config index 694be07..f55def4 100644 --- a/conf/ngc.config +++ b/conf/ngc.config @@ -5,41 +5,49 @@ profiles { outputDir = "results" project = "icope_staging_r" cacheDir = "/ngc/projects2/dp_00005/scratch/nextflow" - fgbio = "java -Xmx4g -XX:+AggressiveHeap -jar /services/tools/fgbio/1.5.1/fgbio.jar" - snpeff = "java -Xmx4g -jar /services/tools/snpeff/5.0e/snpEff.jar" - snpsift = "java -Xmx4g -jar /services/tools/snpeff/5.0e/SnpSift.jar" + reference_genome = "/ngc/projects2/dp_00005/data/predisposed/resources/reference/Homo_sapiens_assembly38_noAlt.fna" + snpeff_db = 'GRCh38.99' snpeff_config = "/ngc/projects2/dp_00005/data/dwf/databases/snpeff/snpEff.config" snpeff_cache = "/ngc/projects2/dp_00005/data/dwf/databases/snpeff/data" - clinvar_db = "/ngc/projects2/dp_00005/data/dwf/databases/clinvar_20230903.vcf.gz" + vep_cache = "/ngc/shared/nf_tools/references/vep_cache" + clinvar_db = "/ngc/projects2/dp_00005/data/dwf/databases/clinvar_20230903.vcf.gz" + danmac_db = "/ngc/projects2/dp_00005/data/dwf/databases/danmac/danmac.vcf.gz" + blacklist_bed = "/ngc/projects2/dp_00005/data/dwf/databases/encode_blacklist/hg38-blacklist.v2.sorted.bed.gz" + repeatmasker_bed = "/ngc/projects2/dp_00005/data/dwf/databases/repeat_masker/repeatmasker.sorted.bed.gz" + gnomad_vcf = "/ngc/projects2/dp_00005/data/dwf/databases/gnomad/4.1/twist_large_gnomad.vcf.bgz" + utr_file = "/ngc/projects2/dp_00005/data/dwf/databases/utr_annotation/uORF_5UTR_GRCh38_PUBLIC.txt" + alphamissense_tsv = "/ngc/projects2/dp_00005/data/dwf/databases/alphamissense/AlphaMissense_hg38.tsv.gz" + loftee_gerp_bw = "/ngc/projects2/dp_00005/data/dwf/databases/loftee/gerp_conservation_scores.homo_sapiens.GRCh38.bw" + loftee_human_ancestor= "/ngc/projects2/dp_00005/data/dwf/databases/loftee/human_ancestor.fa.gz" + loftee_sqlite = "/ngc/projects2/dp_00005/data/dwf/databases/loftee/loftee.sql" } process { - executor = 'pbs' + executor = 'local' scratch = '/scratch/nf_tmp' queueSize = 30 beforeScript = "export _JAVA_OPTIONS=-Djava.io.tmpdir=$params.cacheDir" clusterOptions = "-A $params.project -W group_list=$params.project" - // UMI mapping: withName: DOWNSAMPLE { executor = 'local' - module = ['tools','java/1.8.0','bbmap/38.89'] + // module = ['tools','java/1.8.0','bbmap/38.89'] cpus = 1 memory = 4.GB time = 2.hour } withName: UBAM { - module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] + // module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] cpus = 8 memory = { 20.GB * task.attempt } time = 4.hour } withName: ALIGNMENT_UMI { - module = ['tools', 'intel/perflibs/2020_update4','bwa-mem2/2.2.1','samtools/1.18'] + // module = ['tools', 'intel/perflibs/2020_update4','bwa-mem2/2.2.1','samtools/1.18'] cpus = 40 memory = 160.GB time = 2.hour @@ -47,7 +55,7 @@ profiles { withName: EXTRACT_UMI { executor = 'local' - module = ['tools','java/1.8.0','fgbio/1.5.1'] + // module = ['tools','java/1.8.0','fgbio/1.5.1'] cpus = 2 memory = 4.GB time = 2.hour @@ -55,14 +63,14 @@ profiles { withName: FASTQ { executor = 'local' - module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] + // module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] cpus = 2 memory = 4.GB time = 2.hour } withName: MERGEBAM { - module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] + // module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] cpus = 8 memory = { 20.GB * task.attempt } time = 4.hour @@ -70,7 +78,7 @@ profiles { withName: ALIGNMENT_METRICS { executor = 'local' - module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] + // module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] cpus = 2 memory = 8.GB time = 2.hour @@ -78,7 +86,7 @@ profiles { withName: GC_METRICS { executor = 'local' - module = ['tools','java/1.8.0','gatk/4.3.0.0','gcc/7.4.0','intel/perflibs/2020_update4','R/4.3.0'] + // module = ['tools','java/1.8.0','gatk/4.3.0.0','gcc/7.4.0','intel/perflibs/2020_update4','R/4.3.0'] cpus = 2 memory = 8.GB time = 2.hour @@ -86,7 +94,7 @@ profiles { withName: INSERT_SIZE_METRICS { executor = 'local' - module = ['tools','java/1.8.0','gatk/4.3.0.0','gcc/7.4.0','intel/perflibs/2020_update4','R/4.3.0'] + // module = ['tools','java/1.8.0','gatk/4.3.0.0','gcc/7.4.0','intel/perflibs/2020_update4','R/4.3.0'] cpus = 2 memory = 8.GB time = 2.hour @@ -94,21 +102,15 @@ profiles { withName: HS_METRICS { executor = 'local' - module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] + // module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] cpus = 2 memory = 8.GB time = 2.hour } - withName: DUPLICATE_METRICS { - module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] - cpus = 4 - memory = 16.GB - time = 2.hour - } withName: GROUP_UMI { - module = ['tools','java/1.8.0','fgbio/1.5.1'] + // module = ['tools','java/1.8.0','fgbio/1.5.1'] cpus = 8 memory = { 20.GB * task.attempt } time = 2.hour @@ -116,21 +118,21 @@ profiles { withName: UMI_METRICS { executor = 'local' - module = ['tools','java/1.8.0','fgbio/1.5.1'] + // module = ['tools','java/1.8.0','fgbio/1.5.1'] cpus = 2 memory = 8.GB time = 2.hour } withName: CALL_CONSENSUS { - module = ['tools','java/1.8.0','fgbio/1.5.1'] + // module = ['tools','java/1.8.0','fgbio/1.5.1'] cpus = 2 memory = 8.GB time = 2.hour } withName: MERGE_CONSENSUS { - module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] + // module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] cpus = 8 memory = { 20.GB * task.attempt } time = 4.hour @@ -138,7 +140,7 @@ profiles { withName: FASTQC { executor = 'local' - module = ['tools','perl/5.20.1','fastqc/0.11.8'] + // module = ['tools','perl/5.20.1','fastqc/0.11.8'] cpus = 2 memory = 4.GB time = 4.hour @@ -146,7 +148,7 @@ profiles { withName: MOSDEPTH { executor = 'local' - module = ['tools','htslib/1.16','mosdepth/0.3.3'] + // module = ['tools','htslib/1.16','mosdepth/0.3.3'] cpus = 2 memory = 4.GB time = 2.hour @@ -154,7 +156,7 @@ profiles { withName: FLAGSTAT { executor = 'local' - module = ['tools','samtools/1.18'] + // module = ['tools','samtools/1.18'] cpus = 2 memory = 3.GB time = 2.hour @@ -162,35 +164,35 @@ profiles { withName: MULTIQC { executor = 'local' - module = ['tools','multiqc/1.23'] + // module = ['tools','multiqc/1.23'] cpus = 1 memory = 4.GB time = 4.hour } withName: ALIGNMENT { - module = ['tools', 'intel/perflibs/2020_update4','bwa-mem2/2.2.1','samtools/1.18'] + // module = ['tools', 'intel/perflibs/2020_update4','bwa-mem2/2.2.1','samtools/1.18'] cpus = 40 memory = 160.GB time = 4.hour } withName: MARKDUPLICATES { - module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] + // module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] cpus = 8 memory = { 16.GB * task.attempt } time = 4.hour } withName: MARKDUPLICATES_SPARK { - module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] + // module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] cpus = 8 memory = { 16.GB * task.attempt } time = 4.hour } withName: MARKDUPLICATES_FAST { - module = ['tools','sambamba/0.8.0'] + // module = ['tools','sambamba/0.8.0'] cpus = 8 memory = { 16.GB * task.attempt } time = 4.hour @@ -198,7 +200,7 @@ profiles { withName: CLEAN { executor = 'local' - module = ['tools','samtools/1.18'] + // module = ['tools','samtools/1.18'] cpus = 2 memory = 4.GB time = 2.hour @@ -206,7 +208,7 @@ profiles { withName: ADDREADGROUP { executor = 'local' - module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] + // module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] cpus = 1 memory = 4.GB time = 2.hour @@ -214,7 +216,7 @@ profiles { withName: CRAM { executor = 'local' - module = ['tools','samtools/1.18'] + // module = ['tools','samtools/1.18'] cpus = 2 memory = 4.GB time = 2.hour @@ -222,7 +224,7 @@ profiles { withName: BAM { executor = 'local' - module = ['tools','samtools/1.18'] + // module = ['tools','samtools/1.18'] cpus = 2 memory = 4.GB time = 2.hour @@ -230,7 +232,7 @@ profiles { withName: INDEX { executor = 'local' - module = ['tools','samtools/1.18'] + // module = ['tools','samtools/1.18'] cpus = 2 memory = 4.GB time = 2.hour @@ -238,7 +240,7 @@ profiles { withName: VALIDATE { executor = 'local' - module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] + // module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] cpus = 2 memory = 4.GB time = 2.hour @@ -246,21 +248,21 @@ profiles { withName: INTERVALS { executor = 'local' - module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] + // module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] cpus = 1 memory = 2.GB time = 1.hour } withName: HAPLOTYPECALLER { - module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] + // module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] cpus = 2 memory = 8.GB time = { 5.hour * task.attempt } } withName: HAPLOTYPECALLER_SPLIT { - module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] + // module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] cpus = 2 memory = 8.GB time = 4.hour @@ -268,7 +270,7 @@ profiles { withName: DEEPVARIANT { array = 20 - module = ['tools','singularity/4.0.2'] + // module = ['tools','singularity/4.0.2'] cpus = 2 memory = 8.GB time = 4.hour @@ -276,7 +278,7 @@ profiles { withName: INDELQUAL { executor = 'local' - module = ['tools','anaconda2/4.0.0','lofreq/2.1.3.1'] + // module = ['tools','anaconda2/4.0.0','lofreq/2.1.3.1'] cpus = 2 memory = 8.GB time = 2.hour @@ -284,7 +286,7 @@ profiles { withName: LOFREQ { maxForks = 24 - module = ['tools','anaconda2/4.0.0','lofreq/2.1.3.1'] + // module = ['tools','anaconda2/4.0.0','lofreq/2.1.3.1'] cpus = 2 memory = 8.GB time = 4.hour @@ -292,7 +294,7 @@ profiles { withName: FILTER { executor = 'local' - module = ['tools','bcftools/1.16'] + // module = ['tools','bcftools/1.16'] cpus = 2 memory = 4.GB time = 2.hour @@ -300,14 +302,14 @@ profiles { withName: PILOT_PINPOINT { executor = 'local' - module = ['tools','gcc/7.4.0','intel/perflibs/2020_update4','apache-arrow/11.0.0_CPP','R/4.3.0'] + // module = ['tools','gcc/7.4.0','intel/perflibs/2020_update4','apache-arrow/11.0.0_CPP','R/4.3.0'] cpus = 2 memory = 10.GB time = 1.hour } withName: TEST { - module = ['tools','anaconda3/2023.03'] + // module = ['tools','anaconda3/2023.03'] cpus = 1 memory = 10.GB time = 1.hour @@ -315,7 +317,7 @@ profiles { withName: SUBSET { executor = 'local' - module = ['tools','samtools/1.18'] + // module = ['tools','samtools/1.18'] cpus = 2 memory = 4.GB time = 2.hour @@ -323,7 +325,7 @@ profiles { withName: BQSR { executor = 'local' - module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] + // module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] cpus = 2 memory = 8.GB time = 4.hour @@ -331,14 +333,14 @@ profiles { withName: APPLY_BQSR { executor = 'local' - module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] + // module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] cpus = 2 memory = 8.GB time = 4.hour } withName: HC_TRUTH { - module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] + // module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] cpus = 2 memory = 8.GB time = 4.hour @@ -346,7 +348,7 @@ profiles { withName: HC_TRUTH_JOINT { array = 20 - module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] + // module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] cpus = 2 memory = 8.GB time = 4.hour @@ -354,7 +356,7 @@ profiles { withName: GENOMICSDB { executor = 'local' - module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] + // module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] cpus = 40 memory = 160.GB time = 4.hour @@ -362,7 +364,7 @@ profiles { withName: GENOTYPEGVCF { executor = 'local' - module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] + // module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] cpus = 40 memory = 160.GB time = 4.hour @@ -371,7 +373,7 @@ profiles { withName: GENOTYPEGVCF_SPLIT { executor = 'local' scratch = false - module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] + // module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] cpus = 2 memory = 8.GB time = 4.hour @@ -379,7 +381,7 @@ profiles { withName: CRISP { executor = 'local' - module = ['tools','crisp/20181122'] + // module = ['tools','crisp/20181122'] cpus = 40 memory = 160.GB time = 10.hour @@ -387,7 +389,7 @@ profiles { withName: OCTOPUS { executor = 'local' - module = ['tools','singularity', 'octopus/0.7.4'] + // module = ['tools','singularity', 'octopus/0.7.4'] cpus = 2 memory = 8.GB time = 6.hour @@ -395,7 +397,7 @@ profiles { withName: MERGEVCFS { executor = 'local' - module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] + // module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] cpus = 2 memory = 8.GB time = 4.hour @@ -403,7 +405,7 @@ profiles { withName: VARTABLE { executor = 'local' - module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] + // module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] cpus = 1 memory = 4.GB time = 1.hour @@ -411,7 +413,7 @@ profiles { withName: INDEX_VCF { executor = 'local' - module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] + // module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] cpus = 1 memory = 4.GB time = 1.hour @@ -419,14 +421,14 @@ profiles { withName: NORMALISE_VCF { executor = 'local' - module = ['tools','bcftools/1.16'] + // module = ['tools','bcftools/1.16'] cpus = 2 memory = 4.GB time = 2.hour } withName: FILTER_VARIANTS { - module = ['tools','anaconda3/2023.03'] + // module = ['tools','anaconda3/2023.03'] cpus = 1 memory = 10.GB time = 1.hour @@ -434,14 +436,14 @@ profiles { withName: PINPY { executor = 'local' - module = ['tools','anaconda3/2022.10','bcftools/1.16','htslib/1.16'] + // module = ['tools','anaconda3/2022.10','bcftools/1.16','htslib/1.16'] cpus = 2 memory = 10.GB time = 2.hour } withName: SNPEFF { - module = ['tools', 'anaconda3/2022.10', 'perl/5.20.1', 'jdk/21.0.1', 'ngctools', 'snpeff/5.0e'] + // module = ['tools', 'anaconda3/2022.10', 'perl/5.20.1', 'jdk/21.0.1', 'ngctools', 'snpeff/5.0e'] executor = 'local' cpus = 2 memory = 4.GB @@ -449,7 +451,7 @@ profiles { } withName: SNPSIFT_CLINVAR { - module = ['tools', 'anaconda3/2022.10', 'perl/5.20.1', 'jdk/21.0.1', 'ngctools', 'snpeff/5.0e'] + // module = ['tools', 'anaconda3/2022.10', 'perl/5.20.1', 'jdk/21.0.1', 'ngctools', 'snpeff/5.0e'] executor = 'local' cpus = 2 memory = 4.GB @@ -457,7 +459,7 @@ profiles { } withName: SNPSIFT_FILTER { - module = ['tools', 'anaconda3/2022.10', 'perl/5.20.1', 'jdk/21.0.1', 'ngctools', 'snpeff/5.0e'] + // module = ['tools', 'anaconda3/2022.10', 'perl/5.20.1', 'jdk/21.0.1', 'ngctools', 'snpeff/5.0e'] executor = 'local' cpus = 2 memory = 4.GB @@ -466,7 +468,7 @@ profiles { withName: VARTABLE_PINS { executor = 'local' - module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] + // module = ['tools','oracle_jdk/21.0.2','gatk/4.6.0.0'] cpus = 1 memory = 4.GB time = 1.hour @@ -474,7 +476,7 @@ profiles { withName: BED_ANNOTATE { executor = 'local' - module = ['tools','bcftools/1.16'] + // module = ['tools','bcftools/1.16'] cpus = 1 memory = 4.GB time = 2.hour @@ -482,7 +484,7 @@ profiles { withName: DROP_PL { executor = 'local' - module = ['tools','bcftools/1.16'] + // module = ['tools','bcftools/1.16'] cpus = 1 memory = 4.GB time = 2.hour @@ -490,7 +492,7 @@ profiles { withName: MERGE_PINS { executor = 'local' - module = ['tools','bcftools/1.16','htslib/1.16'] + // module = ['tools','bcftools/1.16','htslib/1.16'] cpus = 1 memory = 4.GB time = 2.hour diff --git a/conf/profiles.config b/conf/profiles.config index 5ba94e8..20872df 100644 --- a/conf/profiles.config +++ b/conf/profiles.config @@ -15,14 +15,16 @@ profiles { params { outputDir = "results" - fgbio = "fgbio" annotate = false spark_markdup = false } process { executor = 'local' - + resourceLimits = [ + cpus: 5, + memory: '30.GB' + ] // Base configuration for unlabelled process cpus = { 1 * task.attempt } memory = { 6.GB * task.attempt } @@ -64,11 +66,16 @@ profiles { memory: '4.GB', time: '1.h' ] + executor = 'local' + + // Base configuration for unlabelled process + cpus = { 2 * task.attempt } + memory = { 4.GB * task.attempt } + time = { 1.h * task.attempt } } params { outputDir = "results" - fgbio = "fgbio" spark_markdup = false filter = true ploidy = 2 @@ -78,16 +85,13 @@ profiles { pooltable = "$projectDir/assets/data/test_data/pooltable.tsv" fqtable = "$projectDir/assets/data/test_data/pooltable.tsv" decodetable = "$projectDir/assets/data/test_data/decodetable.tsv" - reference_genome = "$projectDir/assets/data/reference_genomes/small/small_reference" + reference_genome = "$projectDir/assets/data/reference_genomes/small/small_reference.fna" bedfile = "$projectDir/assets/data/test_data/target_calling.bed" bedfile_bam_extraction = "$projectDir/assets/data/test_data/target_calling.bed" dictfile = "$projectDir/assets/data/test_data/target_calling.dict" snv_list = "$projectDir/assets/data/test_data/snvlist.tsv" mills = "$projectDir/assets/data/test_data/databases/mills.vcf" g1000 = "$projectDir/assets/data/test_data/databases/g1000.vcf" - clinvar_db = "$projectDir/assets/data/test_data/databases/clinvar.vcf" - snpeff_config = "$projectDir/assets/data/test_data/databases/snpeff/snpeff.config" - snpeff_cache = "$projectDir/assets/data/test_data/databases/snpeff/" } } diff --git a/conf/ssi.config b/conf/ssi.config new file mode 100644 index 0000000..6fae7ea --- /dev/null +++ b/conf/ssi.config @@ -0,0 +1,75 @@ +profiles { + + ssi { + + conda.enabled = true + + params { + outputDir = "../results" + pooltable = "pooltable.tsv" + cacheDir = "/srv/data/scratch/nf_tmp/" + reference_genome = "/srv/data/Reference/GATK_bundle/hg38_noScaffoldsMasked/Homo_sapiens_assembly38_noScaffoldsMasked.fasta" + bedfile = "/srv/data/Projects/BabyCompare/vep_optional_annotation/Target_bases_covered_by_probes_SSI_BabyPrecision-v2p1_1X_TE-95845496_hg38_ext151_merged.bed" + runHCParallel = false + max_haplotypes = 128 + snpeff_db = 'GRCh38.99' + snpeff_config = "" + snpeff_cache = "" + annotate_vep = true + vep_cache = "/srv/data/VCF_annotations/VEP/vep" + clinvar_db = "/srv/data/VCF_annotations/CLINVAR/clinvar.vcf.gz" + danmac_db = "" + blacklist_bed = "/srv/data/Projects/BabyCompare/vep_optional_annotation/hg38-blacklist.v2.sorted.bed.gz" + repeatmasker_bed = "/srv/data/Projects/BabyCompare/vep_optional_annotation/repeatmasker.sorted.bed.gz" + gnomad_vcf = "/srv/data/Projects/BabyCompare/vep_optional_annotation/baby_compare_gnomad.vcf.bgz" + utr_file = "/srv/data/Projects/BabyCompare/vep_optional_annotation/uORF_5UTR_GRCh38_PUBLIC.txt" + alphamissense_tsv = "/srv/data/Projects/BabyCompare/vep_optional_annotation/AlphaMissense_hg38.tsv.gz" + loftee_path = "/srv/data/Projects/BabyCompare/vep_optional_annotation/loftee" + loftee_gerp_bw = "/srv/data/Projects/BabyCompare/vep_optional_annotation/gerp_conservation_scores.homo_sapiens.GRCh38.bw" + loftee_human_ancestor= "/srv/data/Projects/BabyCompare/vep_optional_annotation/human_ancestor.fa.gz" + loftee_sqlite = "/srv/data/Projects/BabyCompare/vep_optional_annotation/loftee.sql" + } + + process { + executor = 'local' + beforeScript = "export _JAVA_OPTIONS=-Djava.io.tmpdir=$params.cacheDir" + + resourceLimits = [ + cpus: 12, + memory: '64.GB' + ] + + maxRetries = 2 + + // Base configuration for unlabelled process + cpus = { 2 * task.attempt } + memory = { 6.GB * task.attempt } + time = { 4.h * task.attempt } + + // Labelled configurations + withLabel:process_single { + cpus = { 1 } + memory = { 6.GB * task.attempt } + time = { 4.h * task.attempt } + } + + withLabel:process_low { + cpus = { 2 * task.attempt } + memory = { 6.GB * task.attempt } + time = { 24.h * task.attempt } + } + + withLabel:process_medium { + cpus = { 6 * task.attempt } + memory = { 36.GB * task.attempt } + time = { 8.h * task.attempt } + } + + withLabel:process_multi { + cpus = { 8 * task.attempt } + memory = { 30.GB * task.attempt } + time = { 8.h * task.attempt } + } + } + } +} diff --git a/config.json b/config.json index 189fa33..0debb96 100644 --- a/config.json +++ b/config.json @@ -1,9 +1,8 @@ { "ploidy": 48, - "step": "mapping", + "step": "all", "pooltable": "pooltable.tsv", - "decodetable": "decodetable.tsv", "outputDir": "../results", - "reference_genome": "/ngc/projects2/dp_00005/data/predisposed/resources/reference/Homo_sapiens_assembly38_noAlt", + "reference_genome": "/ngc/projects2/dp_00005/data/predisposed/resources/reference/Homo_sapiens_assembly38_noAlt.fna", "bedfile": "/ngc/projects2/dp_00005/data/predisposed/resources/targets/Probes_merged_ok_SSI_PREDiSPOSED-v1_1X_TE-93129556_hg38.bed" } \ No newline at end of file diff --git a/config_qc.json b/config_qc.json index 969eb3c..64cb398 100644 --- a/config_qc.json +++ b/config_qc.json @@ -2,9 +2,8 @@ "ploidy": 48, "step": "mapping", "downsample": 4, + "outputDir": "../results_qc_150x", "fullQC": "true", "pooltable": "pooltable.tsv", - "decodetable": "decodetable.tsv", - "reference_genome": "/ngc/projects2/dp_00005/data/predisposed/resources/reference/Homo_sapiens_assembly38_noAlt", "bedfile": "/ngc/projects2/dp_00005/data/predisposed/resources/targets/Probes_merged_ok_SSI_PREDiSPOSED-v1_1X_TE-93129556_hg38.bed" } \ No newline at end of file diff --git a/config_test.json b/config_test.json index be7454e..b805818 100644 --- a/config_test.json +++ b/config_test.json @@ -5,7 +5,7 @@ "step": "all", "pooltable": "assets/data/test_data/pooltable.tsv", "decodetable": "assets/data/test_data/decodetable.tsv", - "reference_genome": "assets/data/reference_genomes/small/small_reference", + "reference_genome": "assets/data/reference_genomes/small/small_reference.fna", "bedfile": "assets/data/test_data/target_calling.bed", "bedfile_bam_extraction": "assets/data/test_data/target_calling.bed", "dictfile": "assets/data/test_data/target_calling.dict", diff --git a/envs/ngsep/Dockerfile b/envs/ngsep/Dockerfile new file mode 100644 index 0000000..302ef38 --- /dev/null +++ b/envs/ngsep/Dockerfile @@ -0,0 +1,31 @@ +# Docker container for NGSEP +FROM eclipse-temurin:21-jre + +# Add metadata +LABEL maintainer="bioinformatics" +LABEL description="NGSEP - Next Generation Sequencing Experience Platform" +LABEL version="5.1.0" + +# Install wget for downloading +RUN apt-get update && \ + apt-get install -y wget && \ + rm -rf /var/lib/apt/lists/* + +# Download NGSEP JAR from GitHub +RUN wget -O /usr/local/bin/NGSEPcore_5.1.0.jar \ + https://github.com/NGSEP/NGSEPcore/raw/refs/heads/master/NGSEPcore_5.1.0.jar && \ + chmod 644 /usr/local/bin/NGSEPcore_5.1.0.jar && \ + ln -s /usr/local/bin/NGSEPcore_5.1.0.jar /usr/local/bin/NGSEPcore.jar + +# Create a convenience script for running NGSEP +RUN echo '#!/bin/bash\njava $JAVA_OPTS -jar /usr/local/bin/NGSEPcore.jar "$@"' > /usr/local/bin/ngsep && \ + chmod +x /usr/local/bin/ngsep + +# Create working directory +WORKDIR /data + +# Set environment variable for Java options (can be overridden) +ENV JAVA_OPTS="-Xmx4g" + +# Default command shows help +CMD ["ngsep", "--help"] \ No newline at end of file diff --git a/envs/ngsep/build b/envs/ngsep/build new file mode 100644 index 0000000..4fc9334 --- /dev/null +++ b/envs/ngsep/build @@ -0,0 +1 @@ +docker buildx build --platform linux/amd64,linux/arm64 -t madscort/ngsep:5.1.0 --push . \ No newline at end of file diff --git a/envs/ngsep/environment.yaml b/envs/ngsep/environment.yaml new file mode 100644 index 0000000..bbe76fc --- /dev/null +++ b/envs/ngsep/environment.yaml @@ -0,0 +1,7 @@ +name: ngsep +channels: + - conda-forge + - bioconda + - defaults +dependencies: + - bioconda::ngsep=4.0.1 \ No newline at end of file diff --git a/envs/picard/environment.yaml b/envs/picard/environment.yaml new file mode 100644 index 0000000..13baff6 --- /dev/null +++ b/envs/picard/environment.yaml @@ -0,0 +1,8 @@ +name: picard +channels: + - conda-forge + - bioconda + - defaults +dependencies: + - r-base=4.3.0 + - bioconda::picard=3.4.0 \ No newline at end of file diff --git a/envs/pinpy/build b/envs/pinpy/build index 3418c0b..7efe5aa 100644 --- a/envs/pinpy/build +++ b/envs/pinpy/build @@ -1 +1 @@ -docker buildx build --platform linux/amd64,linux/arm64,windows/amd64 -t madscort/dswf_pinpy:1.0 --push . \ No newline at end of file +docker buildx build --platform linux/amd64,linux/arm64 -t madscort/dswf_pinpy:1.0 --push . \ No newline at end of file diff --git a/envs/pinpy/environment.yaml b/envs/pinpy/environment.yaml index 7e89b5e..c77e65f 100644 --- a/envs/pinpy/environment.yaml +++ b/envs/pinpy/environment.yaml @@ -4,5 +4,6 @@ channels: - bioconda - defaults dependencies: - - bcftools=1.19 - - python=3.9 \ No newline at end of file + - bcftools=1.22 + - python=3.9 + - samtools=1.20 diff --git a/envs/r_env/Dockerfile b/envs/r_env/Dockerfile new file mode 100644 index 0000000..8693d39 --- /dev/null +++ b/envs/r_env/Dockerfile @@ -0,0 +1,13 @@ +FROM mambaorg/micromamba:2.3-ubuntu22.04 + +WORKDIR /work + +COPY environment.yaml /tmp/environment.yml + +USER root + +RUN micromamba install -y -n base -f /tmp/environment.yml && micromamba clean -a -y + +ENV PATH="/opt/conda/bin:${PATH}" + +CMD ["R", "--version"] \ No newline at end of file diff --git a/envs/r_env/build b/envs/r_env/build new file mode 100644 index 0000000..4b98af6 --- /dev/null +++ b/envs/r_env/build @@ -0,0 +1 @@ +docker buildx build --platform linux/amd64,linux/arm64 -t madscort/dswf_r:1.0 --push . \ No newline at end of file diff --git a/envs/r_env/environment.yaml b/envs/r_env/environment.yaml index c18dd66..2e313e2 100644 --- a/envs/r_env/environment.yaml +++ b/envs/r_env/environment.yaml @@ -1,15 +1,28 @@ name: r_env channels: - conda-forge - - bioconda - - defaults +channel_priority: strict dependencies: + # R runtime and HTML utils - r-base=4.3.0 - - r-data.table + - r-rmarkdown + - r-knitr + - pandoc + # R packages - r-dplyr=1.1.4 - r-ggplot2=3.4.4 - r-readr=2.1.5 - r-tidyr=1.3.1 - r-stringr=1.5.1 - r-purrr=1.0.2 - - r-optparse=1.7.4 \ No newline at end of file + - r-data.table + - r-optparse=1.7.4 + - r-tibble + - r-forcats + - r-scales + - r-fs + - r-here + - r-janitor + - r-ggthemes + - r-ggbeeswarm + - r-jsonlite \ No newline at end of file diff --git a/envs/rescue/Dockerfile b/envs/rescue/Dockerfile new file mode 100644 index 0000000..713996f --- /dev/null +++ b/envs/rescue/Dockerfile @@ -0,0 +1,10 @@ +FROM continuumio/miniconda3:latest + +COPY environment.yaml /tmp/environment.yaml + +RUN conda env create --name docker_env --file /tmp/environment.yaml + +ENV PATH /opt/conda/envs/docker_env/bin:$PATH +SHELL ["conda", "run", "-n", "docker_env", "/bin/bash", "-c"] + +WORKDIR /app \ No newline at end of file diff --git a/envs/rescue/build b/envs/rescue/build new file mode 100644 index 0000000..c58e67d --- /dev/null +++ b/envs/rescue/build @@ -0,0 +1 @@ +docker buildx build --platform linux/amd64,linux/arm64 -t madscort/marbl:0.6.0 --push . \ No newline at end of file diff --git a/envs/rescue/environment.yaml b/envs/rescue/environment.yaml new file mode 100644 index 0000000..8db4fb6 --- /dev/null +++ b/envs/rescue/environment.yaml @@ -0,0 +1,18 @@ +name: rescue +channels: + - conda-forge + - bioconda + - defaults +dependencies: + - python=3.11 + - pip + - numpy + - scipy>=1.16.1,<2 + - pandas>=2.3.2 + - pysam=0.23.3 + - htslib>=1.18 + - zlib + - libdeflate + - compilers + - pip: + - marbl-pool==0.6.0 \ No newline at end of file diff --git a/envs/snpsift/environment.yaml b/envs/snpsift/environment.yaml new file mode 100644 index 0000000..23c6529 --- /dev/null +++ b/envs/snpsift/environment.yaml @@ -0,0 +1,7 @@ +name: snpsift +channels: + - conda-forge + - bioconda + - defaults +dependencies: + - snpsift=5.1 \ No newline at end of file diff --git a/envs/vep/Dockerfile b/envs/vep/Dockerfile new file mode 100644 index 0000000..592482f --- /dev/null +++ b/envs/vep/Dockerfile @@ -0,0 +1,69 @@ +FROM ubuntu:20.04 AS build +ENV DEBIAN_FRONTEND=noninteractive + +# Install basic requirements +RUN apt-get update && apt-get -y install \ + wget \ + libncurses5-dev \ + libncursesw5-dev \ + libbz2-dev \ + liblzma-dev \ + build-essential \ + libz-dev \ + git \ + libncurses5 \ + libbz2-1.0 \ + liblzma5 \ + && rm -rf /var/lib/apt/lists/* + +# Install samtools +RUN wget https://github.com/samtools/samtools/releases/download/1.7/samtools-1.7.tar.bz2 && \ + tar xjvf samtools-1.7.tar.bz2 && \ + cd samtools-1.7 && \ + make && \ + make install && \ + cd .. && \ + rm -rf samtools-1.7* + +# Clone LOFTEE +RUN git clone -b grch38 https://github.com/konradjk/loftee.git /loftee + +FROM ensemblorg/ensembl-vep:release_111.0 AS runtime + +USER root + +# Set up Perl environment +ENV PERL_MM_USE_DEFAULT=1 +ENV PERL5LIB=/plugins:/opt/vep/perl5/lib/perl5:$PERL5LIB + +# Install system dependencies +RUN apt-get update && apt-get -y install \ + libncurses5 \ + libtinfo5 \ + libbz2-1.0 \ + liblzma5 \ + libsqlite3-dev \ + zlib1g \ + && rm -rf /var/lib/apt/lists/* + +# Create directory and set permissions +RUN mkdir -p /opt/vep/perl5 && \ + chmod -R 777 /opt/vep/perl5 + +# Copy samtools binary and make it executable +COPY --from=build /usr/local/bin/samtools /usr/local/bin/samtools +RUN chmod +x /usr/local/bin/samtools + +# Copy LOFTEE plugins +COPY --from=build /loftee/ /plugins + +# Install additional Perl dependencies as root +RUN cpanm --local-lib=/opt/vep/perl5 local::lib && \ + cpanm --local-lib=/opt/vep/perl5 --force DBD::SQLite Bio::Perl + +# Verify installations +RUN samtools --version && ls -la /plugins/ && \ + perl -e 'use lib q(/plugins); require q(LoF.pm); require q(utr_splice.pl); print qq(LoFTEE load OK\n)' + +# Set working directory +WORKDIR /opt/vep \ No newline at end of file diff --git a/envs/vep/build b/envs/vep/build new file mode 100644 index 0000000..1d9f1f8 --- /dev/null +++ b/envs/vep/build @@ -0,0 +1 @@ +docker buildx build --platform linux/amd64,linux/arm64 -t madscort/vep:111.0 --push . \ No newline at end of file diff --git a/envs/vep/environment.yaml b/envs/vep/environment.yaml new file mode 100644 index 0000000..5729966 --- /dev/null +++ b/envs/vep/environment.yaml @@ -0,0 +1,8 @@ +name: vep +channels: + - conda-forge + - bioconda + - defaults +dependencies: + - bioconda::ensembl-vep=114.2 + - bioconda::samtools=1.20 diff --git a/lib/Utils.groovy b/lib/Utils.groovy new file mode 100644 index 0000000..b80041e --- /dev/null +++ b/lib/Utils.groovy @@ -0,0 +1,17 @@ +import nextflow.Channel +import static nextflow.Nextflow.* + +class Utils { + /** + * Return a value channel carrying (mainFile, indexOrEmpty). + * indexExt is a single extension like '.bai' or '.tbi'. + */ + static def indexedFileChannel(String filePath, String indexExt = '.tbi') { + if (!filePath) return Channel.value(tuple([], [])) + + def main = file(filePath) + def index = file("${filePath}${indexExt}") + + return Channel.value(index.exists() ? tuple(main, index) : tuple(main, [])) + } +} \ No newline at end of file diff --git a/main.nf b/main.nf index 1037035..550c00e 100644 --- a/main.nf +++ b/main.nf @@ -5,13 +5,6 @@ // Use newest nextflow dsl nextflow.enable.dsl = 2 -log.info """\ - =================================== - D o B S e q - W F - =================================== - """ - .stripIndent() - /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ IMPORT MODULES AND SUBWORKFLOWS @@ -30,9 +23,16 @@ include { ANNOTATION } from './subworkflows/annotation' include { BAM } from './modules/bam' include { INDEX } from './modules/index' +// Method for creating decode table and matrix context +include { MATRIX_CONTEXT } from './modules/matrix_context' + // Pinpoint methods include { PILOT_PINPOINT } from './modules/pilot_pinpoint' -include { PINPOINT } from './subworkflows/pinpoint' +include { PINPOINT_VCF } from './subworkflows/pinpoint_vcf' +include { PINPOINT_TAB } from './subworkflows/pinpoint_tab' + +include { VARIANT_RESCUE } from './subworkflows/variant_rescue' +include { REPORT } from './modules/report' // Test include { TEST } from './modules/test' @@ -47,9 +47,9 @@ include { TEST } from './modules/test' pooltable_ch = Channel .fromPath(params.pooltable) .splitCsv(sep: '\t') - .map { row -> tuple(row[0], [file(row[1]), file(row[2])]) } + .map { row -> tuple(row[0], [file(row[2]), file(row[3])]) } -if (params.step == 'calling') { +if (params.step == 'calling' || (params.step == 'pinpoint' && params.variant_rescue)) { cramtable_ch = Channel .fromPath(params.cramtable) .splitCsv(sep: '\t') @@ -63,13 +63,21 @@ if (params.step == 'pinpoint') { .map { row -> tuple(row[0], file(row[1]), row[2]) } } +rescue_probabilities = Channel.value([]) +report_rmd = Channel.fromPath("$projectDir/bin/report.Rmd", checkIfExists: true) + /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ REFERENCE FILE CHANNELS ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */ -reference_genome_ch = Channel.fromPath(params.reference_genome + "*", checkIfExists: true).collect() +ref_file = file(params.reference_genome) +ref_base = ref_file.parent.resolve(ref_file.baseName) +reference_genome_ch = Channel + .fromPath(["${ref_file}*", "${ref_base}.*"], checkIfExists: true) + .unique { it.name } + .collect() bedfile_ch = Channel.fromPath(params.bedfile, checkIfExists: true).collect() /* @@ -79,13 +87,22 @@ bedfile_ch = Channel.fromPath(params.bedfile, checkIfExists: true).collect() */ workflow { + + log.info """\ + =================================== + D o B S e q - W F + =================================== + """ + .stripIndent() + + if (params.step == 'mapping' || params.step == 'all' || params.step == '') { if (params.umi) { bam_file_w_index_ch = MAPPING_UMI(pooltable_ch, reference_genome_ch, bedfile_ch) } else { bam_file_w_index_ch = MAPPING(pooltable_ch, reference_genome_ch, bedfile_ch) } - } else if (params.step == 'calling') { + } else if (params.step == 'calling' || (params.step == 'pinpoint' && params.variant_rescue)) { BAM(cramtable_ch, reference_genome_ch) bam_file_w_index_ch = INDEX(BAM.out.bam_file) } @@ -110,16 +127,45 @@ workflow { } if (params.step == 'pinpoint' || params.step == 'all' || params.step == '') { + pooltable = file(params.pooltable) if (params.pinpoint_method == 'pilot') { vcf_ch = gatk_ch .mix(lofreq_ch) .map { pool_id, vcf_file -> vcf_file } - PILOT_PINPOINT(vcf_ch.collect(), file(params.pooltable), file(params.decodetable)) + PILOT_PINPOINT(vcf_ch.collect(), pooltable, file(params.decodetable)) pin_ch = PILOT_PINPOINT.out.pinned_variants } else if (params.pinpoint_method == 'new') { - PINPOINT(gatk_ch, file(params.decodetable), reference_genome_ch) - pin_ch = PINPOINT.out.pinned_variants + if (params.decodetable == "") { + MATRIX_CONTEXT(pooltable,[]) + matrix_context = MATRIX_CONTEXT.out.json + decode_table = MATRIX_CONTEXT.out.decodetable + } else { + MATRIX_CONTEXT(pooltable,file(params.decodetable)) + matrix_context = MATRIX_CONTEXT.out.json + decode_table = file(params.decodetable) + } + PINPOINT_VCF(gatk_ch, pooltable, decode_table, matrix_context, reference_genome_ch) + pin_ch = PINPOINT_VCF.out.pinned_variants + + if (params.pinpoint_tab) { + PINPOINT_TAB(gatk_ch, matrix_context, reference_genome_ch) + + if (params.variant_rescue) { + VARIANT_RESCUE(gatk_ch, bam_file_w_index_ch, pooltable, decode_table, reference_genome_ch, bedfile_ch) + rescue_probabilities = VARIANT_RESCUE.out.rescue_probabilities + } + + if (params.pinpoint_report) { + REPORT( + report_rmd, + matrix_context, + PINPOINT_TAB.out.pool_variants, + PINPOINT_TAB.out.pinpoint_variants, + PINPOINT_TAB.out.annotations, + rescue_probabilities) + } + } } } diff --git a/modules/addreadgroup.nf b/modules/addreadgroup.nf index f6c8889..d218d3a 100644 --- a/modules/addreadgroup.nf +++ b/modules/addreadgroup.nf @@ -3,7 +3,7 @@ process ADDREADGROUP { tag "Add read group - $sample_id" conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + container workflow.containerEngine == 'singularity' ? params.container.singularity.gatk : params.container.docker.gatk input: tuple val(sample_id), path(bam_file, stageAs: "raw/*") diff --git a/modules/alignment.nf b/modules/alignment.nf index b2649c4..57f3ee4 100644 --- a/modules/alignment.nf +++ b/modules/alignment.nf @@ -4,8 +4,8 @@ process ALIGNMENT { // Align reads to the reference genome using BWA, convert to BAM and sort. conda "$projectDir/envs/bwa/environment.yaml" - container params.container.bwa - + container workflow.containerEngine == 'singularity' ? params.container.singularity.bwa : params.container.docker.bwa + publishDir "${params.outputDir}/log/mapping/", pattern: "${sample_id}.log", mode:'copy' input: @@ -17,7 +17,7 @@ process ALIGNMENT { path "${sample_id}.log" script: - def db = file(params.reference_genome).getName() + ".fna" + def db = file(params.reference_genome).name // int bwa_cpus = Math.max((task.cpus * 0.75) as int, 1) // int samtools_cpus = Math.max(task.cpus - bwa_cpus, 1) diff --git a/modules/alignment_metrics.nf b/modules/alignment_metrics.nf index f642b22..9c79cfa 100644 --- a/modules/alignment_metrics.nf +++ b/modules/alignment_metrics.nf @@ -4,7 +4,7 @@ process ALIGNMENT_METRICS { // Collect alignment metrics for bam file conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + container workflow.containerEngine == 'singularity' ? params.container.singularity.gatk : params.container.docker.gatk publishDir "${params.outputDir}/log/alignment_metrics/", pattern: "${sample_id}*.align.txt", mode:'copy' @@ -17,7 +17,7 @@ process ALIGNMENT_METRICS { path("${sample_id}*.align.txt"), emit: metrics_file script: - def db = file(params.reference_genome).getName() + ".fna" + def db = file(params.reference_genome).name def log_filename = log_suffix == "" ? "${sample_id}.align.txt" : "${sample_id}_${log_suffix}.align.txt" """ gatk CollectAlignmentSummaryMetrics \ diff --git a/modules/alignment_umi.nf b/modules/alignment_umi.nf index 3809349..7e4e47a 100644 --- a/modules/alignment_umi.nf +++ b/modules/alignment_umi.nf @@ -4,6 +4,7 @@ process ALIGNMENT_UMI { // Align reads to the reference genome using BWA, convert to BAM and sort by QUERY NAME. conda "$projectDir/envs/bwa_umi/environment.yaml" + container workflow.containerEngine == 'singularity' ? params.container.singularity.fgbio : params.container.docker.fgbio publishDir "${params.outputDir}/log/mapping/", pattern: "${sample_id}.log", mode:'copy' @@ -16,7 +17,7 @@ process ALIGNMENT_UMI { path "${sample_id}.log" script: - def db = file(params.reference_genome).getName() + ".fna" + def db = file(params.reference_genome).name // int bwa_cpus = Math.max((task.cpus * 0.75) as int, 1) // int samtools_cpus = Math.max(task.cpus - bwa_cpus, 1) diff --git a/modules/apply_bqsr.nf b/modules/apply_bqsr.nf index 707e7ff..e5b12d8 100644 --- a/modules/apply_bqsr.nf +++ b/modules/apply_bqsr.nf @@ -3,7 +3,7 @@ process APPLY_BQSR { tag "Apply score recalibration - $sample_id" conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + container workflow.containerEngine == 'singularity' ? params.container.singularity.gatk : params.container.docker.gatk input: tuple val(sample_id), path(bam_file, stageAs: "raw/*"), path(bqsr_table) @@ -13,7 +13,7 @@ process APPLY_BQSR { tuple val(sample_id), path("${sample_id}.bam"), emit: corrected_bam_file script: - def db = file(params.reference_genome).getName() + ".fna" + def db = file(params.reference_genome).name """ gatk ApplyBQSR \ -R ${db} \ diff --git a/modules/bam.nf b/modules/bam.nf index f34b97e..67b61e4 100644 --- a/modules/bam.nf +++ b/modules/bam.nf @@ -3,7 +3,7 @@ process BAM { tag "CRAM->BAM - $sample_id" conda "$projectDir/envs/samtools/environment.yaml" - container params.container.samtools + container workflow.containerEngine == 'singularity' ? params.container.singularity.samtools : params.container.docker.samtools input: tuple val(sample_id), path(cram_file) @@ -13,7 +13,7 @@ process BAM { tuple val(sample_id), path("${sample_id}.bam"), emit: bam_file script: - def db = file(params.reference_genome).getName() + ".fna" + def db = file(params.reference_genome).name """ samtools view \ -@ ${task.cpus-1} \ diff --git a/modules/bed_annotate.nf b/modules/bed_annotate.nf index a3d323a..bc3e0a9 100644 --- a/modules/bed_annotate.nf +++ b/modules/bed_annotate.nf @@ -4,7 +4,7 @@ process BED_ANNOTATE { // Annotate VCF by third column in BED file conda "$projectDir/envs/bcftools/environment.yaml" - container params.container.bcftools + container workflow.containerEngine == 'singularity' ? params.container.singularity.bcftools : params.container.docker.bcftools publishDir "${params.outputDir}/log/bed_annotate/${sample_id}/", pattern: "${sample_id}.${caller}.log", mode:'copy' diff --git a/modules/bgzip.nf b/modules/bgzip.nf new file mode 100644 index 0000000..c128374 --- /dev/null +++ b/modules/bgzip.nf @@ -0,0 +1,28 @@ +process BGZIP { + label 'process_low' + tag "VCF->VCF.gz - $sample_id" + + conda "$projectDir/envs/samtools/environment.yaml" + container workflow.containerEngine == 'singularity' ? params.container.singularity.samtools : params.container.docker.samtools + + publishDir "${params.outputDir}/norm_bgzip_idx/", pattern: "*vcf.*", mode:'copy' + + input: + tuple val(sample_id), path(vcf_file) + + output: + tuple val(sample_id), path("${vcf_file}.gz"), path("${vcf_file}.gz.tbi"), emit: bgzf_file + + script: + """ + bgzip ${vcf_file} + tabix -p vcf ${vcf_file}.gz + """ + + stub: + """ + touch "${vcf_file}.gz" + touch "${vcf_file}.gz.tbi" + """ +} + diff --git a/modules/bqsr.nf b/modules/bqsr.nf index 800960d..ffc237a 100644 --- a/modules/bqsr.nf +++ b/modules/bqsr.nf @@ -3,7 +3,7 @@ process BQSR { tag "Calculate score recalibration - $sample_id" conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + container workflow.containerEngine == 'singularity' ? params.container.singularity.gatk : params.container.docker.gatk input: tuple val(sample_id), path(bam_file) @@ -15,7 +15,7 @@ process BQSR { tuple val(sample_id), path(bam_file), path("${sample_id}.BQSR.table"), emit: bqsr_file script: - def db = file(params.reference_genome).getName() + ".fna" + def db = file(params.reference_genome).name def mills_db = file(params.mills).getName() def g1000_db = file(params.g1000).getName() """ diff --git a/modules/call_consensus.nf b/modules/call_consensus.nf index 01217bc..32d054e 100644 --- a/modules/call_consensus.nf +++ b/modules/call_consensus.nf @@ -3,7 +3,7 @@ process CALL_CONSENSUS { tag "Call consensus reads - $sample_id" conda "$projectDir/envs/fgbio/environment.yaml" - container params.container.fgbio + container workflow.containerEngine == 'singularity' ? params.container.singularity.fgbio : params.container.docker.fgbio input: tuple val(sample_id), path(bam_file) @@ -13,7 +13,7 @@ process CALL_CONSENSUS { script: """ - ${params.fgbio} --tmp-dir=. --compression 1 --async-io CallDuplexConsensusReads \ + fgbio --tmp-dir=. --compression 1 --async-io CallDuplexConsensusReads \ --input=${bam_file} \ --output="${sample_id}.ubam" \ --error-rate-pre-umi=45 \ diff --git a/modules/clean.nf b/modules/clean.nf index 6bbbe08..bdb4435 100644 --- a/modules/clean.nf +++ b/modules/clean.nf @@ -3,7 +3,7 @@ process CLEAN { tag "Clean bam file - $sample_id" conda "$projectDir/envs/samtools/environment.yaml" - container params.container.samtools + container workflow.containerEngine == 'singularity' ? params.container.singularity.samtools : params.container.docker.samtools input: tuple val(sample_id), path(bam_file) diff --git a/modules/combinegvcfs.nf b/modules/combinegvcfs.nf index 538c3ce..df4ac8b 100644 --- a/modules/combinegvcfs.nf +++ b/modules/combinegvcfs.nf @@ -4,7 +4,7 @@ process COMBINEGVCFS { // Merge gVCF files into a single GenomicsDB conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + container workflow.containerEngine == 'singularity' ? params.container.singularity.gatk : params.container.docker.gatk publishDir "${params.outputDir}/log/", pattern: "genomicsdb.log", mode:'copy' publishDir "${params.outputDir}/splits/", pattern: "cohort.${interval}.g.vcf.gz", mode:'copy' @@ -21,7 +21,7 @@ process COMBINEGVCFS { path "genomicsdb.log" script: - db = file(params.reference_genome).getName() + ".fna" + db = file(params.reference_genome).name input_files_command = vcfs.collect(){"--variant ${it}"}.join(' ') avail_mem = (task.memory.mega*0.8).intValue() """ diff --git a/modules/cram.nf b/modules/cram.nf index 9908cd8..c984a09 100644 --- a/modules/cram.nf +++ b/modules/cram.nf @@ -3,7 +3,7 @@ process CRAM { tag "BAM->CRAM - $sample_id" conda "$projectDir/envs/samtools/environment.yaml" - container params.container.samtools + container workflow.containerEngine == 'singularity' ? params.container.singularity.samtools : params.container.docker.samtools publishDir "${params.outputDir}/cram/", pattern: "${sample_id}.cram", mode:'copy' @@ -16,7 +16,7 @@ process CRAM { path("${sample_id}_cram.tsv"), emit: cram_info script: - def db = file(params.reference_genome).getName() + ".fna" + def db = file(params.reference_genome).name def publishDir = file(params.outputDir + "/cram/" + sample_id + ".cram") """ samtools view \ diff --git a/modules/crisp.nf b/modules/crisp.nf index 9103112..15c014f 100644 --- a/modules/crisp.nf +++ b/modules/crisp.nf @@ -22,7 +22,7 @@ process CRISP { script: def input_files_command = bam_files.collect(){"--bam ${it}"}.join(' ') - def db = file(params.reference_genome).getName() + ".fna" + def db = file(params.reference_genome).name """ CRISP \ --ref ${db} \ diff --git a/modules/deepvariant.nf b/modules/deepvariant.nf index 76e326b..9a04d01 100644 --- a/modules/deepvariant.nf +++ b/modules/deepvariant.nf @@ -7,7 +7,7 @@ process DEEPVARIANT { // This version works on NGC only. To run locally, comment out "singularity exec" and use "run_deepvariant" directly. conda "$projectDir/envs/deepvar/environment.yaml" - container params.container.deepvariant + container workflow.containerEngine == 'singularity' ? params.container.singularity.deepvariant : params.container.docker.deepvariant publishDir "${params.outputDir}/log/deepvariant/", pattern: "${sample_id}.DV.log", mode:'copy' publishDir "${params.outputDir}/variants/", pattern: "${sample_id}.DV.vcf.gz", mode:'copy' @@ -22,12 +22,11 @@ process DEEPVARIANT { path "${sample_id}.DV.log" script: - def db = file(params.reference_genome).getName() + ".fna" + def db = file(params.reference_genome).name def ref_dir = file(params.reference_genome).getParent() def target_dir = file(params.bedfile).getParent() """ - # singularity exec --bind ${PWD} --bind ${ref_dir} --bind ${target_dir} /services/tools/deepvariant/1.5.0/deepvariant_1.5.0.sif \ - run_deepvariant \ + /opt/deepvariant/bin/run_deepvariant \ --ref=${db} \ --reads=$bam_file \ --output_vcf=${sample_id}.DV.vcf.gz \ diff --git a/modules/downsample.nf b/modules/downsample.nf index f5f7162..87ce742 100644 --- a/modules/downsample.nf +++ b/modules/downsample.nf @@ -9,7 +9,7 @@ process DOWNSAMPLE { // Approximately 40% is on target - so cap set at 400M reads. 200M per file. conda "$projectDir/envs/bbmap/environment.yaml" - container params.container.bbmap + container workflow.containerEngine == 'singularity' ? params.container.singularity.bbmap : params.container.docker.bbmap input: tuple val(sample_id), path(reads, stageAs: 'raw/*') diff --git a/modules/drop_pl.nf b/modules/drop_pl.nf index 086810b..fdfb91d 100644 --- a/modules/drop_pl.nf +++ b/modules/drop_pl.nf @@ -4,7 +4,7 @@ process DROP_PL { // Drop PL tag from VCF file. conda "$projectDir/envs/bcftools/environment.yaml" - container params.container.bcftools + container workflow.containerEngine == 'singularity' ? params.container.singularity.bcftools : params.container.docker.bcftools input: tuple val(sample_id), path(vcf_file, stageAs: "input/*"), path(index, stageAs: "input/*") diff --git a/modules/duplicate_metrics.nf b/modules/duplicate_metrics.nf deleted file mode 100644 index 9d0bfa8..0000000 --- a/modules/duplicate_metrics.nf +++ /dev/null @@ -1,35 +0,0 @@ -process DUPLICATE_METRICS { - label 'process_medium' - tag "$sample_id" - // Mark duplicate reads in BAM files - - conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk - - publishDir "${params.outputDir}/log/duplicate_metrics/", pattern: "${sample_id}*.txt", mode:'copy' - - - input: - tuple val(sample_id), path(bam_file) - val log_suffix - - output: - path "${sample_id}*.txt", emit: metrics_file - - script: - def log_filename = log_suffix == "" ? "${sample_id}.dup.txt" : "${sample_id}_${log_suffix}.dup.txt" - """ - gatk MarkDuplicates \ - --I ${bam_file} \ - --M ${log_filename} \ - --O ${sample_id}.marked.bam - - rm ${sample_id}.marked.bam - """ - - stub: - def log_filename = log_suffix == "" ? "${sample_id}.dup.txt" : "${sample_id}_${log_suffix}.dup.txt" - """ - touch "${log_filename}" - """ -} \ No newline at end of file diff --git a/modules/extract_umi.nf b/modules/extract_umi.nf index b974fd8..9ad4f89 100644 --- a/modules/extract_umi.nf +++ b/modules/extract_umi.nf @@ -10,7 +10,7 @@ process EXTRACT_UMI { // Twist UMI adapters are 5 base pairs with a 2 base pair skip, resulting in the read structure 5M2S+T conda "$projectDir/envs/fgbio/environment.yaml" - container params.container.fgbio + container workflow.containerEngine == 'singularity' ? params.container.singularity.fgbio : params.container.docker.fgbio input: tuple val(sample_id), path(bam_file, stageAs: "raw/*") @@ -20,7 +20,7 @@ process EXTRACT_UMI { script: """ - ${params.fgbio} --tmp-dir=. --compression 1 --async-io ExtractUmisFromBam \ + fgbio --tmp-dir=. --compression 1 --async-io ExtractUmisFromBam \ --input=${bam_file} \ --output=${sample_id}.ubam \ --read-structure=5M2S+T 5M2S+T \ diff --git a/modules/fastq.nf b/modules/fastq.nf index 1566c70..2631777 100644 --- a/modules/fastq.nf +++ b/modules/fastq.nf @@ -4,7 +4,7 @@ process FASTQ { // Convert unaligned uBAM files to FastQ. conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + container workflow.containerEngine == 'singularity' ? params.container.singularity.gatk : params.container.docker.gatk input: tuple val(sample_id), path(bam_file) diff --git a/modules/fastqc.nf b/modules/fastqc.nf index fb8bc94..e800808 100644 --- a/modules/fastqc.nf +++ b/modules/fastqc.nf @@ -3,7 +3,7 @@ process FASTQC { tag "FastQC - ${sample_id}" conda "$projectDir/envs/fastqc/environment.yaml" - container params.container.fastqc + container workflow.containerEngine == 'singularity' ? params.container.singularity.fastqc : params.container.docker.fastqc publishDir "${params.outputDir}/log/fastqc/", pattern: "${sample_id}_fastqc.html", mode:'copy' publishDir "${params.outputDir}/log/fastqc/", pattern: "${sample_id}_fastqc.zip", mode:'copy' diff --git a/modules/filter.nf b/modules/filter.nf index 6bffa31..a9d0aea 100644 --- a/modules/filter.nf +++ b/modules/filter.nf @@ -1,7 +1,7 @@ process FILTER { label 'process_low' conda "$projectDir/envs/bcftools/environment.yaml" - container params.container.bcftools + container workflow.containerEngine == 'singularity' ? params.container.singularity.bcftools : params.container.docker.bcftools publishDir "${params.outputDir}/variants/filtered/", pattern: "${sample_id}.lofreq.vcf.gz", mode:'copy' diff --git a/modules/filter_variants.nf b/modules/filter_variants.nf index c92a19a..0cbbadb 100644 --- a/modules/filter_variants.nf +++ b/modules/filter_variants.nf @@ -1,7 +1,7 @@ process FILTER_VARIANTS { label 'process_low' conda "$projectDir/envs/filter_variants/environment.yaml" - container params.container.filter_variants + container workflow.containerEngine == 'singularity' ? params.container.singularity.filter_variants : params.container.docker.filter_variants publishDir "${params.outputDir}/filtered_variants/", pattern: "*.vcf", mode:'copy' diff --git a/modules/flagstat.nf b/modules/flagstat.nf index 50ca95f..9adbe99 100644 --- a/modules/flagstat.nf +++ b/modules/flagstat.nf @@ -3,7 +3,7 @@ process FLAGSTAT { tag "Flagstat - $sample_id" conda "$projectDir/envs/samtools/environment.yaml" - container params.container.samtools + container workflow.containerEngine == 'singularity' ? params.container.singularity.samtools : params.container.docker.samtools publishDir "${params.outputDir}/log/flagstat/", pattern: "${sample_id}*.flagstat", mode:'copy' diff --git a/modules/freebayes.nf b/modules/freebayes.nf index b0414c9..5cea51c 100644 --- a/modules/freebayes.nf +++ b/modules/freebayes.nf @@ -4,7 +4,7 @@ process FREEBAYES { // Call variants using Lofreq - call parallel conda "$projectDir/envs/freebayes/environment.yaml" - container params.container.freebayes + container workflow.containerEngine == 'singularity' ? params.container.singularity.freebayes : params.container.docker.freebayes publishDir "${params.outputDir}/log/freebayes/", pattern: "${sample_id}.freebayes.log", mode:'copy' publishDir "${params.outputDir}/variants/", pattern: "${sample_id}.freebayes.vcf.gz", mode:'copy' @@ -20,7 +20,7 @@ process FREEBAYES { path "${sample_id}.freebayes.log" script: - def db = file(params.reference_genome).getName() + ".fna" + def db = file(params.reference_genome).name def publishDir = file(params.outputDir + "/variants/" + sample_id + ".freebayes.vcf.gz") """ freebayes \ diff --git a/modules/gc_metrics.nf b/modules/gc_metrics.nf index 710d72b..36cf28d 100644 --- a/modules/gc_metrics.nf +++ b/modules/gc_metrics.nf @@ -3,8 +3,8 @@ process GC_METRICS { tag "$sample_id" // Collect gc bias metrics for bam file - conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + conda "$projectDir/envs/picard/environment.yaml" + container workflow.containerEngine == 'singularity' ? params.container.singularity.picard : params.container.docker.picard publishDir "${params.outputDir}/log/gc_bias_metrics/", pattern: "${sample_id}*.gc.txt", mode:'copy' publishDir "${params.outputDir}/log/gc_bias_metrics/", pattern: "${sample_id}*.gc.summary.txt", mode:'copy' @@ -20,11 +20,11 @@ process GC_METRICS { path("${sample_id}*.gc.summary.txt"), emit: summary_file script: - def db = file(params.reference_genome).getName() + ".fna" + def db = file(params.reference_genome).name def log_filename = log_suffix == "" ? "${sample_id}.gc.txt" : "${sample_id}_${log_suffix}.gc.txt" def log_filename_summary = log_suffix == "" ? "${sample_id}.gc.summary.txt" : "${sample_id}_${log_suffix}.gc.summary.txt" """ - gatk CollectGcBiasMetrics \ + picard CollectGcBiasMetrics \ -I ${bam_file} \ -O ${log_filename} \ -S ${log_filename_summary} \ diff --git a/modules/genomicsdb.nf b/modules/genomicsdb.nf index 975d78d..5760428 100644 --- a/modules/genomicsdb.nf +++ b/modules/genomicsdb.nf @@ -4,7 +4,7 @@ process GENOMICSDB { // Merge gVCF files into a single GenomicsDB conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + container workflow.containerEngine == 'singularity' ? params.container.singularity.gatk : params.container.docker.gatk publishDir "${params.outputDir}/log/", pattern: "genomicsdb.log", mode:'copy' publishDir "${params.outputDir}/", pattern: "./genomicsdb", mode:'copy' @@ -19,9 +19,10 @@ process GENOMICSDB { path "genomicsdb.log" script: + def avail_mem = (task.memory.mega*0.8).intValue() input_files_command = vcfs.collect(){"--variant ${it}"}.join(' ') """ - gatk --java-options "-Xmx8g -XX:-UsePerfData" \ + gatk --java-options "-Xmx${avail_mem}M -XX:-UsePerfData" \ GenomicsDBImport \ $input_files_command \ --genomicsdb-workspace-path genomicsdb \ diff --git a/modules/genotypegvcf.nf b/modules/genotypegvcf.nf index 342821f..29971f8 100644 --- a/modules/genotypegvcf.nf +++ b/modules/genotypegvcf.nf @@ -4,7 +4,7 @@ process GENOTYPEGVCF { // Call variants on multisample gVCF db conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + container workflow.containerEngine == 'singularity' ? params.container.singularity.gatk : params.container.docker.gatk publishDir "${params.outputDir}/log/genotypegvcf/", pattern: "${sample_id}.genotypegvcf.log", mode:'copy' publishDir "${params.outputDir}/variants/", pattern: "${sample_id}.GATK.vcf.gz", mode:'copy' @@ -19,10 +19,11 @@ process GENOTYPEGVCF { path "${sample_id}.genotypegvcf.log" script: - def db = file(params.reference_genome).getName() + ".fna" + def db = file(params.reference_genome).name def publishDir = file(params.outputDir + "/variants/" + sample_id + ".GATK.vcf.gz") + def avail_mem = (task.memory.mega*0.8).intValue() """ - gatk --java-options "-Xmx8g -XX:-UsePerfData" \ + gatk --java-options "-Xmx${avail_mem}M -XX:-UsePerfData" \ GenotypeGVCFs \ -R ${db} \ -V ${gvcf_file} \ diff --git a/modules/genotypegvcf_split.nf b/modules/genotypegvcf_split.nf index c0c760b..d43d82e 100644 --- a/modules/genotypegvcf_split.nf +++ b/modules/genotypegvcf_split.nf @@ -4,7 +4,7 @@ process GENOTYPEGVCF_SPLIT { // Call variants on multisample gVCF db conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + container workflow.containerEngine == 'singularity' ? params.container.singularity.gatk : params.container.docker.gatk publishDir "${params.outputDir}/log/genotypegvcf/${sample_id}/", pattern: "${sample_id}*.genotypegvcf.log", mode:'copy' @@ -17,10 +17,11 @@ process GENOTYPEGVCF_SPLIT { path "${sample_id}*.genotypegvcf.log" script: - def db = file(params.reference_genome).getName() + ".fna" + def db = file(params.reference_genome).name def id = interval ? "${sample_id}.${interval}" : sample_id + def avail_mem = (task.memory.mega*0.8).intValue() """ - gatk --java-options "-Xmx8g -XX:-UsePerfData" \ + gatk --java-options "-Xmx${avail_mem}M -XX:-UsePerfData" \ GenotypeGVCFs \ -R ${db} \ -V ${vcf_file} \ diff --git a/modules/group_umi.nf b/modules/group_umi.nf index a68e542..eb01406 100644 --- a/modules/group_umi.nf +++ b/modules/group_umi.nf @@ -3,7 +3,7 @@ process GROUP_UMI { tag "Group reads by UMI - $sample_id" conda "$projectDir/envs/fgbio/environment.yaml" - container params.container.fgbio + container workflow.containerEngine == 'singularity' ? params.container.singularity.fgbio : params.container.docker.fgbio publishDir "${params.outputDir}/log/group_umi/", pattern: "${sample_id}.family_size_histogram.txt", mode:'copy' @@ -16,7 +16,7 @@ process GROUP_UMI { script: """ - ${params.fgbio} --tmp-dir=. --compression 1 --async-io GroupReadsByUmi \ + fgbio --tmp-dir=. --compression 1 --async-io GroupReadsByUmi \ --strategy=paired \ --input=${bam_file} \ --output="${sample_id}.bam" \ diff --git a/modules/haplotypecaller.nf b/modules/haplotypecaller.nf index cdde1de..77881c0 100644 --- a/modules/haplotypecaller.nf +++ b/modules/haplotypecaller.nf @@ -4,7 +4,7 @@ process HAPLOTYPECALLER { // Call variants using GATK - HaplotypeCaller conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + container workflow.containerEngine == 'singularity' ? params.container.singularity.gatk : params.container.docker.gatk publishDir "${params.outputDir}/log/haplotypecaller/", pattern: "${sample_id}.g.haplotypecaller.log", mode:'copy' publishDir "${params.outputDir}/variants/", pattern: "${sample_id}.GATK.g.vcf.gz", mode:'copy' @@ -19,9 +19,10 @@ process HAPLOTYPECALLER { path "${sample_id}.g.haplotypecaller.log" script: - def db = file(params.reference_genome).getName() + ".fna" + def db = file(params.reference_genome).name + def avail_mem = (task.memory.mega*0.8).intValue() """ - gatk --java-options "-Xmx8g -XX:-UsePerfData" HaplotypeCaller \ + gatk --java-options "-Xmx${avail_mem}M -XX:-UsePerfData" HaplotypeCaller \ -R ${db} \ -I ${bam_file} \ -L ${bedfile} \ @@ -42,7 +43,7 @@ process HAPLOTYPECALLER { -A AlleleFraction \ --disable-read-filter NotDuplicateReadFilter \ --max-alternate-alleles 3 \ - --max-num-haplotypes-in-population 1000 \ + --max-num-haplotypes-in-population ${params.max_haplotypes} \ --max-reads-per-alignment-start 0 \ --create-output-variant-index \ --tmp-dir . \ diff --git a/modules/haplotypecaller_split.nf b/modules/haplotypecaller_split.nf index df310a3..f46cf0c 100644 --- a/modules/haplotypecaller_split.nf +++ b/modules/haplotypecaller_split.nf @@ -4,7 +4,7 @@ process HAPLOTYPECALLER_SPLIT { // Call variants using GATK - HaplotypeCaller conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + container workflow.containerEngine == 'singularity' ? params.container.singularity.gatk : params.container.docker.gatk publishDir "${params.outputDir}/log/haplotypecaller/${sample_id}/", pattern: "${sample_id}*.g.haplotypecaller.log", mode:'copy' publishDir "${params.outputDir}/variants/gvcf/${sample_id}/", pattern: "${sample_id}*.GATK.g.vcf.gz", mode:'copy' @@ -19,12 +19,13 @@ process HAPLOTYPECALLER_SPLIT { path "${sample_id}*.g.haplotypecaller.log" script: - def db = file(params.reference_genome).getName() + ".fna" + def db = file(params.reference_genome).name def id = interval ? "${sample_id}.${interval}" : sample_id def add_intervals = interval ? "-L ${interval}" : "" def add_intersection = interval ? "--interval-set-rule INTERSECTION" : "" + def avail_mem = (task.memory.mega*0.8).intValue() """ - gatk --java-options "-Xmx8g -XX:-UsePerfData" HaplotypeCaller \ + gatk --java-options "-Xmx${avail_mem}M -XX:-UsePerfData" HaplotypeCaller \ -R ${db} \ -I ${bam_file} \ -L ${interval_list} \ diff --git a/modules/haplotypecaller_truth.nf b/modules/haplotypecaller_truth.nf index 73b5808..0d5f4a9 100644 --- a/modules/haplotypecaller_truth.nf +++ b/modules/haplotypecaller_truth.nf @@ -4,7 +4,7 @@ process HC_TRUTH { // Call variants using GATK - HaplotypeCaller conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + container workflow.containerEngine == 'singularity' ? params.container.singularity.gatk : params.container.docker.gatk publishDir "${params.outputDir}/log/haplotypecaller/", pattern: "${sample_id}.haplotypecaller.log", mode:'copy' publishDir "${params.outputDir}/variants/", pattern: "${sample_id}.GATK.vcf.gz", mode:'copy' @@ -19,9 +19,10 @@ process HC_TRUTH { path "${sample_id}.haplotypecaller.log" script: - def db = file(params.reference_genome).getName() + ".fna" + def db = file(params.reference_genome).name + def avail_mem = (task.memory.mega*0.8).intValue() """ - gatk --java-options "-Xmx8g -XX:-UsePerfData" \ + gatk --java-options "-Xmx${avail_mem}M -XX:-UsePerfData" \ HaplotypeCaller \ -R ${db} \ -I ${bam_file} \ diff --git a/modules/haplotypecaller_truth_joint.nf b/modules/haplotypecaller_truth_joint.nf index e99ab2b..15b0711 100644 --- a/modules/haplotypecaller_truth_joint.nf +++ b/modules/haplotypecaller_truth_joint.nf @@ -4,7 +4,7 @@ process HC_TRUTH_JOINT { // Call variants using GATK - HaplotypeCaller conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + container workflow.containerEngine == 'singularity' ? params.container.singularity.gatk : params.container.docker.gatk publishDir "${params.outputDir}/log/haplotypecaller/", pattern: "${sample_id}.g.haplotypecaller.log", mode:'copy' publishDir "${params.outputDir}/variants/", pattern: "${sample_id}.GATK.g.vcf.gz", mode:'copy' @@ -20,9 +20,10 @@ process HC_TRUTH_JOINT { path "${sample_id}.g.haplotypecaller.log" script: - def db = file(params.reference_genome).getName() + ".fna" + def db = file(params.reference_genome).name + def avail_mem = (task.memory.mega*0.8).intValue() """ - gatk --java-options "-Xmx8g -XX:-UsePerfData" \ + gatk --java-options "-Xmx${avail_mem}M -XX:-UsePerfData" \ HaplotypeCaller \ -R ${db} \ -L ${bedfile} \ diff --git a/modules/hs_metrics.nf b/modules/hs_metrics.nf index febd737..e3207f3 100644 --- a/modules/hs_metrics.nf +++ b/modules/hs_metrics.nf @@ -4,7 +4,7 @@ process HS_METRICS { // Collect HSMetrics for bam file conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + container workflow.containerEngine == 'singularity' ? params.container.singularity.gatk : params.container.docker.gatk publishDir "${params.outputDir}/log/hs_metrics/", pattern: "${sample_id}*.hs.txt", mode:'copy' @@ -18,7 +18,7 @@ process HS_METRICS { path("${sample_id}*.hs.txt"), emit: metrics_file script: - def db = file(params.reference_genome).getName() + ".fna" + def db = file(params.reference_genome).name def log_filename = log_suffix == "" ? "${sample_id}.hs.txt" : "${sample_id}_${log_suffix}.hs.txt" """ gatk BedToIntervalList \ diff --git a/modules/indelqual.nf b/modules/indelqual.nf index 5aca450..495e808 100644 --- a/modules/indelqual.nf +++ b/modules/indelqual.nf @@ -4,7 +4,7 @@ process INDELQUAL { // Tag indel quality using Lofreq - indelqual conda "$projectDir/envs/lofreq/environment.yaml" - container params.container.lofreq + container workflow.containerEngine == 'singularity' ? params.container.singularity.lofreq : params.container.docker.lofreq publishDir "${params.outputDir}/log/lofreq_indelqual/", pattern: "${sample_id}.lofreq.iq.log", mode:'copy' @@ -17,7 +17,7 @@ process INDELQUAL { path "${sample_id}.lofreq.iq.log" script: - def db = file(params.reference_genome).getName() + ".fna" + def db = file(params.reference_genome).name """ lofreq indelqual \ --dindel \ diff --git a/modules/index.nf b/modules/index.nf index 8407e1b..fe8cbd5 100644 --- a/modules/index.nf +++ b/modules/index.nf @@ -3,7 +3,7 @@ process INDEX { tag "Index bam file - $sample_id" conda "$projectDir/envs/samtools/environment.yaml" - container params.container.samtools + container workflow.containerEngine == 'singularity' ? params.container.singularity.samtools : params.container.docker.samtools input: tuple val(sample_id), path(bam_file) diff --git a/modules/index_vcf.nf b/modules/index_vcf.nf index f1c65e0..676f423 100644 --- a/modules/index_vcf.nf +++ b/modules/index_vcf.nf @@ -3,7 +3,7 @@ process INDEX_VCF { tag "Index VCF $sample_id" conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + container workflow.containerEngine == 'singularity' ? params.container.singularity.gatk : params.container.docker.gatk publishDir "${params.outputDir}/log/index_vcf/", pattern: "${sample_id}.${caller}.log", mode:'copy' diff --git a/modules/insert_size_metrics.nf b/modules/insert_size_metrics.nf index 6f34489..c7eed54 100644 --- a/modules/insert_size_metrics.nf +++ b/modules/insert_size_metrics.nf @@ -3,8 +3,8 @@ process INSERT_SIZE_METRICS { tag "$sample_id" // Collect insert size metrics for bam file - conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + conda "$projectDir/envs/picard/environment.yaml" + container workflow.containerEngine == 'singularity' ? params.container.singularity.picard : params.container.docker.picard publishDir "${params.outputDir}/log/insert_size_metrics/", pattern: "${sample_id}*.insert.txt", mode:'copy' publishDir "${params.outputDir}/log/insert_size_metrics/", pattern: "${sample_id}.pdf", mode:'copy' @@ -19,10 +19,10 @@ process INSERT_SIZE_METRICS { script: def log_filename = log_suffix == "" ? "${sample_id}.insert.txt" : "${sample_id}_${log_suffix}.insert.txt" """ - gatk CollectInsertSizeMetrics \ + picard CollectInsertSizeMetrics \ -I ${bam_file} \ -O ${log_filename} \ - -H ${sample_id}.pdf + -H plot.pdf """ stub: diff --git a/modules/intervals.nf b/modules/intervals.nf index 2b335f3..2d2f873 100644 --- a/modules/intervals.nf +++ b/modules/intervals.nf @@ -3,7 +3,7 @@ process INTERVALS { tag "BedFileToIntervalList" conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + container workflow.containerEngine == 'singularity' ? params.container.singularity.gatk : params.container.docker.gatk input: path reference_genome @@ -13,7 +13,7 @@ process INTERVALS { path "target_region.interval_list", emit: target_list script: - def db = file(params.reference_genome).getName() + ".fna" + def db = file(params.reference_genome).name """ gatk BedToIntervalList \ I=${bedfile} \ diff --git a/modules/lofreq.nf b/modules/lofreq.nf index 9272533..c600568 100644 --- a/modules/lofreq.nf +++ b/modules/lofreq.nf @@ -4,7 +4,7 @@ process LOFREQ { // Call variants using Lofreq - call parallel conda "$projectDir/envs/lofreq/environment.yaml" - container params.container.lofreq + container workflow.containerEngine == 'singularity' ? params.container.singularity.lofreq : params.container.docker.lofreq publishDir "${params.outputDir}/log/lofreq/", pattern: "${sample_id}.lofreq.log", mode:'copy' publishDir "${params.outputDir}/variants/", pattern: "${sample_id}.lofreq.vcf.gz", mode:'copy' @@ -20,7 +20,7 @@ process LOFREQ { path "${sample_id}.lofreq.log" script: - def db = file(params.reference_genome).getName() + ".fna" + def db = file(params.reference_genome).name def publishDir = file(params.outputDir + "/variants/" + sample_id + ".lofreq.vcf.gz") """ lofreq call-parallel \ diff --git a/modules/markduplicates.nf b/modules/markduplicates.nf index 0a165ac..0936086 100644 --- a/modules/markduplicates.nf +++ b/modules/markduplicates.nf @@ -4,7 +4,7 @@ process MARKDUPLICATES { // Mark duplicate reads in BAM files conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + container workflow.containerEngine == 'singularity' ? params.container.singularity.gatk : params.container.docker.gatk publishDir "${params.outputDir}/log/markduplicates/", pattern: "${sample_id}*.log", mode:'copy' @@ -20,8 +20,9 @@ process MARKDUPLICATES { script: def log_filename = log_suffix == "" ? "${sample_id}.dupMetric.log" : "${sample_id}_${log_suffix}.dupMetric.log" def optical_only = optical_only_tag ? "--duplicate-tagging-policy OpticalOnly" : "" + def avail_mem = (task.memory.mega*0.8).intValue() """ - gatk --java-options -Xmx16g MarkDuplicates \ + gatk --java-options -Xmx${avail_mem}M MarkDuplicates \ --OPTICAL_DUPLICATE_PIXEL_DISTANCE 2500 \ ${optical_only} \ --TMP_DIR . \ diff --git a/modules/markduplicates_fast.nf b/modules/markduplicates_fast.nf index 91b961d..1dea1e0 100644 --- a/modules/markduplicates_fast.nf +++ b/modules/markduplicates_fast.nf @@ -4,7 +4,7 @@ process MARKDUPLICATES_FAST { // Mark duplicate reads in BAM files conda "$projectDir/envs/sambamba/environment.yaml" - container params.container.sambamba + container workflow.containerEngine == 'singularity' ? params.container.singularity.sambamba : params.container.docker.sambamba publishDir "${params.outputDir}/log/markdup_sambamba/", pattern: "${sample_id}*.log", mode:'copy' diff --git a/modules/markduplicates_spark.nf b/modules/markduplicates_spark.nf index acebd17..26c075d 100644 --- a/modules/markduplicates_spark.nf +++ b/modules/markduplicates_spark.nf @@ -4,7 +4,7 @@ process MARKDUPLICATES_SPARK { // Mark duplicate reads in BAM files conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + container workflow.containerEngine == 'singularity' ? params.container.singularity.gatk : params.container.docker.gatk publishDir "${params.outputDir}/log/markduplicates/", pattern: "${sample_id}*.log", mode:'copy' @@ -20,8 +20,10 @@ process MARKDUPLICATES_SPARK { script: def log_filename = log_suffix == "" ? "${sample_id}.dupMetric.log" : "${sample_id}_${log_suffix}.dupMetric.log" def optical_only = optical_only_tag ? "--duplicate-tagging-policy OpticalOnly" : "" + def avail_mem = (task.memory.mega*0.8).intValue() """ - gatk --java-options -Xmx16g MarkDuplicatesSpark \ + gatk --java-options -Xmx${avail_mem}M \ + MarkDuplicatesSpark \ --optical-duplicate-pixel-distance 2500 \ ${optical_only} \ --tmp-dir . \ diff --git a/modules/matrix_context.nf b/modules/matrix_context.nf new file mode 100644 index 0000000..5d7ee64 --- /dev/null +++ b/modules/matrix_context.nf @@ -0,0 +1,35 @@ +process MATRIX_CONTEXT { + label 'process_low' + conda "$projectDir/envs/pinpy/environment.yaml" + container workflow.containerEngine == 'singularity' ? params.container.singularity.pinpy : params.container.docker.pinpy + + publishDir "${params.outputDir}/context/", mode:'copy', pattern: "matrix_context.*" + publishDir "${params.outputDir}/context/", mode:'copy', pattern: "decodetable.tsv" + + input: + path pooltable + path decodetable, stageAs: 'input_table/*' + + output: + path "matrix_context.tsv" + path "matrix_context.json", emit: json + path "decodetable.tsv", emit: decodetable + + script: + def decode = decodetable ? "--decode ${decodetable}" : "" + """ + build_matrix.py \ + --pools ${pooltable} \ + --out-json matrix_context.json \ + --out-tsv matrix_context.tsv \ + --out-decode decodetable.tsv \ + ${decode} + """ + + stub: + """ + touch matrix_context.json + touch matrix_context.tsv + touch decodetable.tsv + """ +} \ No newline at end of file diff --git a/modules/merge_consensus.nf b/modules/merge_consensus.nf index d7707e4..5d9d40a 100644 --- a/modules/merge_consensus.nf +++ b/modules/merge_consensus.nf @@ -3,7 +3,7 @@ process MERGE_CONSENSUS { tag "Merge consensus bam files - $sample_id" conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + container workflow.containerEngine == 'singularity' ? params.container.singularity.gatk : params.container.docker.gatk input: tuple val(sample_id), path(bam_file, stageAs: "raw/*"), path(ubam_file, stageAs: "raw/*") @@ -13,9 +13,10 @@ process MERGE_CONSENSUS { tuple val(sample_id), path("${sample_id}.bam"), emit: bam_file script: - def db = file(params.reference_genome).getName() + ".fna" + def db = file(params.reference_genome).name + def avail_mem = (task.memory.mega*0.8).intValue() """ - gatk --java-options -Xmx20g MergeBamAlignment \ + gatk --java-options -Xmx${avail_mem}M MergeBamAlignment \ TMP_DIR=. \ UNMAPPED=${ubam_file} \ ALIGNED=${bam_file} \ diff --git a/modules/mergebam.nf b/modules/mergebam.nf index e03b84f..d270b94 100644 --- a/modules/mergebam.nf +++ b/modules/mergebam.nf @@ -3,7 +3,7 @@ process MERGEBAM { tag "Merge bam files - $sample_id" conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + container workflow.containerEngine == 'singularity' ? params.container.singularity.gatk : params.container.docker.gatk input: tuple val(sample_id), path(bam_file, stageAs: "raw/*"), path(ubam_file, stageAs: "raw/*") @@ -13,7 +13,7 @@ process MERGEBAM { tuple val(sample_id), path("${sample_id}_raw.bam"), emit: bam_file script: - def db = file(params.reference_genome).getName() + ".fna" + def db = file(params.reference_genome).name """ gatk MergeBamAlignment \ TMP_DIR=. \ diff --git a/modules/mergepins.nf b/modules/mergepins.nf index 1900f33..09600f1 100644 --- a/modules/mergepins.nf +++ b/modules/mergepins.nf @@ -4,7 +4,7 @@ process MERGE_PINS { // Merge pinpointables into a single VCF conda "$projectDir/envs/bcftools/environment.yaml" - container params.container.bcftools + container workflow.containerEngine == 'singularity' ? params.container.singularity.bcftools : params.container.docker.bcftools publishDir "${params.outputDir}/", mode:'copy' diff --git a/modules/mergevcfs.nf b/modules/mergevcfs.nf index 3845fba..4cdceff 100644 --- a/modules/mergevcfs.nf +++ b/modules/mergevcfs.nf @@ -4,7 +4,7 @@ process MERGEVCFS { // Call variants using GATK - HaplotypeCaller conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + container workflow.containerEngine == 'singularity' ? params.container.singularity.gatk : params.container.docker.gatk publishDir "${params.outputDir}/log/mergevcfs/", pattern: "${sample_id}.mergevcfs.log", mode:'copy' publishDir "${params.outputDir}/variants/", pattern: "${sample_id}.GATK.vcf.gz", mode:'copy' diff --git a/modules/mosdepth.nf b/modules/mosdepth.nf index 7cb93fb..6ef7d80 100644 --- a/modules/mosdepth.nf +++ b/modules/mosdepth.nf @@ -6,7 +6,7 @@ process MOSDEPTH { // 5x allele depth (240), 10x allele depth (480), 30x allele depth (1440), 50x allele depth (2400) conda "$projectDir/envs/mosdepth/environment.yaml" - container params.container.mosdepth + container workflow.containerEngine == 'singularity' ? params.container.singularity.mosdepth : params.container.docker.mosdepth publishDir "${params.outputDir}/log/mosdepth/", pattern: "${sample_id}*.per-base.bed.gz", mode:'copy' publishDir "${params.outputDir}/log/mosdepth/", pattern: "${sample_id}*.regions.bed.gz", mode:'copy' diff --git a/modules/mpileup.nf b/modules/mpileup.nf new file mode 100644 index 0000000..c75e4b6 --- /dev/null +++ b/modules/mpileup.nf @@ -0,0 +1,36 @@ +process MPILEUP { + label 'process_low' + tag "BAM->mpileup - $sample_id" + + conda "$projectDir/envs/samtools/environment.yaml" + container workflow.containerEngine == 'singularity' ? params.container.singularity.samtools : params.container.docker.samtools + + publishDir "${params.outputDir}/mpileup/", mode:'copy', pattern: "${sample_id}.mpileup.gz" + + input: + tuple val(sample_id), path(bam_file), path(index) + path(reference) + path(bedfile) + + output: + path("${sample_id}.mpileup.gz"), emit: mpileup_file + + script: + def db = file(params.reference_genome).name + """ + samtools mpileup \ + -l ${bedfile} \ + --max-depth 50000 \ + --min-BQ 0 \ + --output-MQ \ + --fasta-ref ${db} \ + ${bam_file} \ + | gzip > ${sample_id}.mpileup.gz + """ + + stub: + """ + touch "${sample_id}.mpileup.gz" + """ +} + diff --git a/modules/multiqc.nf b/modules/multiqc.nf index beb6a61..42ce000 100644 --- a/modules/multiqc.nf +++ b/modules/multiqc.nf @@ -3,7 +3,7 @@ process MULTIQC { tag "MultiQC everything" conda "$projectDir/envs/multiqc/environment.yaml" - container params.container.multiqc + container workflow.containerEngine == 'singularity' ? params.container.singularity.multiqc : params.container.docker.multiqc publishDir "${params.outputDir}/log/multiqc/", mode:'copy' diff --git a/modules/ngsep.nf b/modules/ngsep.nf new file mode 100644 index 0000000..1a23405 --- /dev/null +++ b/modules/ngsep.nf @@ -0,0 +1,41 @@ +process NGSEP { + label 'process_low' + tag "NGSEP - $sample_id" + // Call variants using NGSEP + + conda "$projectDir/envs/ngsep/environment.yaml" + container workflow.containerEngine == 'singularity' ? params.container.singularity.ngsep : params.container.docker.ngsep + + publishDir "${params.outputDir}/log/ngsep/", pattern: "${sample_id}.ngsep.log", mode:'copy' + publishDir "${params.outputDir}/variants/", pattern: "${sample_id}.ngsep.vcf", mode:'copy' + + input: + tuple val(sample_id), path(bam_file), path(bam_index_file) + path reference_genome + path bedfile + + output: + tuple val(sample_id), path("${sample_id}.ngsep.vcf"), emit: vcf_file + path "${sample_id}.ngsep.log" + + script: + def db = file(params.reference_genome).name + """ + ngsep \ + SingleSampleVariantsDetector \ + -maxAlnsPerStartPos 50000 \ + -h 0.1 \ + -ploidy ${params.ploidy} \ + -r ${db} \ + -i ${bam_file} \ + -o "${sample_id}.ngsep" \ + > >(tee -a "${sample_id}.ngsep.log") \ + 2> >(tee -a "${sample_id}.ngsep.log" >&2) + """ + + stub: + """ + touch "${sample_id}.ngsep.vcf" "${sample_id}.ngsep.log" + """ +} + diff --git a/modules/normalise_vcf.nf b/modules/normalise_vcf.nf index 540649d..78d37c7 100644 --- a/modules/normalise_vcf.nf +++ b/modules/normalise_vcf.nf @@ -6,7 +6,7 @@ process NORMALISE_VCF { // 3. Checks refs against reference genome conda "$projectDir/envs/bcftools/environment.yaml" - container params.container.bcftools + container workflow.containerEngine == 'singularity' ? params.container.singularity.bcftools : params.container.docker.bcftools publishDir "${params.outputDir}/log/normalise_vcf/${sample_id}/", pattern: "${sample_id}.${caller}.log", mode:'copy' publishDir "${params.outputDir}/normalised_variants/", pattern: "${sample_id}.${caller}.vcf", mode:'copy' @@ -21,7 +21,7 @@ process NORMALISE_VCF { path "${sample_id}.${caller}.log" script: - def db = file(params.reference_genome).getName() + ".fna" + def db = file(params.reference_genome).name """ bcftools norm \ --check-ref e \ diff --git a/modules/octopus.nf b/modules/octopus.nf index 00aad12..3cd4b4e 100644 --- a/modules/octopus.nf +++ b/modules/octopus.nf @@ -8,7 +8,7 @@ process OCTOPUS { // This version works on NGC only. To run locally, comment out "singularity exec" and use "octopus" directly. conda "$projectDir/envs/octopus/environment.yaml" - container params.container.octopus + container workflow.containerEngine == 'singularity' ? params.container.singularity.octopus : params.container.docker.octopus publishDir "${params.outputDir}/log/octopus/", pattern: "${sample_id}.octopus.log", mode:'copy' publishDir "${params.outputDir}/variants/", pattern: "${sample_id}.octopus.vcf.gz", mode:'copy' diff --git a/modules/pilot_pinpoint.nf b/modules/pilot_pinpoint.nf index bcdb9ad..a4030fe 100644 --- a/modules/pilot_pinpoint.nf +++ b/modules/pilot_pinpoint.nf @@ -1,7 +1,7 @@ process PILOT_PINPOINT { label 'process_low' conda "$projectDir/envs/r_env/environment.yaml" - container params.container.r_env + container workflow.containerEngine == 'singularity' ? params.container.singularity.r_env : params.container.docker.r_env publishDir "${params.outputDir}/pinned_variants/", pattern: "*.tsv", mode:'copy' publishDir "${params.outputDir}/pinned_variants/outlier_plots/", pattern: "*.pdf", mode: 'copy' diff --git a/modules/pin_basic.nf b/modules/pin_basic.nf new file mode 100644 index 0000000..6ccf3e5 --- /dev/null +++ b/modules/pin_basic.nf @@ -0,0 +1,29 @@ +process PIN_BASIC { + label 'process_low' + conda "$projectDir/envs/filter_variants/environment.yaml" + container workflow.containerEngine == 'singularity' ? params.container.singularity.filter_variants : params.container.docker.filter_variants + + publishDir "${params.outputDir}/variant_compilation/", mode:'copy', pattern: "*variants.tsv" + + input: + path vcf_files + path matrix_context + + output: + path "all_pool_variants.tsv", emit: pool_variants + path "pinpoint_variants.tsv", emit: pinpoint_variants + + script: + """ + pin_basic.py \ + --context ${matrix_context} \ + --vcf-folder . \ + --output . + """ + + stub: + """ + touch all_pool_variants.tsv + touch pinpoint_variants.tsv + """ +} \ No newline at end of file diff --git a/modules/pinpy.nf b/modules/pinpy.nf index 0fdc6b1..ac0b261 100644 --- a/modules/pinpy.nf +++ b/modules/pinpy.nf @@ -1,7 +1,7 @@ process PINPY { label 'process_low' conda "$projectDir/envs/pinpy/environment.yaml" - container params.container.pinpy + container workflow.containerEngine == 'singularity' ? params.container.singularity.pinpy : params.container.docker.pinpy publishDir "${params.outputDir}/", mode:'copy' diff --git a/modules/report.nf b/modules/report.nf new file mode 100644 index 0000000..6a212b4 --- /dev/null +++ b/modules/report.nf @@ -0,0 +1,42 @@ +process REPORT { + label 'process_low' + conda "$projectDir/envs/r_env/environment.yaml" + container workflow.containerEngine == 'singularity' ? params.container.singularity.r_env : params.container.docker.r_env + + publishDir "${params.outputDir}/", mode:'copy', pattern: "variant_report.html" + + input: + path rmd + path matrix_context + path all_variants + path pinpoints + path annotations + path rescue + + output: + path "variant_report.html" + + script: + def annotate = annotations ? "annotations='${annotations}'," : "annotations=NULL," + def rescue_probs = rescue ? "rescue='${rescue}'" : "rescue=NULL" + """ + Rscript -e \ + "rmarkdown::render( + '${rmd}', + params = list( + matrix_context = '${matrix_context}', + pinpoints = '${pinpoints}', + all_variants = '${all_variants}', + ${annotate} + ${rescue_probs} + ), + output_file = 'variant_report.html', + quiet = TRUE + )" + """ + + stub: + """ + touch variant_report.html + """ +} \ No newline at end of file diff --git a/modules/rescue.nf b/modules/rescue.nf new file mode 100644 index 0000000..fe9b3d5 --- /dev/null +++ b/modules/rescue.nf @@ -0,0 +1,33 @@ +process RESCUE { + label 'process_low' + conda "$projectDir/envs/rescue/environment.yaml" + container workflow.containerEngine == 'singularity' ? params.container.singularity.marbl : params.container.docker.marbl + + + publishDir "${params.outputDir}/rescue_probabilities/", mode:'copy', pattern: "predictions.tsv" + + input: + path sampletable + path decodetable + path vcf_files + path mpileup + + output: + path "predictions.tsv", emit: probabilities + + script: + def decode = decodetable ? "--decodetable ${decodetable}" : "" + """ + marbl \ + --vcf-folder . \ + --mpileup-folder . \ + --sampletable ${sampletable} \ + ${decode} \ + --output . + """ + + stub: + """ + touch predictions.tsv + """ +} diff --git a/modules/snpeff.nf b/modules/snpeff.nf index 09745ae..8b8e90b 100644 --- a/modules/snpeff.nf +++ b/modules/snpeff.nf @@ -8,7 +8,7 @@ process SNPEFF { // -noMotif avoids adding motif annotations. conda "$projectDir/envs/snpeff/environment.yaml" - container params.container.snpeff + container workflow.containerEngine == 'singularity' ? params.container.singularity.snpeff : params.container.docker.snpeff publishDir "${params.outputDir}/log/snpeff/${sample_id}/", pattern: "${sample_id}.${caller}.log", mode:'copy' publishDir "${params.outputDir}/log/snpeff/stats/", pattern: "${sample_id}.${caller}_snpeff_stats.csv", mode:'copy' @@ -27,7 +27,7 @@ process SNPEFF { script: """ - ${params.snpeff} \ + snpEff \ ${snpeff_db} \ -c ${snpeff_config} \ -csvStats "${sample_id}.${caller}_snpeff_stats.csv" \ diff --git a/modules/snpsift_clinvar.nf b/modules/snpsift_clinvar.nf index e50c87a..af13739 100644 --- a/modules/snpsift_clinvar.nf +++ b/modules/snpsift_clinvar.nf @@ -4,8 +4,8 @@ process SNPSIFT_CLINVAR { // Add ClinVar annotations to a VCF file using SnpSift. - conda "$projectDir/envs/snpeff/environment.yaml" - container params.container.snpeff + conda "$projectDir/envs/snpsift/environment.yaml" + container workflow.containerEngine == 'singularity' ? params.container.singularity.snpsift : params.container.docker.snpsift publishDir "${params.outputDir}/log/snpsift_clinvar/${sample_id}/", pattern: "${sample_id}.${caller}.log", mode:'copy' publishDir "${params.outputDir}/annotated_variants/", pattern: "${sample_id}.${caller}.vcf", mode:'copy' @@ -22,7 +22,7 @@ process SNPSIFT_CLINVAR { script: def db = file(params.clinvar_db).getName() """ - ${params.snpsift} annotate \ + SnpSift annotate \ ${db} \ ${vcf_file} \ > ${sample_id}.${caller}.vcf \ diff --git a/modules/snpsift_filter.nf b/modules/snpsift_filter.nf index 35720fc..ea2dd31 100644 --- a/modules/snpsift_filter.nf +++ b/modules/snpsift_filter.nf @@ -4,8 +4,8 @@ process SNPSIFT_FILTER { // Filter variants based on ClinVar and SNPEff annotations using SnpSift. - conda "$projectDir/envs/snpeff/environment.yaml" - container params.container.snpeff + conda "$projectDir/envs/snpsift/environment.yaml" + container workflow.containerEngine == 'singularity' ? params.container.singularity.snpsift : params.container.docker.snpsift publishDir "${params.outputDir}/annotated_variants/lof/", pattern: "${sample_id}.${caller}.lof.vcf", mode:'copy' publishDir "${params.outputDir}/annotated_variants/pathogenic/", pattern: "${sample_id}.${caller}.p.vcf", mode:'copy' @@ -22,17 +22,17 @@ process SNPSIFT_FILTER { script: """ - ${params.snpsift} filter \ + SnpSift filter \ "(exists LOF)" \ ${vcf_file} \ > ${sample_id}.${caller}.lof.vcf - ${params.snpsift} filter \ + SnpSift filter \ "CLNSIG =~ 'Pathogenic'" \ ${vcf_file} \ > ${sample_id}.${caller}.p.vcf - ${params.snpsift} filter \ + SnpSift filter \ "((exists LOF) | (CLNSIG =~ 'Pathogenic'))" \ ${vcf_file} \ > ${sample_id}.${caller}.lofp.vcf diff --git a/modules/subset.nf b/modules/subset.nf index e960ab7..4877106 100644 --- a/modules/subset.nf +++ b/modules/subset.nf @@ -3,7 +3,7 @@ process SUBSET { tag "Subset alignment and convert to bam - $sample_id" conda "$projectDir/envs/samtools/environment.yaml" - container params.container.samtools + container workflow.containerEngine == 'singularity' ? params.container.singularity.samtools : params.container.docker.samtools input: tuple val(sample_id), path(bam_file, stageAs: 'raw/*') @@ -14,7 +14,7 @@ process SUBSET { tuple val(sample_id), path("${sample_id}.bam"), emit: bam_file script: - def db = file(params.reference_genome).getName() + ".fna" + def db = file(params.reference_genome).name """ samtools view \ -@ ${task.cpus} \ diff --git a/modules/test.nf b/modules/test.nf index 357b284..eaa4736 100644 --- a/modules/test.nf +++ b/modules/test.nf @@ -1,7 +1,7 @@ process TEST { label 'process_single' conda "$projectDir/envs/pinpy/environment.yaml" - container params.container.pinpy + container workflow.containerEngine == 'singularity' ? params.container.singularity.pinpy : params.container.docker.pinpy publishDir "${params.outputDir}/test/", pattern: "test.*", mode:'copy' @@ -11,8 +11,8 @@ process TEST { val pinpoint_method output: - path "test.passed" optional true - path "test.failed" optional true + path "test.passed", optional: true + path "test.failed", optional: true path "test.log" script: @@ -29,4 +29,4 @@ process TEST { touch test.passed touch test.log """ -} \ No newline at end of file +} diff --git a/modules/ubam.nf b/modules/ubam.nf index fb6cd90..6704ea3 100644 --- a/modules/ubam.nf +++ b/modules/ubam.nf @@ -4,7 +4,7 @@ process UBAM { // Convert FastQ files to unaligned BAM files. conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + container workflow.containerEngine == 'singularity' ? params.container.singularity.gatk : params.container.docker.gatk input: tuple val(sample_id), path(reads) diff --git a/modules/umi_metrics.nf b/modules/umi_metrics.nf index 1200b1b..91d9494 100644 --- a/modules/umi_metrics.nf +++ b/modules/umi_metrics.nf @@ -17,7 +17,7 @@ process UMI_METRICS { // DS: double-stranded tag families where membership is collapsed across single-stranded tag families from the same double-stranded source molecule; i.e. 50/A and 50/B become one family conda "$projectDir/envs/fgbio/environment.yaml" - container params.container.fgbio + container workflow.containerEngine == 'singularity' ? params.container.singularity.fgbio : params.container.docker.fgbio publishDir "${params.outputDir}/log/umi_metrics/logs/", pattern: "${sample_id}.*.txt", mode:'copy' publishDir "${params.outputDir}/log/umi_metrics/plots/", pattern: "${sample_id}.*.pdf", mode:'copy' @@ -34,7 +34,7 @@ process UMI_METRICS { script: """ - ${params.fgbio} --tmp-dir=. --compression 1 --async-io CollectDuplexSeqMetrics \ + fgbio --tmp-dir=. --compression 1 --async-io CollectDuplexSeqMetrics \ --input=${bam_file} \ --output="${sample_id}" """ diff --git a/modules/unique_vcf.nf b/modules/unique_vcf.nf new file mode 100644 index 0000000..8af2bcb --- /dev/null +++ b/modules/unique_vcf.nf @@ -0,0 +1,31 @@ +process UNIQUE_VCF { + label 'process_single' + // Merge, deduplicate, remove sample information, and index a set of VCF files. + + conda "$projectDir/envs/bcftools/environment.yaml" + container workflow.containerEngine == 'singularity' ? params.container.singularity.bcftools : params.container.docker.bcftools + + publishDir "${params.outputDir}/vep_annotate/", pattern: "unique.vcf.*", mode:'copy' + + input: + path vcf_files + + output: + tuple path("unique.vcf.gz"), path("unique.vcf.gz.tbi"), path('variant_keys.txt'), emit: vcf_file + + script: + """ + bcftools merge -m none --force-samples *.vcf.gz \ + | bcftools norm -d exact -O u \ + | bcftools view -G -Oz -o unique.vcf.gz + bcftools index --tbi unique.vcf.gz + bcftools query -f '%CHROM:%POS:%REF:%ALT\n' unique.vcf.gz > variant_keys.txt + """ + + stub: + """ + touch "unique.vcf.gz" + touch "unique.vcf.gz.tbi" + touch "variant_keys.txt" + """ +} diff --git a/modules/validate.nf b/modules/validate.nf index 2a55e29..7084d4e 100644 --- a/modules/validate.nf +++ b/modules/validate.nf @@ -3,7 +3,7 @@ process VALIDATE { tag "Validate bam file - $sample_id" conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + container workflow.containerEngine == 'singularity' ? params.container.singularity.gatk : params.container.docker.gatk publishDir "${params.outputDir}/log/validation/", pattern: "${sample_id}.validation.log", mode:'copy' diff --git a/modules/vartable.nf b/modules/vartable.nf index 1ed860b..203a590 100644 --- a/modules/vartable.nf +++ b/modules/vartable.nf @@ -3,7 +3,7 @@ process VARTABLE { tag "VCF to TSV $sample_id" conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + container workflow.containerEngine == 'singularity' ? params.container.singularity.gatk : params.container.docker.gatk publishDir "${params.outputDir}/log/variant_tables/${sample_id}/", pattern: "${sample_id}.${caller}.log", mode:'copy' publishDir "${params.outputDir}/variant_tables/", pattern: "${sample_id}.${caller}.tsv", mode:'copy' diff --git a/modules/vartable_pins.nf b/modules/vartable_pins.nf index c64c755..2c1509a 100644 --- a/modules/vartable_pins.nf +++ b/modules/vartable_pins.nf @@ -1,7 +1,7 @@ process VARTABLE_PINS { label 'process_single' conda "$projectDir/envs/gatk4/environment.yaml" - container params.container.gatk + container workflow.containerEngine == 'singularity' ? params.container.singularity.gatk : params.container.docker.gatk publishDir "${params.outputDir}/pinpoint_variant_tables/", mode:'copy' diff --git a/modules/vep.nf b/modules/vep.nf new file mode 100644 index 0000000..10e11bf --- /dev/null +++ b/modules/vep.nf @@ -0,0 +1,80 @@ +process VEP { + label 'process_single' + // Annotate using VEP + + conda "$projectDir/envs/vep/environment.yaml" + container workflow.containerEngine == 'singularity' ? params.container.singularity.vep : params.container.docker.vep + + publishDir "${params.outputDir}/vep_annotate/", pattern: "annotation*.tsv", mode:'copy' + + input: + tuple path(vcf_file), path(vcf_index), path(variant_keys) + path reference_genome + path cache_dir // optional cache directory + path utr_annotation // optional UTR annotation file + tuple path(alphamissense), path(am_index) // optional AlphaMissense file + tuple path(clinvar), path(clinvar_index) // optional ClinVar VCF + tuple path(danmac), path(danmar_index) // optional DanMac VCF + tuple path(blacklist), path(blacklist_index) // optional blacklist BED + tuple path(repeatmasks), path(repeat_index) // optional repeat masks BED + tuple path(gnomad), path(gnomad_index) // optional gnomAD VCF + path loftee_plugin_path, stageAs: "loftee" // optional LOFTEE github repo path + path loftee_gerp_bigwig // optional LOFTEE GERP BigWig + path loftee_human_ancestor // optional LOFTEE human ancestor FASTA + path loftee_conservation // optional LOFTEE conservation SQL + + output: + path "annotations.tsv" + path "annotations_w_varid.tsv", emit: annotations_file + + script: + def db = file(params.reference_genome).name + def mode = cache_dir ? "--cache --offline --dir_cache ${cache_dir}" : "--database" + def utr = utr_annotation ? "--plugin UTRAnnotator,file=${utr_annotation}" : "" + def am = alphamissense ? "--plugin AlphaMissense,file=${alphamissense}" : "" + def clinv = clinvar ? "--custom ${clinvar},ClinVar,vcf,exact,0,CLNSIG,CLNREVSTAT,CLNDN" : "" + def danm = danmac ? "--custom ${danmac},DanMac,vcf,exact,0,AC,AN" : "" + def blackl = blacklist ? "--custom ${blacklist},blacklist,bed,overlap,0,4" : "" + def repeatm = repeatmasks ? "--custom ${repeatmasks},repeats,bed,overlap,0,4" : "" + def gnom = gnomad ? "--custom ${gnomad},gnomAD,vcf,exact,0,AC_joint_nfe,AN_joint_nfe,AF_joint_nfe,AC_genomes_nfe,AN_genomes_nfe,AF_genomes_nfe,AC_exomes_nfe,AN_exomes_nfe,AF_exomes_nfe,grpmax_joint,AC_grpmax_joint,AF_grpmax_joint,grpmax_genomes,AC_grpmax_genomes,AF_grpmax_genomes,grpmax_exomes,AC_grpmax_exomes,AF_grpmax_exomes" : "" + def loftee_plugin_dir = workflow.containerEngine ? "/plugins" : "./loftee" + def loft = loftee_gerp_bigwig && loftee_human_ancestor && loftee_conservation ? "--plugin LoF,loftee_path:${loftee_plugin_dir},gerp_bigwig:${loftee_gerp_bigwig},human_ancestor_fa:${loftee_human_ancestor},conservation_file:${loftee_conservation}" : "" + """ + vep \ + --assembly GRCh38 \ + --format vcf \ + --tab \ + --input_file "${vcf_file}" \ + --output_file "annotations.tsv" \ + --fasta "${db}" \ + --check_ref --dont_skip --show_ref_allele \ + --hgvs --canonical --minimal --mane --biotype --numbers --allele_number --symbol \ + --pick --pick_order mane_select,mane_plus_clinical,canonical,appris,tsl,biotype,ccds,rank,length \ + --force_overwrite \ + --plugin SpliceRegion,extended \ + ${mode} \ + ${utr} \ + ${am} \ + ${clinv} \ + ${danm} \ + ${blackl} \ + ${repeatm} \ + ${gnom} \ + ${loft} + + # Header + printf 'varid\t' > annotations_w_varid.tsv + grep '^#[^#]' "annotations.tsv" >> annotations_w_varid.tsv + # Body + paste \ + <(grep -F -v '*' "${variant_keys}") \ + <(grep -v '^#' "annotations.tsv") \ + >> annotations_w_varid.tsv + """ + + stub: + """ + touch "annotations.tsv" + touch "annotations_w_varid.tsv" + """ +} diff --git a/next.pbs b/next.pbs index 78e78b8..a6cfa3b 100644 --- a/next.pbs +++ b/next.pbs @@ -14,12 +14,14 @@ # myqsub [-F "optional arguments to Nexflow"] next.pbs # -#PBS -l nodes=1:ppn=40:thinnode +#PBS -l nodes=1:ppn=40 #PBS -l walltime=24:00:00 #PBS -l mem=160gb #PBS -j oe #PBS -N nextflow +module load nextflow/25.04.6 GitPython + ## exit if any errors or unset variables are encountered set -euo pipefail @@ -29,7 +31,7 @@ version=2.0 root=/ngc/tools/ngctools/${name}/${version} data=/ngc/shared/${name} -export NXF_SINGULARITY_CACHEDIR=${data}/containers +export NXF_SINGULARITY_CACHEDIR=/ngc/projects2/dp_00005/data/dwf/bin/container_cache export NXF_ASSETS=${root}/pipelines export NXF_SCM_FILE=${root}/scm export NXF_PLUGINS_DIR=${data}/plugins @@ -37,8 +39,6 @@ export NXF_PLUGINS_DIR=${data}/plugins # Disable interactive logging #export NXF_ANSI_LOG=false -module load tools ngctools jdk/21.0.1 nextflow/24.08.0-edge - nextflow_essential_params='-profile ngc -resume' nextflow_command="nextflow run main.nf $nextflow_essential_params $@" diff --git a/next.ssi b/next.ssi new file mode 100644 index 0000000..b9bb59d --- /dev/null +++ b/next.ssi @@ -0,0 +1,51 @@ +#!/bin/bash +# +# next.ssi +# +# mads +# 2024-10-20 +# +# This is a simple wrapper that runs the DoBSeq workflow on SSI. +# The pipeline profile uses environment modules and runs all processes as submitted jobs. +# +# +# Usage: +# +# myqsub [-F "optional arguments to Nextflow"] next.pbs +# + +## exit if any errors or unset variables are encountered +set -euo pipefail + +# Disable interactive logging +#export NXF_ANSI_LOG=false + +nextflow_essential_params='-profile ssi -resume' + +nextflow_command="nextflow run main.nf $nextflow_essential_params $@" + +echo "-------- Log generated by next.ssi -------" | tee -a nextflow.log +echo "Command: $0 $@" | tee -a nextflow.log +echo "User:" `whoami` | tee -a nextflow.log +echo "Date:" `date -Is` | tee -a nextflow.log +echo "PBS job ID: ${PBS_JOBID-none}" | tee -a nextflow.log +echo 'DoBSeq version:' `cat VERSION` | tee -a nextflow.log +echo 'git commit:' `(git rev-parse HEAD)` | tee -a nextflow.log +echo 'git status:' | tee -a nextflow.log +(git status -bs) | sed "s/^/ /" | tee -a nextflow.log +echo 'Nextflow version:' `nextflow -version` 2>&1 | tee -a nextflow.log +echo "Nextflow command: ${nextflow_command}" | tee -a nextflow.log + +echo `date -Is` "--------------------- Starting Nextflow ---------------------" | tee -a nextflow.log + +set +eo pipefail + +${nextflow_command} 2>&1 | tee -a nextflow.log + +echo `date -Is` "--------------------- Nextflow finished ---------------------" | tee -a nextflow.log + +if [ -z ${PBS_JOBID+x} ]; then + echo "This was not run as a PBS job" | tee -a nextflow.log +else + echo "This was run as PBS job ${PBS_JOBID}" | tee -a nextflow.log +fi diff --git a/nextflow.config b/nextflow.config index ff24e4a..7b42900 100644 --- a/nextflow.config +++ b/nextflow.config @@ -4,9 +4,6 @@ params { - // ID for VCF output files: - outputId = "" - // Table containing pool ids and corresponding FastQ files pooltable = "" @@ -32,6 +29,9 @@ params { // Experiment ploidy ploidy = 48 + // HaplotypeCaller maximum number of haplotypes to evaluate: + max_haplotypes = 1000 + // Run HaplotypeCaller by intervals. Default is by chromosome. runHCParallel = true intervalList = (1..22).collect { "chr${it}" } + ["chrX"] @@ -69,11 +69,12 @@ params { // Variant calling methods lofreq = false - gatk_joint_calling = false + gatk_joint_calling = false // WGS only + deepvariant = false // WGS only (container only) crisp = false // Pool only (not available through conda) - octopus = false // Pool only (linux only) + octopus = false // Pool only (container only) freebayes = false // Pool only - deepvariant = false // WGS only (linux only) + ngsep = false // Pool only // Variant filtering. Either optimised for sensitivity or balance filter = false @@ -84,14 +85,38 @@ params { // Pinpoint method (pilot or new) pinpoint_method = 'new' + // Use new tabular output method (true or false) + pinpoint_tab = true + + // Use variant rescue model (true or false) + variant_rescue = false + + // Make pinpoint report (true or false) + pinpoint_report = false + // Annotation configurations - annotate = true + annotate = false snpeff_db = 'GRCh38.99' snpeff_config = "" snpeff_cache = "" clinvar_db = "" // Add gene annotation based on third column in bedfile - bedfile_gene_annotation = true + bedfile_gene_annotation = false + + // VEP annotations + annotate_vep = false + vep_cache = "" + clinvar_db = "" + danmac_db = "" + blacklist_bed = "" + repeatmasker_bed = "" + gnomad_vcf = "" + utr_file = "" + alphamissense_tsv = "" + loftee_path = "" + loftee_gerp_bw = "" + loftee_human_ancestor= "" + loftee_sqlite = "" // Misc (never edit here) help = false @@ -107,6 +132,7 @@ process.shell = ['/bin/bash', '-euo', 'pipefail'] includeConfig 'conf/profiles.config' includeConfig 'conf/container.config' includeConfig 'conf/ngc.config' +includeConfig 'conf/ssi.config' // Add additional log and tracing. Trace can be adjusted to specific outputs. diff --git a/subworkflows/annotation.nf b/subworkflows/annotation.nf index 9a7e9e8..862b4d7 100644 --- a/subworkflows/annotation.nf +++ b/subworkflows/annotation.nf @@ -13,9 +13,11 @@ if (params.bedfile_gene_annotation) { bedfile = Channel.fromPath(params.bedfile, checkIfExists: true).collect() } -snpeff_cache_ch = Channel.fromPath(params.snpeff_cache, checkIfExists: true).collect() -clinvardb_ch = Channel.fromPath(params.clinvar_db + "*", checkIfExists: true).collect() -snpeff_config_ch = Channel.fromPath(params.snpeff_config, checkIfExists: true).collect() +if (params.annotate) { + snpeff_cache_ch = Channel.fromPath(params.snpeff_cache, checkIfExists: true).collect() + clinvardb_ch = Channel.fromPath(params.clinvar_db + "*", checkIfExists: true).collect() + snpeff_config_ch = Channel.fromPath(params.snpeff_config, checkIfExists: true).collect() +} workflow ANNOTATION { take: diff --git a/subworkflows/calling.nf b/subworkflows/calling.nf index da4456f..72ded77 100644 --- a/subworkflows/calling.nf +++ b/subworkflows/calling.nf @@ -12,6 +12,7 @@ include { DEEPVARIANT } from '../modules/deepvariant' include { FREEBAYES } from '../modules/freebayes' include { CRISP } from '../modules/crisp' include { OCTOPUS } from '../modules/octopus' +include { NGSEP } from '../modules/ngsep' include { FILTER } from '../modules/filter' include { GENOMICSDB } from '../modules/genomicsdb' include { GENOTYPEGVCF } from '../modules/genotypegvcf' @@ -57,6 +58,10 @@ workflow CALLING { FREEBAYES(bam_file_w_index, reference_genome, bedfile) } + if (params.ngsep) { + NGSEP(bam_file_w_index, reference_genome, bedfile) + } + if (params.runHCParallel) { // Split by interval list to multi-thread calling Channel diff --git a/subworkflows/mapping.nf b/subworkflows/mapping.nf index 41ee236..c6c482a 100644 --- a/subworkflows/mapping.nf +++ b/subworkflows/mapping.nf @@ -27,7 +27,6 @@ include { HS_METRICS } from '../modules/hs_metrics' include { ALIGNMENT_METRICS } from '../modules/alignment_metrics' include { GC_METRICS } from '../modules/gc_metrics' include { INSERT_SIZE_METRICS } from '../modules/insert_size_metrics' -include { DUPLICATE_METRICS } from '../modules/duplicate_metrics' include { MULTIQC } from '../modules/multiqc' @@ -75,7 +74,6 @@ workflow MAPPING { ALIGNMENT_METRICS(ALIGNMENT.out.raw_bam_file, reference_genome, "") GC_METRICS(ALIGNMENT.out.raw_bam_file, reference_genome, "") INSERT_SIZE_METRICS(ALIGNMENT.out.raw_bam_file, "") - DUPLICATE_METRICS(ALIGNMENT.out.raw_bam_file, "") qc_ch = qc_ch.mix( RAW_DEPTH.out.region_dist, diff --git a/subworkflows/mapping_umi.nf b/subworkflows/mapping_umi.nf index c284f0d..551c989 100644 --- a/subworkflows/mapping_umi.nf +++ b/subworkflows/mapping_umi.nf @@ -31,7 +31,6 @@ include { HS_METRICS } from '../modules/hs_metrics' include { ALIGNMENT_METRICS } from '../modules/alignment_metrics' include { GC_METRICS } from '../modules/gc_metrics' include { INSERT_SIZE_METRICS } from '../modules/insert_size_metrics' -include { DUPLICATE_METRICS } from '../modules/duplicate_metrics' include { INDEX as RAW_INDEX } from '../modules/index' include { FLAGSTAT as RAW_FLAGSTAT } from '../modules/flagstat' include { MOSDEPTH as RAW_DEPTH } from '../modules/mosdepth' @@ -137,7 +136,6 @@ workflow MAPPING_UMI { ALIGNMENT_METRICS(MERGEBAM.out.bam_file, reference_genome, "") GC_METRICS(MERGEBAM.out.bam_file, reference_genome, "") INSERT_SIZE_METRICS(MERGEBAM.out.bam_file, "") - DUPLICATE_METRICS(MERGEBAM.out.bam_file, "") UMI_METRICS(GROUP_UMI.out.grouped_bam_file) // After consensus calling (HS metrics) CONS_METRIC(ADDREADGROUP.out.bam_file, reference_genome, bedfile, "deduplicated") @@ -148,7 +146,6 @@ workflow MAPPING_UMI { GC_METRICS.out.metrics_file, GC_METRICS.out.summary_file, INSERT_SIZE_METRICS.out.metrics_file, - DUPLICATE_METRICS.out.metrics_file, CONS_METRIC.out.metrics_file) } diff --git a/subworkflows/pinpoint_tab.nf b/subworkflows/pinpoint_tab.nf new file mode 100644 index 0000000..9d4de33 --- /dev/null +++ b/subworkflows/pinpoint_tab.nf @@ -0,0 +1,67 @@ +/* +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + PINPOINTING WITH TABULAR OUTPUT SUB-WORKFLOW +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +*/ + +import static Utils.indexedFileChannel + +include { INDEX_VCF } from '../modules/index_vcf' +include { DROP_PL } from '../modules/drop_pl' +include { VARTABLE } from '../modules/vartable' +include { NORMALISE_VCF } from '../modules/normalise_vcf' +include { BGZIP } from '../modules/bgzip' +include { PIN_BASIC } from '../modules/pin_basic' +include { UNIQUE_VCF } from '../modules/unique_vcf' +include { VEP } from '../modules/vep' + +alphamissense_ch = indexedFileChannel(params.alphamissense_tsv, '.tbi') +clinvar_ch = indexedFileChannel(params.clinvar_db, '.tbi') +danmac_db_ch = indexedFileChannel(params.danmac_db, '.tbi') +blacklist_bed_ch = indexedFileChannel(params.blacklist_bed, '.tbi') +repeatmasker_bed_ch = indexedFileChannel(params.repeatmasker_bed, '.tbi') +gnomad_vcf_ch = indexedFileChannel(params.gnomad_vcf, '.tbi') + +annotations_ch = Channel.value([]) + +workflow PINPOINT_TAB { + take: + vcf_file + matrix_context + reference_genome + + main: + INDEX_VCF(vcf_file, 'GATK') + DROP_PL(INDEX_VCF.out.vcf_w_index, 'GATK') + NORMALISE_VCF(DROP_PL.out.vcf_file, reference_genome, 'GATK') + BGZIP(NORMALISE_VCF.out.norm_vcf) + BGZIP.out.bgzf_file + .map { sample, vcf, index -> tuple(vcf, index) } + .collect() + .set { basic_input } + PIN_BASIC(basic_input,matrix_context) + UNIQUE_VCF(basic_input) + if (params.annotate_vep) { + VEP( + UNIQUE_VCF.out.vcf_file, + reference_genome, + params.vep_cache ?: [], + params.utr_file ?: [], + alphamissense_ch ?: [], + clinvar_ch, + danmac_db_ch, + blacklist_bed_ch, + repeatmasker_bed_ch, + gnomad_vcf_ch, + params.loftee_path ?: [], + params.loftee_gerp_bw ?: [], + params.loftee_human_ancestor ?: [], + params.loftee_sqlite ?: [] + ) + annotations_ch = VEP.out.annotations_file + } + emit: + pinpoint_variants = PIN_BASIC.out.pinpoint_variants + pool_variants = PIN_BASIC.out.pool_variants + annotations = annotations_ch +} \ No newline at end of file diff --git a/subworkflows/pinpoint.nf b/subworkflows/pinpoint_vcf.nf similarity index 98% rename from subworkflows/pinpoint.nf rename to subworkflows/pinpoint_vcf.nf index d97da61..3e1f74b 100644 --- a/subworkflows/pinpoint.nf +++ b/subworkflows/pinpoint_vcf.nf @@ -15,17 +15,18 @@ include { ANNOTATION } from '../subworkflows/annotation' include { VARTABLE_PINS } from '../modules/vartable_pins' include { MERGE_PINS } from '../modules/mergepins' - if (params.filter) { snv_model = Channel.fromPath("$projectDir/assets/filter/logistic_regression_snv_model.joblib", checkIfExists: true).collect() indel_model = Channel.fromPath("$projectDir/assets/filter/random_forest_indel_model.joblib", checkIfExists: true).collect() model_class = Channel.fromPath("$projectDir/assets/filter/model.py", checkIfExists: true).collect() } -workflow PINPOINT { +workflow PINPOINT_VCF { take: vcf_file + pooltable decode_table + matrix_context reference_genome main: @@ -82,7 +83,6 @@ workflow PINPOINT { vartables, decode_table, 'GATK') - VARTABLE_PINS(PINPY.out.vcf_unique_pins.flatten()) MERGE_PINS(PINPY.out.vcf_unique_2d_pins) diff --git a/subworkflows/variant_rescue.nf b/subworkflows/variant_rescue.nf new file mode 100644 index 0000000..b72c1b1 --- /dev/null +++ b/subworkflows/variant_rescue.nf @@ -0,0 +1,38 @@ +/* +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + VARIANT_RESCUE SUB-WORKFLOW +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ +*/ + +include { INDEX_VCF } from '../modules/index_vcf' +include { DROP_PL } from '../modules/drop_pl' +include { NORMALISE_VCF } from '../modules/normalise_vcf' +include { BGZIP } from '../modules/bgzip' +include { MPILEUP } from '../modules/mpileup' +include { RESCUE } from '../modules/rescue' + +workflow VARIANT_RESCUE { + take: + vcf_file + bam_file_w_index + pooltable + decodetable + reference_genome + bedfile + + main: + INDEX_VCF(vcf_file, 'GATK') + DROP_PL(INDEX_VCF.out.vcf_w_index, 'GATK') + NORMALISE_VCF(DROP_PL.out.vcf_file, reference_genome, 'GATK') + + BGZIP(NORMALISE_VCF.out.norm_vcf) + BGZIP.out.bgzf_file + .map { sample, vcf, index -> tuple(vcf, index) } + .collect() + .set { vcf_files } + MPILEUP(bam_file_w_index, reference_genome, bedfile) + RESCUE(pooltable, decodetable, vcf_files, MPILEUP.out.mpileup_file.collect()) + + emit: + rescue_probabilities = RESCUE.out.probabilities +} \ No newline at end of file