diff --git a/.gitignore b/.gitignore index 7ee95874b46de271aec64960fd553288ba33bce8..b90f8635f8c5b8775e46b5e3eb959404e2c3f6a0 100644 --- a/.gitignore +++ b/.gitignore @@ -25,6 +25,9 @@ _f6_validation/ *.mmi *.bed *.sh +!references/numt/**/*.bed +!src/bin/*.sh core code vscode_cli.tar.gz +desktop.ini diff --git a/README.md b/README.md index 452a31bd3b4d127587124149626e59ce5120a562..e1f05ece660cfed080eb08b2ccbd8b9a96e26a22 100644 --- a/README.md +++ b/README.md @@ -122,6 +122,8 @@ modkit produces the per-sample bedMethyl pileup (bgzip + tabix, `--allow-non-pri - `motif`: per-sample de novo motif discovery for each mod code, plus the GATC 6mA frequency as a tested hypothesis. - mtDNA metrics master table: the per-sample step-2 metrics joined to the samplesheet groups. - Haplogroup identity QC (gated by `--enable_haplogroup`): Clair3 calls homoplasmic mtDNA variants, Haplogrep3 assigns the haplogroup, and the results are merged into `haplogroups_master.tsv`. This is an identity, sample-swap and NUMT-contamination check, not heteroplasmy calling: the same subject should resolve to the same haplogroup across runs. +- Heteroplasmy calling (gated by `--enable_heteroplasmy`): mutserve2 (from the mtDNA-Server 2 container) calls low-frequency mtDNA variants on the chrM reads, both at native depth and at a coverage-matched depth (`--het_downsample_target`), with optional Haplocheck sample-contamination QC (`--enable_haplocheck`). Per-sample variant+AF tables are merged into `heteroplasmy_variants_master.tsv` for downstream group analysis in R. Unlike the haplogroup step, this is purpose-built for low-level (instability) variation. +- Heteroplasmy scan, substitution spectrum and per-molecule phasing (gated by `--enable_mt_variant_scan`, default off): `samtools mpileup` allele counts per strand, heteroplasmic and homoplasmic site calls with hotspot, homopolymer, barcode cross-talk and two-individual mixture flags, the 12- and 6-class spectrum with depth-adjusted group tests, and Fisher phasing of minor alleles on the same molecules. Runs from an R container (`--mt_variant_r_sif`). A cross-sample NanoComp comparison is available behind `--enable_nanopack`. MultiQC then aggregates everything into a single report. @@ -315,6 +317,36 @@ Many of the parameters for this step are based on dorado basecaller, see their [ ``` +```txt +--rescue_mt + +.mt.bam") instead of EXTRACT_MT. Must be set in both the step-2 and the step-3 parameter files. Default: false> +``` + +```txt +--numt_bed + +/numt_loci_.bed", written by src/bin/derive_numt_loci.sh). Required when --rescue_mt is true; the run aborts before any process starts if it is unset. Default: null> +``` + +```txt +--rescue_max_clip + + +``` + +```txt +--rescue_max_numt_frac + + +``` + +```txt +--rescue_force + + +``` + ### Step 3: Methylation Calling, mtDNA Analysis and MultiQC ```txt @@ -347,6 +379,42 @@ Many of the parameters for this step are based on dorado basecaller, see their [ /references/mt_mdp.bed"> ``` +```txt +--rescue_mt + + +``` + +```txt +--modkit_filter_mode + + +``` + +```txt +--modkit_filter_percentile + + +``` + +```txt +--modkit_filter_threshold + + +``` + +```txt +--modkit_num_reads + + +``` + +```txt +--modkit_min_coverage + + +``` + ```txt --enable_haplogroup @@ -371,6 +439,84 @@ Many of the parameters for this step are based on dorado basecaller, see their [ ``` +```txt +--enable_heteroplasmy + + +``` + +```txt +--enable_haplocheck + + +``` + +```txt +--enable_mt_variant_scan + + +``` + +```txt +--mt_variant_r_sif + + +``` + +```txt +--mt_variant_min_depth / --mt_variant_min_vaf / --mt_variant_max_vaf / --mt_variant_min_alt_strand / --mt_variant_min_cover + + +``` + +```txt +--mtdnaserver_sif + + +``` + +```txt +--het_detection_level + +=1%; the 3/5/10% floors are applied downstream in R. Default: "0.01"> +``` + +```txt +--het_downsample_target + + +``` + +```txt +--het_downsample_seed + + +``` + +```txt +--mt_keep_supplementary + + +``` + +```txt +--het_mapq + +=10 filter (effective gate = 20). Default: 20> +``` + +```txt +--het_baseq + + +``` + +```txt +--het_strand_bias + + +``` + > !Note. The following parameters can greatly influence the memory usage. If you're running into `Out of Memory` (OOM) issues, you might want to lower one or more of the default values. Conversely, you can trade memory for cpu overhead by increasing such values. ```txt @@ -395,6 +541,12 @@ Many of the parameters for this step are based on dorado basecaller, see their [ The pipeline also supports running [pre-configured parameter files](https://www.nextflow.io/docs/latest/cli.html#pipeline-parameters). The currently supported analyses are under the `parameters/` subdir and can be used via the option `-params-file "./parameters//.yaml"` in `nextflow run`. All such files make assumptions about the type of data to be used and where they are being stored. +Three project directories are maintained, each holding `basecall.yaml` (step 1), `qc.yaml` (step 2) and `modkit.yaml` (step 3): + +- `parameters/human_blood/`: GRCh38 primary assembly with the Ensembl `MT` contig and the default (human) region BEDs, all from `nextflow.config`; `rescue_mt: true` with `references/numt/grch38/numt_loci_grch38.bed` in both `qc.yaml` and `modkit.yaml`, and `modkit_filter_mode: "fixed"` at 0.75 so the rescued pileups stay comparable with the archived runs. +- `parameters/letizia_mouse/`: GRCm39 (`mice.GRC39m.genome.fa`, contig `chrM`) with the mouse region BEDs (`*_mouse.bed`), `modkit_min_coverage: 3`, the default percentile threshold, `rescue_mt: true` with `references/numt/grcm39/numt_loci_grcm39.bed`; haplogroup, heteroplasmy and haplocheck are off. +- `parameters/nicotine_hippo_mm10/`: mm10 (`mm10/mm10.fa`, contig `chrM`) with the mouse region BEDs, `modkit_min_coverage: 3`, the default percentile threshold and the rescue off in the pipeline (the mm10 hippocampus rescue was run with the standalone `src/bin/rescue_mt.sh`, see `docs/pipeline_notes.md#numt_rescue`); haplogroup, heteroplasmy and haplocheck are off. + [top](#table-of-contents) ## Pipeline output directory @@ -432,6 +584,8 @@ Output is organised per run under three roots: `results//basecalling`, ` 1. `motif`: per-sample de novo motif discovery (plus the GATC 6mA check). 1. `nanocomp` *(gated by `--enable_nanopack`)*: cross-sample/run NanoComp comparison. 1. `haplogroup` *(gated by `--enable_haplogroup`)*: per-sample `_haplogroup.txt` and `haplogroups_master.tsv` (mtDNA identity QC), with the Clair3 variant calls under `haplogroup/clair3/`. +1. `heteroplasmy` *(gated by `--enable_heteroplasmy`)*: `heteroplasmy_variants_master.tsv` plus subdirs `mt_bams/` (extracted chrM), `matched_bams/` and `logs/` (coverage-matched BAMs + per-sample downsample fractions), `mutserve/{native,matched}/` (per-sample VCFs), and `haplocheck.txt` (sample contamination, gated by `--enable_haplocheck`). +1. `mt_variants` *(gated by `--enable_mt_variant_scan`)*: `counts/` per-sample allele counts, `sites_called.tsv`, `hotspots.tsv`, `spectrum_per_sample.tsv`, `spectrum12_per_group.tsv`, `spectrum6_per_group.tsv`, `group_tests.tsv`, `sites/` per-sample BEDs, `readbases/`, `phasing_pairs.tsv`, `phasing_summary.tsv` and the PNG figures. 1. `multiqc`: the final `multiQC_report.html` and its data directory, aggregating all of the QC inputs above. [top](#table-of-contents) @@ -572,6 +726,9 @@ The pipeline can be executed in an HPC environment using [Slurm](https://slurm.s - Variant Calling and Haplogroup - [Clair3](https://github.com/HKU-BAL/Clair3) (run from `images/clair3.sif`, not bundled in the main container) - [Haplogrep3](https://github.com/genepi/haplogrep3) (phylotree `phylotree-rcrs@17.2`) + - [mtDNA-Server 2](https://github.com/genepi/mtdna-server-2) (heteroplasmy; run from `images/mtdna-server-2.sif` = `quay.io/genepi/mtdna-server-2:v2.1.16`, not bundled in the main container) + - [mutserve](https://github.com/seppinho/mutserve) (low-frequency mtDNA variant caller, bundled in the mtDNA-Server 2 image) + - [Haplocheck](https://github.com/genepi/haplocheck) (sample contamination, bundled in the mtDNA-Server 2 image) - Other Genomics Tools - [Samtools](https://github.com/samtools/samtools) - [bedtools](https://github.com/arq5x/bedtools2) diff --git a/containers/debian-nanopore.def b/containers/debian-nanopore.def index 91bd55143ec823af1d39db077d95ae546c4f615c..58364fd4b792b04d08560b6e3b1a9581ef244411 100644 --- a/containers/debian-nanopore.def +++ b/containers/debian-nanopore.def @@ -2,8 +2,6 @@ Bootstrap: docker From: debian:12 %post - # Install all basic packages and get update - # using clean and rm at the end just to clean some temporary files apt-get update && apt-get install -y \ locales \ wget \ @@ -20,47 +18,37 @@ From: debian:12 default-jre-headless \ && apt-get clean && rm -rf /var/lib/apt/lists/* - # Install latest jq from GitHub JQ_URL="https://github.com/jqlang/jq/releases/download/jq-{{ JQ_VERSION }}/jq-linux-amd64" wget -O /usr/local/bin/jq "$JQ_URL" chmod +x /usr/local/bin/jq - # Set timezone and language for container ## ln -fs /usr/share/zoneinfo/America/Sao_Paulo /etc/localtime locale-gen en_US.UTF-8 echo 'export LANGUAGE="en_US.UTF8"' >> "$SINGULARITY_ENVIRONMENT" echo 'export LANG="en_US.UTF8"' >> "$SINGULARITY_ENVIRONMENT" echo 'export LC_ALL="en_US.UTF8"' >> "$SINGULARITY_ENVIRONMENT" - # Install python packages - # using no-cache-dir so we dont keep a copy of the downloaded package - # and break-system-packages to override PEP 668 (which blocks pip installs) - # 1) NumPy/pandas compatíveis com pycoQC + # see docs/pipeline_notes.md#containers python3 -m pip install --no-cache-dir --break-system-packages \ "numpy<2.0.0" "pandas<2.2.0" "scipy=={{ SCIPY_VERSION }}" - # 2) demais libs python3 -m pip install --no-cache-dir --break-system-packages \ pod5=={{ POD5_VERSION }} \ multiqc=={{ MULTIQC_VERSION }} \ plotly=={{ PLOTLY_VERSION }} \ ont-fast5-api - # Clone and install pycoQC from duceppemo fork - # (change plotly version to match MultiQC plotly version >= 5.18) + # see docs/pipeline_notes.md#containers cd /opt - # clone from the CIACD fork git clone https://gmapsrv.pucrs.br/gitlab/ccd-public/pycoQC.git cd pycoQC python3 -m pip install --no-cache-dir --break-system-packages . - # 3) NanoPack (NanoPlot/NanoComp) -- repeat the numpy/pandas caps so pip can't upgrade them past - # what pycoQC/MultiQC need. Unpinned; pin later via `pip show NanoPlot` if you want reproducibility. + # see docs/pipeline_notes.md#containers python3 -m pip install --no-cache-dir --break-system-packages \ "numpy<2.0.0" "pandas<2.2.0" \ NanoPlot NanoComp - # Install modkit cd /opt mkdir -p modkit cd modkit @@ -71,7 +59,6 @@ From: debian:12 test -x "$MODKIT_BIN" ln -sf "$MODKIT_BIN" /usr/local/bin/modkit - # Install Dorado cd /opt mkdir -p dorado cd dorado @@ -82,32 +69,21 @@ From: debian:12 test -x "$DORADO_BIN" ln -sf "$DORADO_BIN" /usr/local/bin/dorado - # Install minimap2 (standalone, used by INDEX_REFERENCE to build the .mmi; - # same version dorado's embedded aligner uses, so the index is compatible) + # see docs/pipeline_notes.md#containers cd /opt wget -O /tmp/minimap2.tar.bz2 "https://github.com/lh3/minimap2/releases/download/v{{ MINIMAP2_VERSION }}/minimap2-{{ MINIMAP2_VERSION }}_x64-linux.tar.bz2" - # --no-same-owner: the minimap2 tarball preserves the maintainer's uid/gid, which a - # fakeroot build cannot map -> "Cannot change ownership ... Invalid argument". Extract - # as the build user and ignore the archived ownership. tar --no-same-owner -xjf /tmp/minimap2.tar.bz2 -C /tmp mv "/tmp/minimap2-{{ MINIMAP2_VERSION }}_x64-linux/minimap2" /usr/local/bin/minimap2 chmod +x /usr/local/bin/minimap2 rm -rf /tmp/minimap2.tar.bz2 "/tmp/minimap2-{{ MINIMAP2_VERSION }}_x64-linux" - # Install mosdepth (fast coverage; binary release) wget -O /usr/local/bin/mosdepth "https://github.com/brentp/mosdepth/releases/download/v{{ MOSDEPTH_VERSION }}/mosdepth" chmod +x /usr/local/bin/mosdepth - # Install cramino (fast long-read BAM QC; binary release) - # NB: cramino tags have NO leading "v" (e.g. 1.4.0). Use the STATIC MUSL build: - # the glibc asset (cramino-linux) is built against GLIBC 2.38 > debian 12's 2.36 and - # fails at runtime ("GLIBC_2.38 not found"); cramino-linux-musl is statically linked. + # see docs/pipeline_notes.md#containers wget -O /usr/local/bin/cramino "https://github.com/wdecoster/cramino/releases/download/{{ CRAMINO_VERSION }}/cramino-linux-musl" chmod +x /usr/local/bin/cramino - # Install Haplogrep3 (mtDNA haplogroup classification, F6; Java app -> uses the JRE installed above). - # NB: Clair3 is intentionally NOT here -- it ships TensorFlow + a conda stack and would bloat this - # image; it runs from its own official image (images/clair3.sif) via a per-process container override. cd /opt wget -O hg3.zip "https://github.com/genepi/haplogrep3/releases/download/v{{ HAPLOGREP3_VERSION }}/haplogrep3-{{ HAPLOGREP3_VERSION }}-linux.zip" unzip hg3.zip -d haplogrep3 @@ -115,16 +91,11 @@ From: debian:12 HG3_BIN="$(find /opt/haplogrep3 -type f -name haplogrep3 | head -n 1)" test -n "$HG3_BIN" chmod +x "$HG3_BIN" - # wrapper (not a symlink) so the launcher resolves its own libs by its real path + # see docs/pipeline_notes.md#containers printf '#!/bin/bash\nexec "%s" "$@"\n' "$HG3_BIN" > /usr/local/bin/haplogrep3 chmod +x /usr/local/bin/haplogrep3 - # Pre-cache the phylotree INTO the image (F6): at runtime /opt is read-only, so Haplogrep3 cannot - # download trees (it caches to /opt/haplogrep3/trees/, the jar's install dir). The build node IS - # writable + online, so this warm-up classify bakes {{ HAPLOGREP_TREE }} into the image -> offline-safe - # AND version-pinned. Must run from the install dir (it reads data/ + haplogrep3.yaml). NO `|| true`: - # if the tree can't be fetched, FAIL THE BUILD here instead of shipping an image that silently sinks - # F6 at runtime. test -s on the warm-up output proves the tree actually resolved AND classify works. + # see docs/pipeline_notes.md#containers cd /opt/haplogrep3 ./haplogrep3 classify --tree {{ HAPLOGREP_TREE }} --in data/examples/example-wgs.vcf --out /tmp/_warm.txt test -s /tmp/_warm.txt @@ -136,7 +107,6 @@ From: debian:12 export LC_ALL="en_US.UTF-8" %test - # Check if installations are on path and display their versions dorado --version modkit --version minimap2 --version diff --git a/docs/pipeline_notes.md b/docs/pipeline_notes.md index 6b10e646ff3b5af3dfa7305e91aec290625e3d95..3117e933568adf7e48e489ccb526024016a11f90 100644 --- a/docs/pipeline_notes.md +++ b/docs/pipeline_notes.md @@ -7,6 +7,7 @@ This file documents the behavior, inputs, gating and gotchas of each Nextflow sc - [Entry and configuration](#entry-and-configuration) - [main](#main) - [nextflow](#nextflow) + - [containers](#containers) - [Step 1 — basecalling](#step-1--basecalling) - [basecalling](#basecalling) - [index_reference](#index_reference) @@ -21,7 +22,10 @@ This file documents the behavior, inputs, gating and gotchas of each Nextflow sc - [samtools_qc](#samtools_qc) - [mosdepth](#mosdepth) - [nanoplot](#nanoplot) + - [mt_coverage_profile](#mt_coverage_profile) + - [numt_rescue](#numt_rescue) - [Step 3 — methylation, mtDNA analysis and report](#step-3--methylation-mtdna-analysis-and-report) + - [extract_mt](#extract_mt) - [modkit_and_multiqc](#modkit_and_multiqc) - [modkit](#modkit) - [modkit_dmr](#modkit_dmr) @@ -31,6 +35,8 @@ This file documents the behavior, inputs, gating and gotchas of each Nextflow sc - [mtdna_metrics](#mtdna_metrics) - [nanocomp](#nanocomp) - [haplogroup](#haplogroup) + - [heteroplasmy](#heteroplasmy) + - [mt_variant_scan](#mt_variant_scan) ## Entry and configuration @@ -44,6 +50,8 @@ Step 1 (basecalling and alignment) needs a run identifier: it uses `--run_id` (e Step 3 (methylation calling and MultiQC) reads filtered BAMs from `bam_filtering/*-Filtered*.bam` and their indexes from `*-Filtered*.bam.bai` in the step-3 input directory. If those globs match nothing the corresponding channels are empty and MODKIT_AND_MULTIQC receives no work. The `.bai` ids are normalised by stripping a trailing `.bam` so they line up with the BAM baseName. The mito-only reference for modkit is reconstructed as `/.fa`, the same file INDEX_REFERENCE published in step 1, which avoids loading the full genome for the methylation context. F4/F5 grouping (DMR and entropy) is driven by `--samplesheet`, which also supplies the labels for the metrics join; the per-sample mtDNA metric TSVs come from `intermediate_qc_reports/mtdna_metrics/` published by step 2. +Two rescue-related behaviours also live in this file. In step 2, if `rescue_mt` is set and `numt_bed` is null, the run aborts with an error before any process starts. In step 3, when `rescue_mt` is on, the rescued MT-only BAMs published by step 2 are read from `bam_filtering_mt/*.mt.bam` with the file `simpleName` as key (the real_barcode, since the file is `.mt.bam`) and joined with their `.bai` by that key rather than by emission order, so a missing or extra index cannot shift the pairing; `MODKIT_AND_MULTIQC` then uses that channel as `mt_bam` instead of running `EXTRACT_MT` (`#numt_rescue`). + Load-bearing gotcha: BAMs and BAIs are paired with a keyed `join(..., by: 0)` on the id, not by zipping two separately-sorted channels in emission order. If this is changed to rely on emission/sort order, a BAM can be paired with the wrong sample's index whenever the `.bam` and `.bam.bai` listings sort differently, which silently produces incorrect methylation results rather than failing. ### nextflow @@ -52,6 +60,8 @@ File: `src/nextflow.config` Main Nextflow configuration for the nanopore mtDNA pipeline. Per-step queue settings live in `src/configs` and are pulled in via `includeConfig` based on `params.step` (step 1 uses `queue-basecalling.config`, everything else `queue-default.config`). +Comment policy for this repository: the code carries no comments. A non-obvious choice is marked in the code by a one-line pointer of the form `see docs/pipeline_notes.md#`, and the rationale behind the choice lives here, under the section the pointer names; a pointer is never expanded back into an inline comment, and the section here is updated in the same change that alters the choice it explains. + `cleanup = true` is the production/full-run setting (500 GB disk cap): Nextflow removes `work/` as each step succeeds, so intermediates do not pile up across step1->step2->step3. This is safe because every `publishDir` uses `mode "copy"`, so `results/` are real copies independent of `work/`. The trade-off is that `-resume` of an already-successful step re-runs it (a failed step keeps its `work/`, so resume-on-failure still works). Set it `false` when debugging a run. Path layout: when `run_id` is set, inputs come from `data_root/raw/` and outputs go under `results_root//{basecalling,qc,modkit}`; otherwise inputs are `data_root/raw` and outputs go under `results_/...`. `steps_2_input_directory`/`steps_3_input_directory` are only populated when `run_id` is set (null otherwise), wiring each step to the previous step's output. @@ -62,7 +72,9 @@ Gating flags: `enable_nanopack` (F5 NanoPlot/NanoComp) stays `false` until the c Other parameters: `mapq` is the MAPQ filter threshold for BAMs (0 = no filtering); `qscore_thresh` is the quality-score cutoff; `basecall_speed` selects fast/hac/sup (`@latest` for the newest available); `barcoding_kit` is the barcoding kit (null skips `--kit-name`); `basecall_config` is an explicit model path or null to derive from speed+mods; `basecall_trim` is the basecalling trim mode ("all"/"primers"/"adapters"/"none", set "none" to not trim during basecalling); with `basecall_demux = true`, `trimmed_barcodes = true` makes demux only separate files while `false` makes demux trim then separate. The `gpu` process label sets `containerOptions = '--nv --nvccli'`: `--nv` is the default flag for CUDA apps (can be enabled by default in `/etc/apptainer/apptainer.conf`), and `--nvccli` uses the nvidia-container-cli. -NUMT control: rather than a methylation-based classifier (since removed), reads are aligned to the whole human genome (GRCh38) so NUMT-derived reads compete for their true nuclear loci; only authentic mtDNA stays on the mito contig, which is what modkit (`--region ${mt_contig}`) analyses. +NUMT control: rather than a methylation-based classifier (since removed), reads are aligned to the whole genome so NUMT-derived reads compete for their true nuclear loci; only authentic mtDNA stays on the mito contig, which is what modkit (`--region ${mt_contig}`) analyses. **Known failure mode (found 2026-09-05 on the mouse hippocampus data):** where a nuclear copy is near-identical to mtDNA, the competition goes the other way. Authentic mtDNA reads contained in that segment tie with the nuclear copy, get MAPQ 0, and the step-2 MAPQ filter removes them from the mito contig. In mouse the chr1 copy of chrM 6,393-11,042 (99.96 percent identity on GRCm39; mm10 chr1:24,611,534-24,616,177) erased a quarter of chrM (CO2 to ND4L) in every hippocampus sample and cut the Letizia cohort to half depth over CO2 to ATP6. The per-window chrM coverage profile (`#mt_coverage_profile`) is therefore a mandatory QC after step 2, and the opt-in step-2 rescue (`#numt_rescue`, `params.rescue_mt`) rebuilds the mito contig from the unfiltered BAM. Human GRCh38 has no copy close enough to matter (best 98.55 percent, 1:629,083-634,924); the check was run and is recorded in the human blood compendium. + +Rescue parameters (`#numt_rescue`): `rescue_mt` (opt-in), `numt_bed` (padded NUMT loci for the build), `rescue_max_clip` (300 bp of unaligned sequence drops a read), `rescue_max_numt_frac` (0.01, stop rule) and `rescue_force`. `modkit_num_reads` (50000) makes every MT read enter the modkit threshold estimate, and `modkit_filter_mode` with `modkit_filter_threshold` selects the percentile or the fixed pass threshold (`#modkit`). Load-bearing gotchas: - `mt_contig` must match both the reference FASTA and the dorado-aligned BAM contig name (`MT` in GRCh38/Ensembl, `chrM` in UCSC). It is consumed by modkit `--region`, INDEX_REFERENCE (MT.fa extraction) and coverage; a mismatch silently yields empty mtDNA output. @@ -70,6 +82,18 @@ Load-bearing gotchas: - `haplogrep_tree` must match the `HAPLOGREP_TREE` value baked into the container (`containers/versions.txt`); drift breaks or mis-calls F6 haplogroup QC. - `basecall_mods` cannot list more than one modification per nucleotide (valid examples: `4mC_5mC`, `5mCG_5hmCG`, `5mC_5hmC`, `6mA`). +### containers + +File: `containers/debian-nanopore.def`, built with `--build-arg-file containers/versions.txt` (README, Getting Started); the `{{ NAME }}` placeholders in the recipe are the variables of `versions.txt`, which is also where every tool version is pinned. + +Python stack. The first pip call installs `numpy<2.0.0`, `pandas<2.2.0` and `scipy=={{ SCIPY_VERSION }}` together: pycoQC, and the MultiQC stack it shares plotly with, do not work with numpy 2 or pandas 2.2, and pinning scipy in the same call makes pip resolve the three against each other instead of upgrading one of them in a later call. Every pip call in the recipe uses `--no-cache-dir`, so no copy of the downloaded wheels is kept in the image, and `--break-system-packages`, which overrides PEP 668; without it pip refuses to install into Debian 12's system Python. pycoQC is not installed from PyPI: the CIACD GitLab fork (`https://gmapsrv.pucrs.br/gitlab/ccd-public/pycoQC.git`, derived from the duceppemo fork) is cloned into `/opt` and pip-installed from source, because that fork raises pycoQC's plotly requirement to `>= 5.18` so it matches the plotly version MultiQC needs (`PLOTLY_VERSION`); the README names the fork but not this reason. NanoPlot and NanoComp (NanoPack) are installed in a separate, later pip call that repeats the `numpy<2.0.0` and `pandas<2.2.0` caps, so that resolving NanoPack's dependencies cannot upgrade those two past what pycoQC and MultiQC need. NanoPlot and NanoComp themselves are unpinned (the README's version note says so); if reproducibility requires it, pin them later from `pip show NanoPlot` and `pip show NanoComp` inside the built image. + +Binaries. A standalone minimap2 binary is installed even though dorado embeds its own aligner: `INDEX_REFERENCE` (`#index_reference`) uses it to build the `.mmi`, and `MINIMAP2_VERSION` is kept at the same version as dorado's embedded minimap2 so the index stays compatible with the aligner that reads it. The minimap2 tarball preserves the maintainer's uid/gid, which a fakeroot build cannot map and fails on with `Cannot change ownership ... Invalid argument`; tar is therefore run with `--no-same-owner`, also applied to the modkit and dorado tarballs, so archives are extracted as the build user and the archived ownership is ignored. cramino release tags carry no leading `v` (`1.4.0`, not `v1.4.0`), and the static musl asset (`cramino-linux-musl`) is downloaded instead of the glibc asset (`cramino-linux`) because the latter is built against GLIBC 2.38, newer than Debian 12's 2.36, and fails at runtime with `GLIBC_2.38 not found`. + +Haplogrep3. It is a Java application, which is why `default-jre-headless` is in the apt package list and why `%test` runs `java -version`. `/usr/local/bin/haplogrep3` is a bash wrapper script that `exec`s the real launcher path, not a symlink like modkit and dorado, because the Haplogrep3 launcher resolves its own libraries relative to its real path and a symlink would not preserve that. The phylotree named by `HAPLOGREP_TREE` is baked into the image by running a warm-up `haplogrep3 classify` on the bundled `data/examples/example-wgs.vcf` during `%post`: at runtime `/opt` is read-only, so Haplogrep3 cannot download trees then (it caches them to `/opt/haplogrep3/trees/`, its install directory), whereas the build node is writable and online. The warm-up must run from the install directory, because Haplogrep3 reads `data/` and `haplogrep3.yaml` relative to it, and it is deliberately not guarded with `|| true`: an unreachable tree fails the build instead of shipping an image that silently sinks the F6 haplogroup step at runtime, and `test -s` on the warm-up output plus `test -d /opt/haplogrep3/trees` prove that the tree resolved and that classify works. This is what lets `HAPLOGROUP` run offline with a pinned tree version; `params.haplogrep_tree` must match `HAPLOGREP_TREE` (`#nextflow`, `#haplogroup`). + +Not in this image. Clair3 is intentionally not built in, because it ships TensorFlow plus a conda stack that would bloat the image; it runs from its own official image (`images/clair3.sif`, `params.clair3_sif`) through a per-process container override (`#haplogroup`, and the README build section). The heteroplasmy arm likewise runs mutserve and haplocheck from the mtDNA-Server 2 image (`params.mtdnaserver_sif`, `#heteroplasmy`), and the R stages of the variant scan run from `params.mt_variant_r_sif` (`#mt_variant_scan`), because the main image has no R. + ## Step 1 — basecalling ### basecalling @@ -137,6 +161,8 @@ NanoPlot is gated behind `params.enable_nanopack`: it only runs when that flag i Channel pairing gotcha: the filtered/total BAM and BAI channels are each re-keyed by the sample's real_barcode (the first token of the filename before `-`) and combined with `.join()` rather than relying on emission order. This protects the QC fan-out from misalignment: if any sample is dropped by the `min_mapped_reads` filter, positional/order-based channel combination would silently shift labels so an id pairs with the wrong sample's BAM. Keep the per-file keying plus `.join()`; do not replace it with positional zipping. +The sub-workflow also takes `mt_reference` (the `.fa` written by `INDEX_REFERENCE`) and `numt_bed`. When `params.rescue_mt` is true it feeds `RESCUE_MT` with the unfiltered total BAM and BAI (`tbam_keyed.join(tbai_keyed)`, the same real_barcode keying), not the MAPQ-filtered pair, because the step-2 MAPQ filter is exactly what removes the NUMT-tied authentic mtDNA reads the rescue rebuilds the mitochondrial contig from; `RESCUE_MT_SUMMARY` then collects the per-sample summaries (`#numt_rescue`). + ### filtering_and_qc_from_minknow File: src/sub_workflows/FILTERING_AND_QC_FROM_MINKNOW.nf @@ -149,6 +175,8 @@ The FILTER_BAM call pairs each BAM with its summary text by sorting both channel No load-bearing inline comments were present. +That statement and the input list above predate the rescue branch. The sub-workflow now also takes `mt_reference` and `numt_bed`, and when `params.rescue_mt` is true it branches `FILTER_BAM`'s unfiltered BAM and BAI (`total_bam`/`total_bai`, each keyed by the token before the first `-` of the file name and joined on that key) into `RESCUE_MT`, with `RESCUE_MT_SUMMARY` collecting the per-sample summaries (`#numt_rescue`). + ### convert_input_from_minknow File: `src/modules/convert_input_from_minknow.nf` @@ -189,6 +217,8 @@ Two processes that run pycoQC over the aligned BAMs to produce per-sample HTML/J A sample whose pycoQC run fails produces no QC plot for that sample, but the sample is not dropped from the pipeline: it still flows on to mtdna_metrics/modkit. This is enforced by `errorStrategy 'ignore'`, which is load-bearing. A near-empty barcode yields an empty sequencing_summary, and pycoQC's pandas `read_csv` then throws an `EmptyDataError`. With `'ignore'` that failure removes only that one sample from the pycoQC chain instead of aborting the whole step-2 QC. Do not change `errorStrategy` to `'terminate'` or `'retry'` here, or a single empty/near-empty sample will sink QC for the entire run. +Both processes, PYCOQC_NO_FILTER and PYCOQC_FILTER, carry that `errorStrategy 'ignore'` for the same reason: a bad or empty sample skips its unfiltered or filtered plot without sinking the step. + ### samtools_qc File: `src/modules/samtools_qc.nf` @@ -221,8 +251,312 @@ Runs NanoPlot on a per-sample filtered BAM to produce a read-length and quality The process is gated by `params.enable_nanopack`. It needs NanoPlot installed in `containers/debian-nanopore.def`, which is not yet present, so the step is effectively unavailable until that container dependency is added. The output is declared `optional: true`, so when nothing is produced the channel is simply empty rather than failing the run. +### mt_coverage_profile + +File: `src/bin/mt_coverage_profile.R`. Standalone Rscript, run after step 2 on every mtDNA +project; not wired into Nextflow yet. It reads every `*_modkit_pileup.bed.gz` in a directory +through tabix, averages Nvalid_cov of one modification code (default `a`, one row per covered +adenine per strand) per window on the mitochondrial contig, writes a per-window table, a +per-sample gap table and a figure, and exits with status 2 when some window is below the +minimum in every sample. A gap shared by all samples is the signature of a near-identical +NUMT dropping authentic mtDNA reads at the MAPQ filter (mouse chr1 copy, 2026-09-05); a gap +in one sample is a sample problem. A partial loss shows as a dip that never reaches the +minimum, so the figure must be read as well as the exit status: on the Letizia cohort the +tool exits 0 while the profile halves over CO2 to ATP6 (`#numt_rescue`). The sample id is +the file name before `_modkit_pileup.bed.gz`, truncated at the first dash, which is the +pipeline convention (`-Filtered_primary_mapq_10`), so the samplesheet join +works on pipeline output as well as on plain `_modkit_pileup.bed.gz`; `--run-id` filters +the samplesheet to one run when `real_barcode` values repeat across runs. Runs in the rocker +container with the study library (data.table, Rsamtools, GenomicRanges, ggplot2). Examples: + + Rscript src/bin/mt_coverage_profile.R --dir --contig chrM --length 16299 --samplesheet references/samplesheet_letizia.csv --id-col real_barcode --group-col group + Rscript src/bin/mt_coverage_profile.R --dir --contig MT --length 16569 --samplesheet references/samplesheet.csv --id-col real_barcode --group-col group --run-id run01 + +Optional `--features ` draws the gene track; `--window` (250), `--min` (3) +and `--mod` (`a`) are tunable. Results so far: mouse hippocampus MAPQ-filtered pileups report +7,001-9,750 (exit 2) and nothing on the rescued ones; Letizia and human blood run01 exit 0 +before and after the rescue (`#numt_rescue`). Still to do: add it as a step-3 pre-check that +fails the run on exit 2. + +### numt_rescue + +Files: `src/bin/derive_numt_loci.sh`, `src/bin/rescue_mt.sh`, and the verification tools +`src/bin/mt_depth_compare.sh`, `mt_site_compare.sh`, `mt_supplementary_check.sh`, +`mt_read_accounting.sh`, `mt_pileup_check.sh`; reference data in `references/numt//`. +All run inside `images/debian-nanopore.sif` and take their settings from environment +variables documented in each script header, so a Slurm launcher or a Nextflow process can +drive them without editing. Written 2026-09-05 for the mouse hippocampus data and rolled +out on 2026-09-06 to the Letizia mouse and human blood projects. + +#### The failure mode, restated + +Whole-genome alignment is the pipeline's NUMT defence (`#nextflow`, NUMT control): a nuclear +copy of mtDNA attracts its own reads, and only authentic mitochondrial reads stay on the +mitochondrial contig. The defence inverts when a copy is near-identical. A read that lies +entirely inside the copied segment aligns equally well to both places, gets MAPQ 0, and is +removed by the step-2 filter (`mapq = 10`). The loss is largest in the middle of the copied +segment, because reads near its ends reach unique mitochondrial sequence and are placed with +confidence; the coverage profile therefore shows a V-shaped dip whose depth grows with the +copy's length and identity and shrinks with read length. Found on 2026-09-05 in the mouse +hippocampus data (mm10, chr1 copy, a hole from 7,001 to 9,750 in every animal) and on +2026-09-06 in the Letizia mouse data (GRCm39, same copy, a dip to about half the flanking +depth over CO2 to ATP6). + +#### Deriving the nuclear copies (`derive_numt_loci.sh`) + +The mitochondrial contig is extracted from the genome fasta itself (`samtools faidx`), so the +sequence the rescue realigns to is exactly the one the reads were aligned to, and aligned back +against the whole genome: + + minimap2 -cx asm20 -I 500M -N 500 -p 0.01 --secondary=yes genome.fa .fa + +`-c` gives base-level identity instead of minimizer counts; `-p 0.01` and `-N 500` keep the +secondary hits (the full-length self-hit would otherwise suppress every copy); `-I 500M` +indexes the genome in batches so the job runs in 6 GB of memory (peak 5.4 GB on GRCm39) in +about one minute. Hits of at least 200 aligned bases are kept, sorted, and written three ways: +`numt_hits_.tsv` (chromosome interval, mitochondrial interval, strand, matches, +aligned length, MAPQ, identity), `numt_loci_.bed` (hits padded by 5,000 bp on each +side and merged: the candidate regions the rescue reads), and `numt_core_.bed` +(unpadded). `numt_near_identical_.tsv` lists copies of at least 1,000 bp at 98 percent +identity or more, which are the only ones able to take MAPQ from authentic reads. The md5 of +the mitochondrial sequence is compared across every fasta copy in `references/` +(`mt_contig_md5.tsv`); mouse `chrM.fa`, `mm10/chrM.fa` and the GRCm39 contig are identical +(16,299 bp), as are human `MT.fa`, `chrMT.fa` and the GRCh38 contig (16,569 bp). + +| build | hits >= 200 bp | candidate intervals (span) | near-identical copies | +|---|---|---|---| +| GRCm39 | 11 | 11 (124,926 bp) | chr1:24,650,615-24,655,265 = chrM 6,393-11,042, 99.96 percent | +| GRCh38 | 24 | 21 (269,391 bp) | 1:629,083-634,924 = MT 3,913-9,755, 98.55 percent | + +The GRCm39 list is shorter than the 21-interval mm10 list used for the hippocampus rescue +because that list came from minimizer mapping without base-level alignment, which also +reports copies below 90 percent identity. Such copies cannot capture nanopore mtDNA reads +(a read at 98 to 99 percent accuracy scores far better on the mitochondrial contig), and +including them only adds nuclear reads to the candidate set, so their absence is intended. +The GRCh38 copies overlapping ND4L/ND4 (5:134,923,308 and X:126,472,730) are at 94 percent +and do not affect that locus; the human before-profile is flat. + +#### The rescue (`rescue_mt.sh`) + +Input is the step-2 `-Unfiltered.bam` (the whole-genome alignment before the MAPQ filter, +published to `qc/bam_filtering/` and archived on Drive). Three stages, each reusing its output +when newer than its input: + +1. **Candidates.** Primary records, any MAPQ, that lie on the mitochondrial contig or on any + padded NUMT locus (`samtools view -F 0x904 -L`). Primary records only, on purpose: a + primary record carries the whole read (soft-clipped) and its MM/ML/MN tags, so each read + is extracted exactly once and the realignment regenerates its supplementary records. + Per-sample counts go to `qc/candidate_summary.tsv`: candidates from the mitochondrial + contig, from NUMT loci, from NUMT cores (unpadded), and with MAPQ 0. +2. **Realignment.** Candidates are split by origin and realigned to the mitochondrial contig + alone, keeping the origin in the read group (`.mt`, `.numt`; `SM:` as + `EXTRACT_MT` sets it): + + samtools fastq -T MM,ML,MN | minimap2 -y -Y -ax lr:hq --secondary=no -R @RG... .fa | samtools view -F 4 | samtools sort + + `samtools fastq` returns the read in its original orientation, which is the frame the + modification tags are defined in, and `-y` re-emits the tags verbatim. `-Y` keeps + supplementary records soft-clipped, which `modkit --allow-non-primary` needs for the reads + that cross the origin of the circular contig. Then, per read, aligned bases are summed over + its primary and supplementary records and the remainder counted as unaligned; reads with + `MAX_CLIP` (300) or more unaligned bases are dropped (see the filter below), the rest form + `qc/bam_filtering/-Filtered__realigned.bam`. Depth is profiled in 250-bp windows + (`qc/mt_depth/`) and the per-read table and the dropped reads are kept in `qc/candidates/`. + `qc/realign_summary.tsv` holds the diagnostics per sample. +3. **modkit**, mt scope, on the realigned BAM: `sample-probs --force`, `pileup`, `summary`, + all with `--region --num-reads 50000 --allow-non-primary` (see `#modkit` for why + `--num-reads` is required; `--force` because a cancelled run leaves `thresholds.tsv` + behind). The pass-threshold policy is a parameter, `FILTER_MODE`: `percentile` (default, + `--filter-percentile 0.1`, the pipeline's current rule, used for the mouse projects) or + `fixed` (`--filter-threshold 0.75`). The fixed mode exists because every archived human + blood run (`REBASECALL_dorado2_hac6/run01` to `run09b`, 2026-06) was called with the fixed + 0.75 threshold, and a rescued pileup can only be compared with its archived counterpart + under the same policy: on run01 the percentile rule sets thresholds of 0.86 to 0.90 and + fails 7.0 reads per adenine site against 2.9, so valid coverage drops 5 to 9 percent even + though the BAM gained reads. Outputs land in `modkit/modkit/` with the pipeline's file + names, so step 3 tools and `mt_coverage_profile.R` read them unchanged. `PROVENANCE.txt` + records tool versions, references, parameters, the filter policy and the code revision. + +**Why a filter on unaligned bases, and why 300.** Realigning NUMT-locus reads to the +mitochondrial contig alone turns a true nuclear read into a partial mitochondrial alignment: +the copied segment aligns, the nuclear flank is soft-clipped. Adapters and barcodes also +soft-clip (this pipeline basecalls with `--no-trim`), but they are short. Measured on the +realigned reads (2026-09-06): in the Letizia mice, 213 of the 219 mitochondrial-origin reads +with less than 80 percent aligned had fewer than 300 unaligned bases and 191 were shorter than +500 bp (adapters on short reads), while in human blood 49 of the 52 NUMT-origin reads below +80 percent had 300 to several thousand unaligned bases (nuclear flanks). A fraction-based +cut would therefore discard short mitochondrial reads and keep long nuclear ones; the +absolute cut removes 16 of 10,850 mouse reads (0.15 percent) and 71 of 3,380 human reads, +50 of them from NUMT loci. Reads entirely contained in a copied segment cannot be told apart +by any method; they are rare because nuclear depth at a single locus is a small fraction of +mitochondrial depth, and they are the residual risk the decision rule below is about. + +What the filter removes on the mitochondrial side was checked on the ten human runs +(2026-09-06, `mt_chimera_check.sh` and `mt_tail_detail.sh` on the whole-genome records kept in +`qc/candidates/.candidates.bam`): 278 MT-origin reads of 30,574, median length 5 to 16 kb +against 1 to 3 kb for MT reads at large, mitochondrial segments at MAPQ 60 and 94 to 99 percent +identity, tails that carry further mitochondrial pieces (163), align nowhere (73) or hit a +nuclear coordinate with a 93 to 341 bp alignment at MAPQ 1 to 20 and about 20 percent +mismatches (42), recurring at the same low-complexity positions across unrelated participants +(chr9:78,744,236, chr22:40,007,216, chr8:111,249,051). They are mitochondrial reads with noise +tails, not NUMT reads; a NUMT insertion would show a long nuclear flank at MAPQ 60. Dropping +them costs 1 to 2 percent of bases and no site-level information; whether to scope the filter +to NUMT-origin reads only is a project decision. + +**Decision rule (gate).** The standalone script stops before modkit, with exit status 3, when +the candidates from NUMT loci exceed `MAX_NUMT_FRAC` (0.01) times the candidates from the +mitochondrial contig in any sample, unless `FORCE=1`. This is the rule written after the +hippocampus rescue. The raw ratio counts the rescued mitochondrial reads together with nuclear +reads, so it cannot pass on any project where the rescue matters (hippocampus 19 percent, +Letizia 11 percent, human blood 16 percent); it works as a forced stop that puts +`realign_summary.tsv` in front of the analyst. The columns that separate the two read classes +are `from_numt_alnfrac_low` and `clip_dropped_numt` against their `_mt` counterparts: in the +mice the NUMT-origin reads distribute like the mitochondrial ones (excess of partially aligned +reads 17 of 10,850), in blood they do not (39 percent partially aligned against 1 percent). +The Nextflow process therefore gates on the second measure (`params.rescue_max_numt_frac`, +clip-dropped NUMT-origin reads as a fraction of primaries: 0.07 percent in the mice, 1.5 +percent in blood). + +**Supplementary records and the origin.** Verified with `mt_supplementary_check.sh`: Letizia +step 2 kept 616 split reads on chrM, the rescued BAM holds 609, 598 of them crossing the +origin; only 3 reads whose primary record lies outside the candidate regions (chimeras with a +mitochondrial supplementary) are not carried. Human run01: 469 before, 455 after, 451 origin +spanners, 9 lost at MAPQ 10. + +#### Results per project + +| project | before | after | claim sites | +|---|---|---|---| +| hippocampus nicotine (mm10, 8 mice) | hole 7,001-9,750 in every mouse, exit 2 | no window below 3 reads, 10,308 of 10,310 adenines callable | all 38 features testable | +| Letizia (GRCm39, 19 animals, two flowcells pooled) | V dip, exit 0 (median depth 5 to 30 per strand-site), 1,823 of 11,050 candidates at MAPQ 0 | depth inside 6,393-11,042 x1.44 (centre x2.0), outside x1.00; every adenine and cytosine covered on both strands in every animal; thresholds equal to the pipeline run | valid reads at 9,529 (minus strand) 144 to 160, 9,658 117 to 130, 12,262 unchanged; ATP6 7.6 to 15.3, CO2 8.3 to 15.2, ND3 14.3 to 16.0, ND4 15.6 to 18.3 | +| human blood run01 (GRCh38, 6 participants) | flat, exit 0, 23 of 3,764 candidates at MAPQ 0 | BAM depth flat at x0.97 to 1.00, no localised gain; 50 nuclear-like reads dropped, 83 NUMT-locus reads kept; under the fixed 0.75 threshold valid coverage within 3 percent of the archive | ND4L/ND4 unchanged | + +In human blood the best copy (98.55 percent) is resolved by whole-genome alignment at these +read accuracies, so there is no dip to repair; the rescue is neutral there (the reads it adds +are offset by the reads the clip filter removes) and is run on every run, one at a time, with +the fixed 0.75 threshold to match the archive, so that the whole mtDNA set exists in one +consistent rescued form next to the mouse projects. Acceptance per run: before- and +after-profile exit 0, every read carrying MM tags, no localised gain in the window depth +ratio, supplementary-only losses in the single digits, and thresholds equal to the archived +ones. Rescued outputs go to `gdrive:DCNL/nanopore_analyses/human_blood_chrM_rescue//`; +the archived `REBASECALL_dorado2_hac6/` results are never touched. + +#### Nextflow processes (`src/modules/rescue_mt.nf`) + +`RESCUE_MT` runs in step 2 (both entry points) when `params.rescue_mt` is true, once per +sample on the `-Unfiltered.bam` that `FILTER_BAM` publishes, keyed by real_barcode like +the other step-2 QC processes. It calls `rescue_mt.sh` from `bin/` with `STAGE=realign` +(candidates, realignment by origin, clip filter, depth profile; no modkit, because step 3's +`MODKIT` does that with the pipeline's own threshold policy) and publishes +`qc/bam_filtering_mt/.mt.bam` plus `qc/mt_rescue/.rescue_summary.tsv`, +`.depth250.tsv` and `.clip_dropped.tsv`. The mitochondrial fasta is the +`.fa` that `INDEX_REFERENCE` writes next to the genome in step 1; the NUMT loci +come from `params.numt_bed`. `RESCUE_MT_SUMMARY` concatenates the per-sample rows into +`qc/mt_rescue/rescue_summary.tsv`, prints one line per sample, and fails the run when +clip-dropped NUMT-origin reads exceed `params.rescue_max_numt_frac` of the primaries in any +sample, unless `params.rescue_force`; the stop is the moment to read the table and decide. + +In step 3, `main.nf` reads `bam_filtering_mt/*.mt.bam` (key = `simpleName`, joined with the +`.bai` by key) and `MODKIT_AND_MULTIQC` uses that channel as `mt_bam` instead of running +`EXTRACT_MT`; everything downstream (modkit, entropy, heteroplasmy, DMR) is unchanged. The +read groups carry `SM:`, as `EXTRACT_MT` sets them. Off by default; enabled for +`letizia_mouse` (percentile threshold) and for `human_blood` (`modkit_filter_mode: fixed`, +0.75, the policy of its archive), in both the `qc.yaml` and the `modkit.yaml` of each. +Load-bearing: `rescue_mt` must be set in both the step-2 and the step-3 parameter files, +or step 3 silently falls back to `EXTRACT_MT` on the filtered BAM. + +`numt_bed` is mandatory when `rescue_mt` is true: `main.nf` refuses the step-2 run with an +error before any process starts if `rescue_mt` is set and `numt_bed` is null, so a missing loci +file cannot turn into a rescue that silently realigns the mitochondrial reads alone. +`rescue_max_clip` set to 0 disables the clip filter entirely: `rescue_mt.sh` gates the filter on +`MAX_CLIP -gt 0`, copies the realigned BAM through unchanged and writes an empty +`clip_dropped.tsv`. In both step-2 entry points the rescue is fed the unfiltered BAM, never the +MAPQ-filtered one, because the MAPQ filter is what removes the reads the rescue exists to +recover. + +#### Verification tools + +- `mt_depth_compare.sh`: per-window and per-site depth from BAMs, before (unfiltered BAM at + the step-2 MAPQ) and after (realigned BAM), so the gain can be measured before any pileup + exists. +- `mt_site_compare.sh`: valid coverage and percent modified at named sites and regions from + the before and after pileups (sample ids matched by the token before the first dash). +- `mt_supplementary_check.sh`: split reads, origin spanners, and reads the rescue cannot + carry. +- `mt_read_accounting.sh`: what the unfiltered BAMs contain (mapped records, mitochondrial + fraction, read groups) and the aligned-fraction and clip-length distributions by origin, + which is how the 300-bp cut was chosen. +- `mt_pileup_check.sh`: pileup rows per modification code and strand, thresholds against the + archived summaries, per-site read accounting (valid, fail, nocall) before and after, and + the share of reads carrying MM tags. +- `mt_coverage_profile.R` (`#mt_coverage_profile`) on the rescued pileups: expect exit 0. +- `mt_deletion_scan.sh`: first-pass scan for mtDNA deletions in the realigned BAMs. A molecule + with a deletion aligns as two same-strand segments with a reference gap and a query gap near + zero, or as one alignment with a long D operation; candidates are clustered per sample on both + breakpoints and reported only with two or more supporting reads, because one junction read is + indistinguishable from a chimera, and a rough deletion fraction is derived from the depth at + the breakpoints. On 2026-09-06 neither cohort showed a cluster: the human runs gave 64 + single-read candidates over 56 barcodes (50 to 522 bp CIGAR deletions scattered over the + genome, plus six split reads whose query gaps of 175 to 297 bp mark alignment breaks in + low-quality stretches, not junctions), the Letizia mice 7; nothing at the common deletion + (m.8470 to m.13447). Detection floor is about two reads over the local depth, so 1 to 5 + percent heteroplasmy per sample at these depths. A dedicated long-read caller (Sniffles2 or + cuteSV) would add breakpoint refinement and the direct-repeat check but needs a container + addition; the scan answers first whether there is anything to call. + +#### Heteroplasmy scan, mutational spectrum and phasing (`mt_allele_counts.sh`, `mt_read_bases.sh`, `mt_spectrum_phasing.R`) + +Replaces, for now, the heteroplasmy arm whose container (`mtdna-server-2.sif`) is no longer on +the cluster, without any container change. `mt_allele_counts.sh` runs `samtools mpileup` (base +quality at least 10, BAQ off, secondary and duplicate records excluded) on each realigned BAM +and writes per-position, per-strand counts of A, C, G, T, deletions and insertions. +`mt_spectrum_phasing.R --stage call` calls a heteroplasmic site where the most frequent +non-reference allele has depth at least 20, a fraction between 2 and 90 percent and at least 3 +reads on each strand (90 percent or more is homoplasmic); flags systematic hotspots (a minor +allele at 0.5 percent or more in at least half of the cohort), homopolymer context (runs of 4 +plus one flank) and barcode cross-talk (a minor allele below 10 percent that is homoplasmic in +another sample of the same run); flags a sample as a mixture of two individuals when it carries +10 or more clean sites at 20 to 80 percent and excludes it from group statistics; tabulates the +12 strand-specific and 6 collapsed substitution classes; compares groups with Wilcoxon or +Kruskal-Wallis, a linear model with log depth, and a permutation chi-square on the pooled +spectrum; writes one BED of called sites per sample. `mt_read_bases.sh` extracts per-read bases +at those sites (`mpileup --output-QNAME`), and `--stage phase` tests every pair of sites in a +sample for co-occurrence of the minor alleles on the same molecules (Fisher's exact test on the +2 by 2 table of reads covering both; linked if p below 0.01, odds ratio above 1, at least 3 +shared reads). + +Load-bearing lessons from the first run (2026-09-06, human blood and Letizia): sites at 85 to 90 +percent are haplogroup polymorphisms depressed by error (m.1438A>G, m.2706A>G, the HVS1 poly-C +sites), so spectra of somatic change should be read on sites below 20 percent; a sample with +dozens of sites at 40 to 60 percent whose minor alleles are all linked on the molecules is a +mixture of two people (two run04 controls, both assigned haplogroup root R by the haplogroup +QC), and the haplogroup step should flag root-only calls and identical calls within a run; at 13 +to 67x (the mice) the both-strand rule puts the detection floor at 9 to 45 percent and no +heteroplasmy is callable, while the absence of homoplasmic differences is a strain identity +check. + +Launchers used on Atria live in `~/slurm_jobs/nanopore/` (`numt_loci.sbatch`, +`fetch_unfiltered_bams.sbatch`, `rescue_mt.sbatch`, `mt_coverage_profile.sbatch`, +`tool.sbatch`, `upload_results.sbatch`, `human_run_driver.sh`); they only set variables and +call the tools inside the container. Storage discipline: unfiltered BAMs are staged one run +at a time and deleted after the rescued outputs are verified on Drive (`rclone check`). + ## Step 3 — methylation, mtDNA analysis and report +### extract_mt + +File: `src/modules/extract_mt.nf` + +Process `EXTRACT_MT` pulls the mitochondrial contig (`${params.mt_contig}`) out of the whole-genome filtered BAM once per sample, sets `@RG ID/SM:`, sorts and indexes, and emits `(id, mt.bam, mt.bai)`. It runs in step 3 core (NOT gated) because the MT-only BAM is shared by the methylation modkit steps (MODKIT pileup/extract/summary + entropy) AND, when enabled, the heteroplasmy caller. Extracting once gives every modkit threshold estimate an MT-restricted input (no nuclear contamination — the original goal), avoids extracting twice, and shrinks the BAMs the downstream steps read. The whole-genome alignment is the NUMT defense, but it fails for near-identical nuclear copies (see `#nextflow`, NUMT control): check the per-window chrM coverage profile (`#mt_coverage_profile`) before trusting the extract, and use the rescue (`#numt_rescue`) when it dips. + +`mt_keep_supplementary` (default true) keeps supplementary alignments (reads spanning the circular 16569->1 origin become primary+supplementary), needed for F7-A D-loop recovery (`modkit --allow-non-primary`) and heteroplasmy origin coverage. Because the BAM is shared, this flag is GLOBAL: `false` adds `-F 0x800` and removes supplementary from methylation too. + +`samtools addreplacerg` sets both `@RG ID` and `SM` to the real barcode (the process `id`) so that the sample name carried by the MT BAM is correct in everything downstream that reads it from the read group rather than from the file name: the mutserve per-sample VCF sample column, the `bcftools merge` of the native-arm VCFs into the multi-sample VCF, and the haplocheck contamination report keyed by that sample name. Without a correct `SM` those outputs would be labelled by whatever read group dorado wrote, not by the sample. The rescue writes the same `SM:` into its read groups (`#numt_rescue`), so both sources of `mt_bam` label samples identically. + +Load-bearing gotchas: +- `errorStrategy 'ignore'` keeps one bad sample from sinking step 3, but EXTRACT_MT is now upstream of EVERYTHING MT-based, so the script guards against a `params.mt_contig` mismatch (e.g. `MT` vs `chrM`): it counts reads after `samtools view` and exits 1 if zero, surfacing the failure in the execution report instead of silently zeroing out the whole step. +- The extracted BAM keeps the whole-genome `@SQ` header. Do not "clean" non-MT `@SQ` lines: `samtools reheader` only rewrites header text and does not remap read tids, so dropping `@SQ` lines corrupts the BAM. modkit/mutserve handle the extra @SQ fine. + ### modkit_and_multiqc File: src/sub_workflows/MODKIT_AND_MULTIQC.nf @@ -235,21 +569,25 @@ MERGE_MTDNA_METRICS builds the per-sample mtDNA metrics master table and MultiQC Group labels come from the samplesheet (real_barcode -> group). The region BEDs (mt_macro, mt_genes, mt_control_elements, mt_mdp) stratify DMR and entropy; the control-elements BED gives localize its feature centers. -DMR uses the bgzip+tabix pileups keyed by real_barcode and joined to group, excluding user_unknown because that group has no trauma contrast. Entropy uses the filtered mod-BAMs, also keyed by real_barcode and joined to group, but keeps all groups (descriptive). Localize and motif search run per-sample and consume MODKIT.out.pileup directly. +DMR uses the bgzip+tabix pileups keyed by real_barcode and joined to group. Groups are binary (control vs user): trauma is no longer a grouping factor, so the former user_trauma/user_notrauma/user_unknown all collapse to 'user' via the samplesheet, and DMR runs the single control-vs-user pair. Entropy uses the filtered mod-BAMs, also keyed by real_barcode and joined to group (both groups). Localize and motif search run per-sample and consume MODKIT.out.pileup directly. + +Wiring contract. `bam_bai` pairs each filtered BAM with its BAI by a keyed join on the real_barcode token (the first `-`-delimited field of the BAM id and of the `.bai` file name), not by channel emission order, so the pairing does not depend on the order in which the two input channels emit. When `params.rescue_mt` is set, `mt_bam` is the `rescued_mt_bams` channel that step 2's `RESCUE_MT` published and `EXTRACT_MT` is not run; otherwise `EXTRACT_MT` runs on `bam_bai`. Both branches must yield `mt_bam` with the same tuple shape, `(real_barcode, mt.bam, mt.bai)`, because every downstream consumer (MODKIT, the entropy join, the heteroplasmy arms, mt_variant_scan) destructures it that way. The sentence above saying that entropy uses the filtered mod-BAMs is stale: entropy runs on the MT-only `mt_bam` channel (from `EXTRACT_MT` or the rescue), keyed by real_barcode and joined to `sample_groups`, not on the whole-genome filtered mod-BAMs (`#modkit_entropy`). Gating flags: NanoComp runs only when params.enable_nanopack is set (it is gated until the container ships NanoComp). The haplogroup identity QC (Clair3 -> Haplogrep3 -> merge) runs only when params.enable_haplogroup is set. Load-bearing notes: -- NUMT control is handled upstream via whole-genome alignment plus modkit --region ${mt_contig}. +- NUMT control is handled upstream via whole-genome alignment plus modkit --region ${mt_contig}, with the near-identical-copy caveat documented under `#nextflow` (mouse chr1 copy, 2026-09-05); with `params.rescue_mt` the MT-only BAM is the step-2 rescue output (`#numt_rescue`) and `EXTRACT_MT` is skipped. - Keep the .unique() on sample_groups. A real_barcode can appear on more than one samplesheet row for dual-barcode samples (always with the same group); without dedup the downstream .join() fans out and double-counts samples. -- The pileup-to-group join is by the real_barcode token (id.tokenize('-')[0]), not by emission order, and the user_unknown filter is intentional gating, not cosmetic. +- The pileup-to-group join is by the real_barcode token (id.tokenize('-')[0]), not by emission order. Group comes straight from the samplesheet `group` column; with binary control/user there is no user_unknown to filter (DMR runs the single control-vs-user pair). - Keep .ifEmpty([]) on the HAPLOGROUP collect so MERGE_HAPLOGROUPS always runs. If every Haplogrep3 result is empty, this still produces a header-only haplogroups_master.tsv instead of emitting no file at all. +- Strand collapsing (modkit `--combine-strands` or any equivalent) must not be re-added anywhere in the modkit chain invoked here (pileup, DMR, entropy): mtDNA 6mA/5mC is heavy/light-strand asymmetric (`#modkit_dmr`, `#modkit_entropy`). +- The native-arm VCF and TBI collects feeding MERGE_VCFS_HET are bare `.collect()` with no `.ifEmpty([])`, on purpose; `#heteroplasmy` explains why an empty native arm must skip Haplocheck rather than crash it. ### modkit File: `src/modules/modkit.nf` -Process `MODKIT` runs the modkit toolchain over a single coordinate-aware BAM and emits modification tables for the mitochondrial contig. Inputs: a `(id, bam)` tuple, the matching `bai`, the reference FASTA, and three tuning values (`modkit_threads`, `modkit_isize` interval size, `modkit_extract_qsize` extract queue size). The reference is indexed in-process with `samtools faidx` before any modkit call, so the small mito FASTA does not need a pre-built `.fai`. Region is fixed to `params.mt_contig` for every step and the filter threshold is `0.75` throughout. Output to `params.modkit_out_dir/modkit/`: the bgzipped+tabixed pileup (`*_modkit_pileup.bed.gz` plus `.tbi`, emitted as `pileup`), the read-level extract calls, the summary text, and all produced files via the `allfiles` glob. +Process `MODKIT` runs the modkit toolchain over the **MT-only BAM** from `EXTRACT_MT` (so every threshold estimate is MT-restricted, never contaminated by nuclear reads) and emits modification tables for the mitochondrial contig. Inputs: a `(id, bam, bai)` tuple, the reference FASTA, and three tuning values (`modkit_threads`, `modkit_isize` interval size, `modkit_extract_qsize` extract queue size). The reference is indexed in-process with `samtools faidx` before any modkit call. Region is still set to `params.mt_contig` on every step (belt-and-suspenders over the already-MT BAM). **Filtering: the pass threshold is estimated automatically at the 10th percentile** (`--filter-percentile ${params.modkit_filter_percentile}` on pileup/extract; `--filter-quantile` on summary — different flag name, same param/semantics) — **not the old fixed 0.75**. No `--seed`: in the human blood data every barcode has far fewer reads than modkit's default sampling size (~100-500 vs 10042 in Nanopore long-read), so the threshold is estimated from ALL reads and is deterministic without a seed. That premise does not hold for every project: the mouse hippocampus mice carry 3,000 to 16,000 chrM reads, and above 10,042 modkit would subsample the threshold estimate non-reproducibly. `sample-probs`, `pileup` and `summary` therefore also pass `--num-reads ${params.modkit_num_reads}` (default 50000), so every MT read enters the estimate whatever the depth. The policy itself is a parameter: `modkit_filter_mode = "percentile"` (default, the rule above) or `"fixed"`, which passes `--filter-threshold ${params.modkit_filter_threshold}` (0.75) to extract, pileup and summary instead; `fixed` exists so that the human blood project can be reprocessed under the same threshold as its archive (next note), because valid coverage and percent modified are only comparable under one policy. Passing `--seed` would make modkit assume fractional sampling and require `--sampling-frac` (this is what broke the first full run). Provenance note (2026-09-06): the archived human blood runs on Drive (`REBASECALL_dorado2_hac6/run01` to `run09b`, processed 2026-06-05) were produced by the previous main-line code with `--filter-threshold 0.75` (the cached task dumps in `~/proveniencia/nanopore/runs/` show 148 modkit tasks at 0.75, 18 at 0.90 and 10 at the percentile), so a re-analysis under the percentile rule changes valid coverage and percent modified for reasons unrelated to the data; compare like with like. Output to `params.modkit_out_dir/modkit/`: the bgzipped+tabixed pileup (`*_modkit_pileup.bed.gz` plus `.tbi`, emitted as `pileup`), the read-level extract calls, the summary text, and all produced files via the `allfiles` glob. Steps run in order: `sample-probs` (histogram of base-mod probabilities), `extract calls` (per-read TSV, bgzf, mapped-only and allow-non-primary), `pileup` (site-level bedMethyl), and `summary` (global per-mod/base stats). The pileup is written with `--bgzf` and then tabix-indexed so the DMR, localize, and motif steps can consume it directly. @@ -284,7 +622,7 @@ File: src/modules/modkit_entropy.nf Computes per-group methylation entropy (heterogeneity within sliding windows) as a dysregulation readout that needs no cross-sample absolute calibration. Entropy is computed per strand because mtDNA 6mA/5mC could be heavy/light asymmetric. `--base A` measures 6mA; `--base C` measures all cytosine mods together (entropy cannot split 5mC from 5hmC). Results are stratified over the supplied region BEDs. -Inputs: all group-eligible step-2 filtered mod-BAMs with MM/ML tags (collected list) plus their `.bai` indexes; a `manifest` TSV with one line per BAM as `\t`; a small mito reference `MT.fa`; and a collected list of region BEDs. Distinct groups are derived from column 2 of the manifest, and one `-s ` is passed per BAM in that group so modkit aggregates entropy across all of the group's mod-BAMs. The reference `.fai` is built at the start with `samtools faidx` because modkit reads it alongside the reference. +Inputs: all group-eligible **MT-only mod-BAMs from `EXTRACT_MT`** (MM/ML tags preserved; collected list) plus their `.bai` indexes; a `manifest` TSV with one line per BAM as `\t`; a small mito reference `MT.fa`; and a collected list of region BEDs. Distinct groups are derived from column 2 of the manifest, and one `-s ` is passed per BAM in that group so modkit aggregates entropy across all of the group's mod-BAMs. The reference `.fai` is built at the start with `samtools faidx` because modkit reads it alongside the reference. Output: `entropy/**` (one output directory per group/region/base) and `*.log`, both declared `optional: true`. A group/region with insufficient coverage produces no entropy output for that combination. @@ -295,7 +633,7 @@ Load-bearing gotchas: - The shell array must be named `GROUP_LIST`, never `GROUPS`. `GROUPS` is a special bash variable whose assignments have no effect (bash manual), so `GROUPS=(...)` is silently ignored and the per-group loop would run over nothing with no error. Command substitution is used rather than process substitution `< <(...)` because Apptainer has no `/dev/fd`, and `|| true` guards an empty manifest under `set -e`. - Do not add `--combine-strands`. The heavy and light strands must stay separate, since they are biologically asymmetric for 6mA/5mC. - `--regions` requires `--out-bed` to be a directory, which is why each group/region/base gets its own `outdir`. Do not repoint `--out-bed` at a file. -- `--filter-threshold 0.75` must match the value used in the pileup step so entropy and pileup share the same call set. `--min-coverage 3` is the modkit default. +- Filtering uses the auto 10th-percentile threshold (`--filter-percentile ${params.modkit_filter_percentile}`), matching the pileup step. Entropy has no `--seed`/`--region`, but its input is the MT-only BAM whose read count is below modkit's default sampling size, so the estimate is MT-restricted and effectively deterministic (uses all reads). `--min-coverage ${params.modkit_min_coverage}` (=10, raised from the modkit default of 3 because mtDNA depth is ~100-800x). ### modkit_localize @@ -307,6 +645,8 @@ Outputs are both optional. `${id}_localize.tsv` (emit `localize`) and `*.log` (e The script first runs `samtools faidx` and derives a genome-sizes TSV (`\t`) from the resulting `.fai` via `cut -f1,2`, because `modkit localize` requires that file. The `--window` is set to 200 bp instead of the modkit default of 2000 bp: the whole mtDNA is only ~16.5 kb, so a small window around each regulatory-feature midpoint is appropriate. The window is tunable. Output is per-sample; aggregation across samples is done downstream in R. +`--min-coverage ${params.modkit_min_coverage}` (10) is passed as well, raised from the modkit default of 3 because mtDNA depth is about 100-800x; the mouse parameter files lower it to 3. + Load-bearing gotchas: - The `.tbi` is taken as an explicit `path(tbi)` input on purpose. `modkit localize` expects the tabix index next to the input as `.tbi`, so both files must land in the same task work dir. Dropping the explicit input stops Nextflow from staging the index and breaks the step. - `errorStrategy 'ignore'` is intentional. This is an exploratory per-sample fan-out where a low-coverage sample can legitimately exit non-zero; ignoring keeps one weak sample from aborting the whole run, and the failure is still recorded in the Nextflow execution report. @@ -319,6 +659,8 @@ Process `MODKIT_MOTIF_SEARCH` runs `modkit motif search` per sample (fan-out). I The search runs over all modification codes present (6mA/5mC/5hmC) because no `--mod-code` is passed; modkit then processes every code rather than hiding 5mC/5hmC motifs. `--known-motif GATC 1 a` is passed as a hypothesis so the 6mA frequency at GATC is always reported without imposing GATC as a filter. The published GATC motif came from a DpnI assay that only sees GATC (instrument bias); Nanopore is context-agnostic, so motifs are discovered de novo and GATC is merely checked. +`--min-coverage ${params.modkit_min_coverage}` (10) is passed as well, raised from the modkit default of 3 because mtDNA depth is about 100-800x; the mouse parameter files lower it to 3. + Gotchas, in order of how easy they are to reintroduce as bugs: - The `.tbi` is taken as an explicit `path(tbi)` input so Nextflow stages it alongside the `.bed.gz` in the work directory. `--contig` requires bgzip+tabix input, and modkit assumes the index lives next to the input as `.tbi`. Dropping the explicit `path(tbi)` input breaks `--contig` at runtime even though the pipeline graph looks fine. @@ -374,3 +716,61 @@ MERGE_HAPLOGROUPS joins each per-sample rank-1 row to the samplesheet by real_ba Load-bearing gotchas: - HAPLOGROUP always exits 0 and writes a published log. `publishDir` does not run for a failed task, so under `errorStrategy 'ignore'` a crashed task left no trace, which is how haplogroup originally errored. The task now always exits 0 and publishes its log; the `.txt` exists only on success, so a bad sample neither sinks step 3 nor hides its cause. - MERGE strips Haplogrep3's quoted fields first (`gsub(/"/, "")`) and keeps only rank-1 rows (`$3 == "1"`). It re-splits the samplesheet on `,` because the per-sample TSVs are tab-separated. +- haplogrep3 has no `--version` flag (invoking it errors with exit 2). The HAPLOGROUP script therefore runs under `set -uo pipefail` without `-e` and deliberately contains no `--version` probe: an earlier version probed `--version` under `set -e`, which aborted the script before `classify` ran, and because the failed task published nothing this was the silent run01 F6 haplogroup failure. Do not reintroduce `set -e` or a `--version` check in this block. + +### heteroplasmy + +File: src/modules/heteroplasmy.nf + +Feature: mtDNA heteroplasmy calling, optional and gated by `params.enable_heteroplasmy` (default `false`). This is a *separate* caller from the `haplogroup` Clair3 step and serves a different goal: quantifying low-level intra-individual variation (instability) rather than the haploid consensus. Clair3 runs in `--haploid_precise`, the mode designed to *suppress* low-AF events, so its "heteroplasmy" is carry-over and unfit for a group claim; mutserve2 (Bayesian model: 1000G prior + base-quality likelihood) is the purpose-built low-frequency caller. The full scientific design (cohort coverage confounder, why coverage-matching, the 3/5/10% floors, the N reconciliation) lives in `DESIGN_heteroplasmia_caller_2026-06-14.md`; this section documents the code only. + +Key premise that simplifies everything: the BAMs are already aligned to whole-genome GRCh38 (see `#nextflow` NUMT control, including its near-identical-copy caveat), so the NuMT defense is already in place. The caller therefore only needs to *extract* chrM from the existing filtered BAMs — there is no re-alignment and no need for POD5/FASTQ. + +Integration choice: rather than running the genepi/mtdna-server-2 Nextflow pipeline nested inside ours (its `min_mean_coverage=50` QC gate would reject our 40x-matched samples), we run its components directly from the official container `params.mtdnaserver_sif` (Apptainer of `quay.io/genepi/mtdna-server-2:v2.1.16`, which ships mutserve 2.0.3, haplocheck 1.3.3, samtools/bcftools 1.19, java). The `MUTSERVE`, `MERGE_VCFS_HET` and `HAPLOCHECK_HET` processes use that image; `EXTRACT_MT`, `DOWNSAMPLE_MATCH`, `PARSE_HETEROPLASMY` and `MERGE_HETEROPLASMY` use the main image (samtools/awk). Every process is `errorStrategy 'ignore'` because it is a per-sample fan-out that must not sink step 3. + +Flow (wired in `MODKIT_AND_MULTIQC.nf` under `if (params.enable_heteroplasmy)`): EXTRACT_MT now runs ONCE in step 3 core (shared with the methylation modkit steps — see #extract_mt), and heteroplasmy reuses its `mt_bam` output (no second extraction) -> two arms. The **native arm** runs MUTSERVE on the full-depth mito BAM (sensitivity check + the input to Haplocheck). The **matched arm** runs DOWNSAMPLE_MATCH then MUTSERVE on the coverage-matched BAM (the primary outcome). Both arms feed PARSE_HETEROPLASMY -> MERGE_HETEROPLASMY. MUTSERVE is included twice via `MUTSERVE as MUTSERVE_NATIVE` / `as MUTSERVE_MATCHED`; the `arm` value ('native'/'matched') drives both the output filenames and the `mutserve//` publish subdir. + +EXTRACT_MT (now in `../modules/extract_mt.nf`, shared — see #extract_mt) pulls only `${params.mt_contig}`, sets `@RG ID/SM:`, sorts and indexes. Supplementary alignments are kept by default (`mt_keep_supplementary=true`): reads spanning the circular 16569->1 origin become primary+supplementary, and dropping them would lose D-loop/origin coverage. **IMPORTANT — this MT BAM is now shared with the methylation steps**, so `mt_keep_supplementary` is a GLOBAL switch: `false` adds `-F 0x800` at extraction and therefore removes supplementary from the modkit pileup/extract/entropy too (degrading F7-A D-loop recovery via `--allow-non-primary`), not just from the caller. Treat it as an intentional reproducibility tradeoff. mutserve includes supplementary by default and has no flag to exclude them (verified in `BamAnalyser.java`), so this extraction is the only control point. + +DOWNSAMPLE_MATCH equalizes depth at the source rather than adjusting for it later (count data; a covariate does not fix detection bias). It mirrors genepi's `subsampling.nf`: raw meandepth from `samtools coverage` field 7, fraction = `T/mean`, then `samtools view --subsample-seed --subsample `. `T = het_downsample_target` (40x, cohort-derived in the DESIGN doc, not guessed). Samples below T are dropped from the matched arm (strict matching) but survive in the native arm; every sample writes a `*.downsample.log` recording mean/T/fraction/status. + +`T` is derived in section 2.2 of `DESIGN_heteroplasmia_caller_2026-06-14.md`, and `het_downsample_seed` (42) is a fixed `samtools view --subsample-seed` so that the coverage-matched subsample is reproducible run to run. Inclusion in the matched arm is decided by a floating-point comparison in awk (`mean >= T` and `mean > 0`), not by an integer `-gt` on a rounded value, so a sample sitting exactly at the target depth is kept and downsampled with fraction 1.0 (a no-op); because only samples at or above `T` pass, the fraction `T/mean` is always at most 1.0. A sample whose mean depth resolves to 0 (an empty MT BAM, or an `mt_contig` mismatch that makes `samtools coverage` print nothing, in which case `mean` falls back to 0 via `${mean:-0}`) is excluded from the matched arm exactly like a below-target sample and logged with `status=excluded_below_target`; it still runs in the native arm. + +MUTSERVE mirrors genepi's `modules/local/mutserve.nf` verbatim: `java -jar /opt/mutserve/mutserve.jar call --level ${het_detection_level} --reference --mapQ ${het_mapq} --baseQ ${het_baseq} --no-ansi --strand-bias ${het_strand_bias} --write-raw`, then `bcftools norm -m-any` to split multiallelics and `tabix`. `--level 0.01` captures everything >=1%; the 3/5/10% floors are applied downstream in R, because `--level` is a reporting threshold (not a model change), so filtering by AF in R is equivalent and lets all three floors come from one call. The reference is `reference_file`, which in step 3 is the small mito `MT.fa` (GRCh38 chrM = rCRS), the same file modkit uses. + +`het_strand_bias` is 1.6 because that is the value the genepi mtDNA-Server 2 pipeline passes; the mutserve tool default is 1.2, so leaving the flag off would silently change the strand-bias filter. The `bcftools norm -m-any` after the call splits multiallelic records so that each alternate allele becomes its own row, which keeps per-site counting consistent with the genepi pipeline. + +PARSE_HETEROPLASMY parses the mutserve **variant table** (`${id}_${arm}.txt`, written automatically beside the VCF) rather than the VCF: with `awk` it keeps `Filter==PASS` rows and emits per `(sample, arm)` a TSV of `pos/ref/variant/variant_level/coverage/cov_fwd/cov_rev/type`. The table is preferred over `bcftools query [%AF]` for two reasons confirmed at the smoke-test: it carries per-strand coverage (the §5e strand control) plus an explicit `Type` (1=homoplasmy, 2=heteroplasmy/low-level) and `Filter` (drops `BLACKLISTED` sites), and it sidesteps a bgzip gotcha — mutserve writes the `.vcf.gz` as plain gzip, which `bcftools query` cannot read until the `norm -Oz` re-block. It deliberately does NOT classify or apply floors — that stays in R. MERGE_HETEROPLASMY concatenates all per-sample TSVs and joins the samplesheet by real_barcode (same two-pass awk as MERGE_MTDNA_METRICS) into `heteroplasmy_variants_master.tsv` (columns: sample_id, arm, real_barcode, run_id, group, pos, ref, variant, variant_level, coverage, cov_fwd, cov_rev, type). + +Column positions in the mutserve table consumed by the awk parser are 1=ID, 2=Filter, 3=Pos, 4=Ref, 5=Variant, 6=VariantLevel, 11=Coverage, 12=CoverageFWD, 13=CoverageREV and 14=Type; the parser emits `id`, `arm` and then columns 3, 4, 5, 6, 11, 12, 13, 14 in that order. MERGE_HETEROPLASMY reads `run_id` from samplesheet column 1, `real_barcode` from column 3 and `group` from column 7 (the same 1/3/7 layout `MERGE_HAPLOGROUPS` uses); the join key is the `sample_id` prefix before the first `-`, and a barcode absent from the samplesheet gets `run_id` and `group` = `NA` rather than dropping the row. + +Contamination QC: when `params.enable_haplocheck` is set, MERGE_VCFS_HET `bcftools merge`s the native-arm VCFs into one multi-sample VCF and HAPLOCHECK_HET runs `haplocheck.jar --out haplocheck.txt --raw` on it (mirrors genepi `haplogroups_contamination.nf`). Haplocheck detects *sample* contamination via phylogeny (a different problem from NuMT); it is used here as a QC gate, because contamination inflates apparent heteroplasmy and could be a differential confounder. It runs on the native arm because more depth gives a better contamination signal. + +The distinct sample columns of the merged VCF come from the per-sample `@RG SM` set at extraction (`#extract_mt`). `bcftools merge` requires at least two input VCFs, so when only one native-arm VCF survives, MERGE_VCFS_HET copies it through as `heteroplasmy_native_merged.vcf.gz` instead of merging (Haplocheck accepts a single-sample VCF). In `MODKIT_AND_MULTIQC.nf` the native-arm VCF and TBI collects feeding MERGE_VCFS_HET are deliberately bare `.collect()` with no `.ifEmpty([])`: `bcftools merge` rejects empty input, so if every native mutserve call fails the Haplocheck QC is skipped (the process never fires) rather than crashing on an empty merge; a run in that state is already degraded, so skipping the contamination check is acceptable. + +Load-bearing gotchas: +- The extracted BAM keeps the whole-genome `@SQ` header. Do NOT "clean" it by deleting non-MT `@SQ` lines: `samtools reheader` only rewrites header text and does not remap read `tid`s, so dropping `@SQ` lines corrupts the BAM. This is why `samtools coverage` is called with `-r ${params.mt_contig}` (otherwise `NR==2` would be the first whole-genome contig, not MT). +- `het_mapq=20` (mutserve default) is applied on top of the step-2 MAPQ>=10 filter, so the effective MAPQ gate is 20. That is intentional (more conservative for low-frequency calls). +- `--level` is fixed at 0.01 explicitly because the genepi nextflow.config default (0.01) and the web UI default (0.02) diverge. + +Validation status after the b41 (run07) smoke-test (`slurm_jobs/05_smoke_heteroplasmy.sbatch`): (1) CONFIRMED — bundled mutserve 2.0.3 accepts `--strand-bias`, `--write-raw`, `--no-ansi`, `--contig-name`, `--level`, etc. (the public README is stale; the JAR `--help` lists them). (2) CONFIRMED — mutserve accepts `MT.fa` (header `>MT`) as `--reference` against the `MT`-contig BAM with no `--contig-name` (call exit 0; heteroplasmy detected at pos 204/2814/11762/16189). (3) CONFIRMED — AF is in the VCF FORMAT (`GT:AF:BQ:DP`); now moot since the parser reads the `.txt`. Gotcha recorded: mutserve writes `.vcf.gz` as plain gzip, so any `bcftools` use needs the `norm -Oz` re-block first (MUTSERVE does it for the Haplocheck VCF). The `${id}_${arm}.txt` and `*_raw.txt` output names are CONFIRMED (b41 produced `b41.txt` + `b41_raw.txt`). (4) STILL TO CONFIRM at the pilot run: `samtools addreplacerg -r 'ID:..' -r 'SM:..'` (field form) on the container's samtools — the smoke extracted MT without it. + +### mt_variant_scan + +Opt-in arm (`enable_mt_variant_scan`, default off) that runs the tools described under +[numt_rescue](#numt_rescue) (`mt_allele_counts.sh`, `mt_spectrum_phasing.R`, `mt_read_bases.sh`) +on the shared MT-only BAM of step 3: per-sample allele counts, one cohort call of sites, spectrum +and group tests, per-read bases at each sample's called sites, and one cohort phasing pass. +Off by default because the thesis dropped heteroplasmy as an axis, and because at the depths of +these projects (blood 96 to 208x, mice 13 to 67x) the arm is a sample-integrity check (mixtures, +contamination, strain identity) rather than a group comparison; see the first-run lessons in the +numt_rescue section. The R stages run in a separate container (`mt_variant_r_sif`, R with +data.table and ggplot2; `mt_variant_r_libs` sets `R_LIBS_USER` when the packages live outside the +image) because the main image has no R. Each per-sample process links the staged BAM under the +`-...realigned.bam` name the tools glob, so the tools keep their standalone interface. The +cohort here is the run being processed (`run_id`), so group tests are skipped when a run holds +fewer than four samples or one group, and phasing is skipped for samples with fewer than two +called sites; the cohort-wide pass over all human runs at once was done with the standalone +launchers (`~/slurm_jobs/nanopore/variants_pipeline.sh`). + +The arm runs on the shared MT-only `mt_bam` of step 3 rather than extracting again. `MT_VARIANT_CALL` writes one `.sites.bed` per sample, and those site BEDs are re-keyed to the sample id by stripping the `.sites.bed` suffix before being joined back to `mt_bam` for `MT_READ_BASES`, so the per-read pass is keyed by real_barcode like everything else in step 3 and never by emission order. diff --git a/parameters/human_blood/modkit.yaml b/parameters/human_blood/modkit.yaml index 6913f1ed08dde2a9fbc8d3293efda41bb0bd6554..44fecc041103c7ff5eba8e2472dccd2b2eb00733 100644 --- a/parameters/human_blood/modkit.yaml +++ b/parameters/human_blood/modkit.yaml @@ -5,4 +5,8 @@ modkit_threads: 8 modkit_isize: 10000 modkit_extract_qsize: 1000 samtools_threads: 12 -multiqc_config: "./references/multiqc_config.yaml" \ No newline at end of file +multiqc_config: "./references/multiqc_config.yaml" +rescue_mt: true +modkit_filter_mode: "fixed" +modkit_filter_threshold: "0.75" +enable_mt_variant_scan: false diff --git a/parameters/human_blood/qc.yaml b/parameters/human_blood/qc.yaml index d7cfb6d8d6d3d9519c7539cd24a93a65ea0cfaa7..4322ffca775a913004f10264c15cd14bb26b2b24 100644 --- a/parameters/human_blood/qc.yaml +++ b/parameters/human_blood/qc.yaml @@ -4,4 +4,6 @@ run_id: "run01" mapq: "10" qscore_thresh: "9" min_mapped_reads_thresh: 500 -samtools_threads: 12 \ No newline at end of file +samtools_threads: 12 +rescue_mt: true +numt_bed: "/home/joaochrusciel/bio_projects/nanopore/references/numt/grch38/numt_loci_grch38.bed" diff --git a/parameters/letizia_mouse/basecall.yaml b/parameters/letizia_mouse/basecall.yaml new file mode 100644 index 0000000000000000000000000000000000000000..c2c49672f42259f1f970beb6d60fdb878b2659fd --- /dev/null +++ b/parameters/letizia_mouse/basecall.yaml @@ -0,0 +1,10 @@ +project_name: "letizia_mouse" +step: 1 +run_id: "letizia_all" +reference_file: "/home/joaochrusciel/bio_projects/nanopore/references/mice.GRC39m.genome.fa" +mt_contig: "chrM" +basecall_speed: "hac@v6.0.0" +basecall_mods: "5mC_5hmC,6mA" +barcoding_kit: "SQK-RBK114-24" +samplesheet: "/home/joaochrusciel/bio_projects/nanopore/references/samplesheet_letizia.csv" +fast5_batch_size: 25 diff --git a/parameters/letizia_mouse/modkit.yaml b/parameters/letizia_mouse/modkit.yaml new file mode 100644 index 0000000000000000000000000000000000000000..75a2f292145dc0e2f41b6b46d1acd96a9be0433b --- /dev/null +++ b/parameters/letizia_mouse/modkit.yaml @@ -0,0 +1,20 @@ +project_name: "letizia_mouse" +step: 3 +run_id: "letizia_all" +reference_file: "/home/joaochrusciel/bio_projects/nanopore/references/mice.GRC39m.genome.fa" +mt_contig: "chrM" +samplesheet: "/home/joaochrusciel/bio_projects/nanopore/references/samplesheet_letizia.csv" +enable_haplogroup: false +enable_haplocheck: false +enable_heteroplasmy: false +mt_macro_bed: "/home/joaochrusciel/bio_projects/nanopore/references/mt_macro_mouse.bed" +mt_genes_bed: "/home/joaochrusciel/bio_projects/nanopore/references/mt_genes_mouse.bed" +mt_control_elements_bed: "/home/joaochrusciel/bio_projects/nanopore/references/mt_control_elements_mouse.bed" +mt_mdp_bed: "/home/joaochrusciel/bio_projects/nanopore/references/mt_mdp_mouse.bed" +modkit_min_coverage: 3 +modkit_threads: 8 +modkit_isize: 10000 +modkit_extract_qsize: 1000 +multiqc_config: "./references/multiqc_config.yaml" +rescue_mt: true +enable_mt_variant_scan: false diff --git a/parameters/letizia_mouse/qc.yaml b/parameters/letizia_mouse/qc.yaml new file mode 100644 index 0000000000000000000000000000000000000000..5f8e9dbd484bba790c90363a275d28f5f7c3ff36 --- /dev/null +++ b/parameters/letizia_mouse/qc.yaml @@ -0,0 +1,10 @@ +project_name: "letizia_mouse" +step: "2_from_step_1" +run_id: "letizia_all" +mt_contig: "chrM" +mapq: "10" +qscore_thresh: "9" +min_mapped_reads_thresh: 150 +samtools_threads: 12 +rescue_mt: true +numt_bed: "/home/joaochrusciel/bio_projects/nanopore/references/numt/grcm39/numt_loci_grcm39.bed" diff --git a/parameters/nicotine_hippo_mm10/basecall.yaml b/parameters/nicotine_hippo_mm10/basecall.yaml new file mode 100644 index 0000000000000000000000000000000000000000..43f67268e630ddea60d0b89b99e2369619fcdc10 --- /dev/null +++ b/parameters/nicotine_hippo_mm10/basecall.yaml @@ -0,0 +1,10 @@ +project_name: "nicotine_hippo_mm10" +step: 1 +run_id: "nic_hippo_run2" +reference_file: "/home/joaochrusciel/bio_projects/nanopore/references/mm10/mm10.fa" +mt_contig: "chrM" +basecall_speed: "hac@v6.0.0" +basecall_mods: "5mC_5hmC,6mA" +barcoding_kit: "SQK-RBK114-24" +samplesheet: "/home/joaochrusciel/bio_projects/nanopore/references/samplesheet_hippo_mm10.csv" +fast5_batch_size: 25 diff --git a/parameters/nicotine_hippo_mm10/modkit.yaml b/parameters/nicotine_hippo_mm10/modkit.yaml new file mode 100644 index 0000000000000000000000000000000000000000..abe3194a2a1965af8f1d16bf7dbdede745dbbdc4 --- /dev/null +++ b/parameters/nicotine_hippo_mm10/modkit.yaml @@ -0,0 +1,20 @@ +project_name: "nicotine_hippo_mm10" +step: 3 +run_id: "nic_hippo_run2" +reference_file: "/home/joaochrusciel/bio_projects/nanopore/references/mm10/mm10.fa" +mt_contig: "chrM" +methyl_scope: "genome" +samplesheet: "/home/joaochrusciel/bio_projects/nanopore/references/samplesheet_hippo_mm10.csv" +enable_haplogroup: false +enable_haplocheck: false +enable_heteroplasmy: false +mt_macro_bed: "/home/joaochrusciel/bio_projects/nanopore/references/mt_macro_mouse.bed" +mt_genes_bed: "/home/joaochrusciel/bio_projects/nanopore/references/mt_genes_mouse.bed" +mt_control_elements_bed: "/home/joaochrusciel/bio_projects/nanopore/references/mt_control_elements_mouse.bed" +mt_mdp_bed: "/home/joaochrusciel/bio_projects/nanopore/references/mt_mdp_mouse.bed" +modkit_min_coverage: 3 +modkit_threads: 8 +modkit_isize: 10000 +modkit_extract_qsize: 1000 +multiqc_config: "./references/multiqc_config.yaml" +enable_mt_variant_scan: false diff --git a/parameters/nicotine_hippo_mm10/qc.yaml b/parameters/nicotine_hippo_mm10/qc.yaml new file mode 100644 index 0000000000000000000000000000000000000000..db0dc5105ce3705b34e06fa2ec694c38c0c2aebb --- /dev/null +++ b/parameters/nicotine_hippo_mm10/qc.yaml @@ -0,0 +1,9 @@ +project_name: "nicotine_hippo_mm10" +step: "2_from_step_1" +run_id: "nic_hippo_run2" +reference_file: "/home/joaochrusciel/bio_projects/nanopore/references/mm10/mm10.fa" +mt_contig: "chrM" +mapq: "10" +qscore_thresh: "9" +min_mapped_reads_thresh: 150 +samtools_threads: 12 diff --git a/references/numt/grch38/mt_contig_md5.tsv b/references/numt/grch38/mt_contig_md5.tsv new file mode 100644 index 0000000000000000000000000000000000000000..15083d6a23f3bef1a063762525a677499efc8a3f --- /dev/null +++ b/references/numt/grch38/mt_contig_md5.tsv @@ -0,0 +1,4 @@ +source length md5_sequence path +genome:MT 16569 c68f52674c9fb33aef52dcf399755519 /home/joaochrusciel/bio_projects/nanopore/references/numt/grch38/MT.fa +MT.fa 16569 c68f52674c9fb33aef52dcf399755519 /home/joaochrusciel/bio_projects/nanopore/references/MT.fa +chrMT.fa 16569 c68f52674c9fb33aef52dcf399755519 /home/joaochrusciel/bio_projects/nanopore/references/chrMT.fa diff --git a/references/numt/grch38/numt_core_grch38.bed b/references/numt/grch38/numt_core_grch38.bed new file mode 100644 index 0000000000000000000000000000000000000000..63112978d61def2d45d4c5e46f2d89cec2b3dc13 --- /dev/null +++ b/references/numt/grch38/numt_core_grch38.bed @@ -0,0 +1,23 @@ +1 629083 634924 +11 10507886 10510336 +13 109424124 109424380 +14 32484097 32485118 +17 22519207 22520853 +17 22522098 22531468 +17 53105732 53106385 +2 49229632 49229899 +2 87824889 87825350 +3 96617187 96618510 +4 12640293 12640635 +4 155463824 155466145 +5 80650021 80652368 +5 94567455 94570918 +5 100045937 100055045 +5 123760825 123761738 +5 134923308 134928527 +6 61574112 61574629 +7 142665257 142667693 +9 33656687 33659130 +X 55179796 55183964 +X 126471703 126472452 +X 126472466 126473284 diff --git a/references/numt/grch38/numt_hits_grch38.tsv b/references/numt/grch38/numt_hits_grch38.tsv new file mode 100644 index 0000000000000000000000000000000000000000..ae62d816bd87c7e9de28395f5e2fff8921d5980a --- /dev/null +++ b/references/numt/grch38/numt_hits_grch38.tsv @@ -0,0 +1,25 @@ +chrom start end mt_contig mt_start mt_end strand matches aln_len mapq identity type +1 629083 634924 MT 3913 9755 + 5758 5843 60 0.9855 P +11 10507886 10510336 MT 520 2972 - 2310 2452 60 0.9421 P +13 109424124 109424380 MT 981 1237 - 254 256 60 0.9922 P +14 32484097 32485118 MT 5582 6606 - 953 1024 60 0.9307 P +17 22519207 22520853 MT 14377 16024 + 1361 1648 60 0.8258 P +17 22522098 22531468 MT 636 10064 + 7932 9441 60 0.8402 P +17 53105732 53106385 MT 6817 7470 + 628 653 0 0.9617 S +2 49229632 49229899 MT 6740 7007 - 254 267 0 0.9513 S +2 87824889 87825350 MT 8313 8774 + 442 461 0 0.9588 S +3 96617187 96618510 MT 1395 2718 - 1262 1323 0 0.9539 S +4 12640293 12640635 MT 9337 9679 - 321 342 0 0.9386 S +4 155463824 155466145 MT 989 3305 - 1913 2342 60 0.8168 P +5 80650021 80652368 MT 340 2697 - 2223 2358 60 0.9427 P +5 94567455 94570918 MT 12661 16124 + 3024 3465 0 0.8727 S +5 100045937 100055045 MT 6116 15183 - 8057 9118 1 0.8836 P +5 123760825 123761738 MT 597 1510 + 782 920 0 0.8500 S +5 134923308 134928527 MT 10268 15487 - 4909 5219 0 0.9406 S +6 61574112 61574629 MT 2417 2934 + 471 517 0 0.9110 S +7 142665257 142667693 MT 654 3093 + 2060 2455 0 0.8391 S +9 33656687 33659130 MT 650 3093 + 2067 2459 1 0.8406 P +X 55179796 55183964 MT 638 4831 - 3361 4225 0 0.7955 S +X 126471703 126472452 MT 6552 7302 + 698 750 0 0.9307 S +X 126472466 126472732 MT 686 953 - 258 267 0 0.9663 S +X 126472730 126473284 MT 10605 11159 + 519 554 0 0.9368 S diff --git a/references/numt/grch38/numt_loci_grch38.bed b/references/numt/grch38/numt_loci_grch38.bed new file mode 100644 index 0000000000000000000000000000000000000000..340d9b75a190403331018d8e37439b63a05a4bf4 --- /dev/null +++ b/references/numt/grch38/numt_loci_grch38.bed @@ -0,0 +1,21 @@ +1 624083 639924 +11 10502886 10515336 +13 109419124 109429380 +14 32479097 32490118 +17 22514207 22536468 +17 53100732 53111385 +2 49224632 49234899 +2 87819889 87830350 +3 96612187 96623510 +4 12635293 12645635 +4 155458824 155471145 +5 80645021 80657368 +5 94562455 94575918 +5 100040937 100060045 +5 123755825 123766738 +5 134918308 134933527 +6 61569112 61579629 +7 142660257 142672693 +9 33651687 33664130 +X 55174796 55188964 +X 126466703 126478284 diff --git a/references/numt/grch38/numt_near_identical_grch38.tsv b/references/numt/grch38/numt_near_identical_grch38.tsv new file mode 100644 index 0000000000000000000000000000000000000000..b1e79781fd7bbe968adf230302791e2b3d02ce89 --- /dev/null +++ b/references/numt/grch38/numt_near_identical_grch38.tsv @@ -0,0 +1,2 @@ +chrom start end mt_contig mt_start mt_end strand matches aln_len mapq identity type +1 629083 634924 MT 3913 9755 + 5758 5843 60 0.9855 P diff --git a/references/numt/grcm39/mt_contig_md5.tsv b/references/numt/grcm39/mt_contig_md5.tsv new file mode 100644 index 0000000000000000000000000000000000000000..1c71fdf8eff54cb3b85a06ee466b833c43fb33c7 --- /dev/null +++ b/references/numt/grcm39/mt_contig_md5.tsv @@ -0,0 +1,4 @@ +source length md5_sequence path +genome:chrM 16299 11c8af2a2528b25f2c080ab7da42edda /home/joaochrusciel/bio_projects/nanopore/references/numt/grcm39/chrM.fa +chrM.fa 16299 11c8af2a2528b25f2c080ab7da42edda /home/joaochrusciel/bio_projects/nanopore/references/chrM.fa +chrM.fa 16299 11c8af2a2528b25f2c080ab7da42edda /home/joaochrusciel/bio_projects/nanopore/references/mm10/chrM.fa diff --git a/references/numt/grcm39/numt_core_grcm39.bed b/references/numt/grcm39/numt_core_grcm39.bed new file mode 100644 index 0000000000000000000000000000000000000000..ff40e8e289da8df61935ef09f429c31ca486db76 --- /dev/null +++ b/references/numt/grcm39/numt_core_grcm39.bed @@ -0,0 +1,11 @@ +chr1 24650615 24655265 +chr1 132718643 132719096 +chr10 84861831 84862504 +chr10 95913386 95913756 +chr11 90429307 90429963 +chr12 97028201 97028566 +chr13 85274680 85275634 +chr2 22477299 22480546 +chr4 79920566 79923437 +chr5 60200060 60200442 +chr6 9889892 9890197 diff --git a/references/numt/grcm39/numt_hits_grcm39.tsv b/references/numt/grcm39/numt_hits_grcm39.tsv new file mode 100644 index 0000000000000000000000000000000000000000..419c32f828d7723c5140bba4f119819c3270f315 --- /dev/null +++ b/references/numt/grcm39/numt_hits_grcm39.tsv @@ -0,0 +1,12 @@ +chrom start end mt_contig mt_start mt_end strand matches aln_len mapq identity type +chr1 24650615 24655265 chrM 6393 11042 - 4648 4650 60 0.9996 P +chr1 132718643 132719096 chrM 15600 16054 + 416 454 37 0.9163 P +chr10 84861831 84862504 chrM 15320 15991 - 581 682 60 0.8519 P +chr10 95913386 95913756 chrM 1511 1881 - 362 370 60 0.9784 P +chr11 90429307 90429963 chrM 6215 6871 - 595 656 60 0.9070 P +chr12 97028201 97028566 chrM 15674 16039 + 364 365 60 0.9973 P +chr13 85274680 85275634 chrM 12445 13400 - 909 955 60 0.9518 P +chr2 22477299 22480546 chrM 4440 7699 - 3174 3260 60 0.9736 P +chr4 79920566 79923437 chrM 12487 15356 + 2671 2871 60 0.9303 P +chr5 60200060 60200442 chrM 4849 5231 + 375 383 60 0.9791 P +chr6 9889892 9890197 chrM 15611 15916 - 292 305 60 0.9574 P diff --git a/references/numt/grcm39/numt_loci_grcm39.bed b/references/numt/grcm39/numt_loci_grcm39.bed new file mode 100644 index 0000000000000000000000000000000000000000..9d204b25f81910600583094771a3f75ff79e66f5 --- /dev/null +++ b/references/numt/grcm39/numt_loci_grcm39.bed @@ -0,0 +1,11 @@ +chr1 24645615 24660265 +chr1 132713643 132724096 +chr10 84856831 84867504 +chr10 95908386 95918756 +chr11 90424307 90434963 +chr12 97023201 97033566 +chr13 85269680 85280634 +chr2 22472299 22485546 +chr4 79915566 79928437 +chr5 60195060 60205442 +chr6 9884892 9895197 diff --git a/references/numt/grcm39/numt_near_identical_grcm39.tsv b/references/numt/grcm39/numt_near_identical_grcm39.tsv new file mode 100644 index 0000000000000000000000000000000000000000..204788064ed86ac235cc21b069a3bc555d4ed0e5 --- /dev/null +++ b/references/numt/grcm39/numt_near_identical_grcm39.tsv @@ -0,0 +1,2 @@ +chrom start end mt_contig mt_start mt_end strand matches aln_len mapq identity type +chr1 24650615 24655265 chrM 6393 11042 - 4648 4650 60 0.9996 P diff --git a/src/bin/derive_numt_loci.sh b/src/bin/derive_numt_loci.sh new file mode 100755 index 0000000000000000000000000000000000000000..66e3cafa2e09f48f671f7cbc763244b858d7e5a9 --- /dev/null +++ b/src/bin/derive_numt_loci.sh @@ -0,0 +1,81 @@ +#!/usr/bin/env bash +#' Derive the nuclear loci homologous to the mitochondrial contig (NUMTs) for one genome build. +#' +#' The mitochondrial contig is extracted from the genome itself and aligned back against the +#' whole genome with minimap2 in assembly mode, every hit of at least MIN_LEN aligned bases is +#' kept, and three files are written: a hit table with identity and the mitochondrial interval +#' each copy corresponds to, a flank-padded merged bed (the candidate regions read by +#' rescue_mt.sh) and a core bed without flanks. The mitochondrial sequence md5 is compared with +#' every fasta listed in MT_FA. Runs inside images/debian-nanopore.sif. +#' Rationale: docs/pipeline_notes.md#numt-rescue +#' +#' Environment: GENOME MT_CONTIG BUILD OUT_DIR [MT_FA colon-separated fasta paths to compare] +#' [MIN_LEN=200] [FLANK=5000] [THREADS=8] [BATCH=500M] [NEAR_IDENT=0.98] [NEAR_LEN=1000] +set -euo pipefail +export LC_ALL=C +: "${GENOME:?}" "${MT_CONTIG:?}" "${BUILD:?}" "${OUT_DIR:?}" +MT_FA="${MT_FA:-}" +MIN_LEN="${MIN_LEN:-200}" +FLANK="${FLANK:-5000}" +THREADS="${THREADS:-8}" +BATCH="${BATCH:-500M}" +NEAR_IDENT="${NEAR_IDENT:-0.98}" +NEAR_LEN="${NEAR_LEN:-1000}" +mkdir -p "$OUT_DIR" +seqmd5() { grep -v "^>" "$1" | tr -d "\n\r" | tr "a-z" "A-Z" | md5sum | cut -d " " -f1; } +seqlen() { grep -v "^>" "$1" | tr -d "\n\r" | wc -c; } + +echo "[info] build=$BUILD genome=$GENOME contig=$MT_CONTIG threads=$THREADS min_len=$MIN_LEN flank=$FLANK" +if [ ! -s "$GENOME.fai" ]; then echo "[info] indexing genome"; samtools faidx "$GENOME"; fi +if ! awk -v c="$MT_CONTIG" '$1==c{f=1} END{exit !f}' "$GENOME.fai"; then + echo "[error] contig $MT_CONTIG not in $GENOME.fai; first contigs:"; cut -f1 "$GENOME.fai" | head -30; exit 1 +fi +MT_OUT="$OUT_DIR/$MT_CONTIG.fa" +samtools faidx "$GENOME" "$MT_CONTIG" > "$MT_OUT" +samtools faidx "$MT_OUT" +MT_LEN=$(cut -f2 "$MT_OUT.fai") +MD5="$OUT_DIR/mt_contig_md5.tsv" +printf 'source\tlength\tmd5_sequence\tpath\n' > "$MD5" +printf '%s\t%s\t%s\t%s\n' "genome:$MT_CONTIG" "$MT_LEN" "$(seqmd5 "$MT_OUT")" "$MT_OUT" >> "$MD5" +for f in $(echo "$MT_FA" | tr ":" " "); do + if [ -s "$f" ]; then printf '%s\t%s\t%s\t%s\n' "$(basename "$f")" "$(seqlen "$f")" "$(seqmd5 "$f")" "$f" >> "$MD5"; else echo "[warn] $f missing"; fi +done +echo "[info] mitochondrial sequence per fasta:"; cat "$MD5" +n_md5=$(tail -n +2 "$MD5" | cut -f3 | sort -u | wc -l) +if [ "$n_md5" -eq 1 ]; then echo "[ok] every mitochondrial fasta carries the same sequence"; else echo "[warn] mitochondrial fasta copies differ; the rescue must use $MT_OUT"; fi + +PAF="$OUT_DIR/numt_raw_$BUILD.paf" +echo "[info] minimap2 -cx asm20 -I $BATCH -N 500 -p 0.01 --secondary=yes genome $MT_CONTIG" +minimap2 -t "$THREADS" -cx asm20 -I "$BATCH" -N 500 -p 0.01 --secondary=yes "$GENOME" "$MT_OUT" > "$PAF" 2> "$OUT_DIR/minimap2_$BUILD.log" +echo "[info] raw alignments: $(wc -l < "$PAF")" + +HITS="$OUT_DIR/numt_hits_$BUILD.tsv" +printf 'chrom\tstart\tend\tmt_contig\tmt_start\tmt_end\tstrand\tmatches\taln_len\tmapq\tidentity\ttype\n' > "$HITS" +awk -F"\t" -v mt="$MT_CONTIG" -v minlen="$MIN_LEN" 'BEGIN{OFS="\t"} $6!=mt && $11>=minlen {tp="NA"; for(i=13;i<=NF;i++) if($i ~ /^tp:A:/) tp=substr($i,6,1); print $6,$8,$9,$1,$3,$4,$5,$10,$11,$12,sprintf("%.4f",$10/$11),tp}' "$PAF" \ + | sort -t "$(printf '\t')" -k1,1 -k2,2n >> "$HITS" + +mkbed() { + awk -F"\t" -v fl="$1" 'NR==FNR{len[$1]=$2; next} FNR>1{s=$2-fl; if(s<0)s=0; e=$3+fl; if(e>len[$1])e=len[$1]; print $1"\t"s"\t"e}' "$GENOME.fai" "$HITS" \ + | sort -t "$(printf '\t')" -k1,1 -k2,2n \ + | awk -F"\t" 'BEGIN{OFS="\t"} NR==1{c=$1;s=$2;e=$3;next} $1==c && $2<=e {if($3>e)e=$3; next} {print c,s,e; c=$1;s=$2;e=$3} END{if(NR)print c,s,e}' +} +LOCI="$OUT_DIR/numt_loci_$BUILD.bed" +CORE="$OUT_DIR/numt_core_$BUILD.bed" +mkbed "$FLANK" > "$LOCI" +mkbed 0 > "$CORE" + +NEAR="$OUT_DIR/numt_near_identical_$BUILD.tsv" +head -1 "$HITS" > "$NEAR" +awk -F"\t" -v id="$NEAR_IDENT" -v ln="$NEAR_LEN" 'NR>1 && $11>=id && $9>=ln' "$HITS" | sort -t "$(printf '\t')" -k9,9nr >> "$NEAR" + +n_hits=$(($(wc -l < "$HITS") - 1)) +n_near=$(($(wc -l < "$NEAR") - 1)) +span_loci=$(awk -F"\t" '{s+=$3-$2} END{print s+0}' "$LOCI") +span_core=$(awk -F"\t" '{s+=$3-$2} END{print s+0}' "$CORE") +echo "[summary] $BUILD: $n_hits hits of at least $MIN_LEN bp; $n_near near-identical (identity >= $NEAR_IDENT and >= $NEAR_LEN bp)" +echo "[summary] candidate intervals with $FLANK-bp flanks: $(wc -l < "$LOCI") spanning $span_loci bp; core intervals: $(wc -l < "$CORE") spanning $span_core bp" +echo "[summary] near-identical copies (these are the ones that steal MAPQ from authentic mtDNA reads):" +cat "$NEAR" +echo "[summary] longest 15 hits:" +head -1 "$HITS"; tail -n +2 "$HITS" | sort -t "$(printf '\t')" -k9,9nr | head -15 +echo "[ok] wrote $HITS, $LOCI, $CORE, $NEAR, $MD5" diff --git a/src/bin/mt_allele_counts.sh b/src/bin/mt_allele_counts.sh new file mode 100755 index 0000000000000000000000000000000000000000..16430b10d74094ac2ad9e146017d9a1b4bc8e6bd --- /dev/null +++ b/src/bin/mt_allele_counts.sh @@ -0,0 +1,42 @@ +#!/usr/bin/env bash +#' Per-position, per-strand allele counts on the mitochondrial contig for every realigned BAM. +#' +#' samtools mpileup with base quality >= MIN_BQ, BAQ disabled, unmapped, secondary, QC-fail and +#' duplicate records excluded; the pileup string is parsed into forward (upper case) and reverse +#' (lower case) counts of A, C, G, T, plus deletions (*) and insertion events (+). Writes one +#' .counts.tsv per sample into OUT_COUNTS with columns +#' pos ref depth Af Cf Gf Tf Ar Cr Gr Tr del ins. Runs inside images/debian-nanopore.sif. +#' +#' Environment: BAM_DIRS (colon-separated dirs holding *_realigned.bam) REF_MT MT_CONTIG +#' OUT_COUNTS [MIN_BQ=10] [MIN_MQ=0] +set -euo pipefail +export LC_ALL=C +: "${BAM_DIRS:?}" "${REF_MT:?}" "${MT_CONTIG:?}" "${OUT_COUNTS:?}" +MIN_BQ="${MIN_BQ:-10}"; MIN_MQ="${MIN_MQ:-0}" +mkdir -p "$OUT_COUNTS" +[ -s "$REF_MT.fai" ] || samtools faidx "$REF_MT" +n=0 +for d in $(echo "$BAM_DIRS" | tr ":" " "); do + for bam in "$d"/*_realigned.bam; do + [ -s "$bam" ] || continue + id=$(basename "$bam"); id=${id%%-*} + n=$((n + 1)) + samtools mpileup -f "$REF_MT" -r "$MT_CONTIG" -Q "$MIN_BQ" -q "$MIN_MQ" -B -d 1000000 -a \ + --ff UNMAP,SECONDARY,QCFAIL,DUP "$bam" 2>/dev/null \ + | awk -F"\t" -v OFS="\t" 'BEGIN{print "pos","ref","depth","Af","Cf","Gf","Tf","Ar","Cr","Gr","Tr","del","ins"} + { s=$5; ref=toupper($3); delete c; del=0; ins=0; i=1; L=length(s) + while (i<=L) { ch=substr(s,i,1) + if (ch=="^") { i+=2; continue } + if (ch=="$") { i++; continue } + if (ch=="+" || ch=="-") { i++; num=""; while (substr(s,i,1) ~ /[0-9]/) { num=num substr(s,i,1); i++ } if (ch=="+") ins++; i+=num+0; continue } + if (ch=="*" || ch=="#") { del++; i++; continue } + if (ch==".") c[ref "f"]++; else if (ch==",") c[ref "r"]++ + else if (ch ~ /[ACGT]/) c[ch "f"]++; else if (ch ~ /[acgt]/) c[toupper(ch) "r"]++ + i++ } + dp=c["Af"]+c["Cf"]+c["Gf"]+c["Tf"]+c["Ar"]+c["Cr"]+c["Gr"]+c["Tr"] + print $2, ref, dp, c["Af"]+0, c["Cf"]+0, c["Gf"]+0, c["Tf"]+0, c["Ar"]+0, c["Cr"]+0, c["Gr"]+0, c["Tr"]+0, del, ins }' \ + > "$OUT_COUNTS/$id.counts.tsv" + echo "[$id] $(awk 'NR>1{s+=$3; n++} END{printf "mean depth %.1f over %d positions", s/n, n}' "$OUT_COUNTS/$id.counts.tsv")" + done +done +echo "[ok] $n samples counted into $OUT_COUNTS" diff --git a/src/bin/mt_chimera_check.sh b/src/bin/mt_chimera_check.sh new file mode 100755 index 0000000000000000000000000000000000000000..a25ff94ce81618dbd3341245b5b71fffb95898a1 --- /dev/null +++ b/src/bin/mt_chimera_check.sh @@ -0,0 +1,54 @@ +#!/usr/bin/env bash +#' Characterise the MT-origin reads that the clip filter dropped: what is their unaligned tail? +#' +#' For every rescued run, the reads listed in qc/candidates/.clip_dropped.tsv with origin +#' `mt` are looked up in qc/candidates/.candidates.bam, which still holds their original +#' whole-genome primary records. From those records: the identity of the mitochondrial segment +#' (1 - NM / aligned length), the soft-clipped bases, and the SA tag, whose chromosomes say +#' where the rest of the read aligned in the whole genome (a nuclear chromosome = an +#' inter-molecular chimera; the mt contig = a second mitochondrial fragment; none = a tail +#' that aligns nowhere). The same identity statistic is computed for all MT-origin primaries +#' of the sample as the reference distribution. Runs inside images/debian-nanopore.sif. +#' +#' Environment: RUNS_ROOT (holds _chrM/) RUNS (colon-separated run ids) MT_CONTIG OUT_CHIM +set -euo pipefail +export LC_ALL=C +: "${RUNS_ROOT:?}" "${RUNS:?}" "${MT_CONTIG:?}" "${OUT_CHIM:?}" +mkdir -p "$(dirname "$OUT_CHIM")" +READS="${OUT_CHIM%.tsv}_reads.tsv" +printf 'run\tsample\tqname\tread_len\tmt_aligned\tmt_identity\tclip_bp\tsa_targets\n' > "$READS" +printf 'run\tsample\tdropped_mt_reads\tsa_nuclear\tsa_mt_only\tsa_none\tmedian_identity_dropped\tmedian_identity_all_mt\tmedian_read_len_dropped\tmedian_read_len_all_mt\n' > "$OUT_CHIM" +TMP=$(mktemp -d); trap 'rm -rf "$TMP"' EXIT +for run in $(echo "$RUNS" | tr ":" " "); do + for cd in "$RUNS_ROOT/${run}_chrM"/qc/candidates/*.clip_dropped.tsv; do + id=$(basename "$cd" .clip_dropped.tsv) + cand="$RUNS_ROOT/${run}_chrM/qc/candidates/$id.candidates.bam" + [ -s "$cand" ] || continue + awk -F"\t" '$2=="mt"{print $1}' "$cd" | sort -u > "$TMP/names" + n_drop=$(wc -l < "$TMP/names") + samtools view -F 0x904 "$cand" "$MT_CONTIG" | awk -F"\t" -v OFS="\t" -v run="$run" -v s="$id" -v names="$TMP/names" -v mt="$MT_CONTIG" ' + BEGIN { while ((getline q < names) > 0) drop[q]=1 } + { cig=$6; aln=0; clip=0; nm=0; sa="" + while (match(cig, /^[0-9]+[MIDNSHP=X]/)) { n=substr(cig,1,RLENGTH-1)+0; op=substr(cig,RLENGTH,1); if (op=="M"||op=="I"||op=="D"||op=="="||op=="X") aln+=n; if (op=="S"||op=="H") clip+=n; cig=substr(cig,RLENGTH+1) } + for (i=12;i<=NF;i++) { if ($i ~ /^NM:i:/) nm=substr($i,6)+0; if ($i ~ /^SA:Z:/) sa=substr($i,6) } + ident=(aln>0 ? 1-nm/aln : 0); rl=length($10) + allid[++na]=ident; alllen[na]=rl + if ($1 in drop) { + tg="none"; nuc=0; mto=0 + if (sa != "") { tg=""; n2=split(sa, parts, ";"); for (k=1;k<=n2;k++) { if (parts[k]=="") continue; split(parts[k], f, ","); tg=tg f[1] ":" int(f[2]/100000) ","; if (f[1]==mt) mto=1; else nuc=1 } } + if (nuc) cls="nuclear"; else if (mto) cls="mt_only"; else cls="none" + c[cls]++; dropid[++nd]=ident; droplen[nd]=rl + print run, s, $1, rl, aln, sprintf("%.4f", ident), clip, tg > "'"$TMP"'/reads_part" + } } + function med(a, n, i, j, t) { for (i=2;i<=n;i++) { t=a[i]; j=i-1; while (j>0 && a[j]>t) { a[j+1]=a[j]; j-- } a[j+1]=t } return (n ? a[int((n+1)/2)] : 0) } + END { printf "%s\t%s\t%d\t%d\t%d\t%d\t%.4f\t%.4f\t%d\t%d\n", run, s, nd+0, c["nuclear"]+0, c["mt_only"]+0, c["none"]+0, med(dropid, nd+0), med(allid, na+0), med(droplen, nd+0), med(alllen, na+0) }' >> "$OUT_CHIM" + [ -f "$TMP/reads_part" ] && cat "$TMP/reads_part" >> "$READS" && rm -f "$TMP/reads_part" + echo "[$run $id] dropped MT-origin reads $n_drop" + done +done +echo +echo "[summary]" +cat "$OUT_CHIM" +echo +awk -F"\t" 'NR>1{d+=$3; nu+=$4; mo+=$5; no+=$6} END{printf "all runs: %d dropped MT-origin reads: %d with a supplementary alignment on a nuclear chromosome, %d with supplementary alignments on the mt contig only, %d with no supplementary alignment at all\n", d, nu, mo, no}' "$OUT_CHIM" +echo "[ok] wrote $OUT_CHIM and $READS" diff --git a/src/bin/mt_coverage_profile.R b/src/bin/mt_coverage_profile.R new file mode 100755 index 0000000000000000000000000000000000000000..0eb040d29e87c8bd4484aff9390bbd8d829fcb13 --- /dev/null +++ b/src/bin/mt_coverage_profile.R @@ -0,0 +1,105 @@ +#!/usr/bin/env Rscript +#' Per-window depth profile of a mitochondrial contig from modkit bedMethyl pileups. +#' +#' Reads every `*_modkit_pileup.bed.gz` in --dir through its tabix index, keeps the rows of +#' one modification code on the contig (default 6mA, one row per covered adenine per strand), +#' averages Nvalid_cov per window per sample, and reports windows below --min reads in every +#' sample (a hole shared by all samples is the signature of a NUMT dropout; see the nanopore +#' pipeline notes, #numt-rescue). Exit status 2 when such a hole exists, so a launcher can +#' stop. The sample id is the file name before `_modkit_pileup.bed.gz`, truncated at the first +#' dash, which is the pipeline convention (`-Filtered_primary_mapq_10`), so the +#' samplesheet join works on pipeline output as well as on plain `_modkit_pileup.bed.gz`. +#' Runs in the rocker container with the Bioconductor library used by the study analyses; +#' needs data.table, Rsamtools, GenomicRanges, ggplot2. +#' +#' Usage: +#' Rscript mt_coverage_profile.R --dir --contig chrM --length 16299 \ +#' [--window 250] [--min 3] [--mod a] [--out ] \ +#' [--samplesheet --id-col real_barcode --group-col group [--run-id ]] \ +#' [--features ] [--label "text for the subtitle"] +#' +#' Outputs in --out (default: --dir/mt_coverage_profile): mt_coverage_profile.tsv (one row per +#' window, one column per sample), mt_coverage_gaps.tsv (per-sample windows below --min and +#' the all-sample gap span), mt_coverage_profile.png. + +suppressPackageStartupMessages({ library(data.table); library(Rsamtools); library(GenomicRanges); library(ggplot2) }) + +parse_args <- function(args) { + defaults <- list(dir = NULL, contig = "chrM", length = NULL, window = 250L, min = 3, mod = "a", out = NULL, + samplesheet = NULL, `id-col` = "real_barcode", `group-col` = "group", `run-id` = NULL, features = NULL, label = "") + i <- 1L + while (i <= length(args)) { + key <- sub("^--", "", args[i]) + stopifnot("unknown argument" = key %in% names(defaults), "missing value" = i + 1L <= length(args)) + defaults[[key]] <- args[i + 1L] + i <- i + 2L + } + stopifnot("--dir is required" = !is.null(defaults$dir), "--length (contig length in bp) is required" = !is.null(defaults$length)) + defaults$length <- as.integer(defaults$length); defaults$window <- as.integer(defaults$window); defaults$min <- as.numeric(defaults$min) + if (is.null(defaults$out)) defaults$out <- file.path(defaults$dir, "mt_coverage_profile") + defaults +} + +opt <- parse_args(commandArgs(trailingOnly = TRUE)) +dir.create(opt$out, showWarnings = FALSE, recursive = TRUE) +files <- list.files(opt$dir, pattern = "_modkit_pileup[.]bed[.]gz$", full.names = TRUE) +stopifnot("no *_modkit_pileup.bed.gz in --dir" = length(files) > 0, + "every pileup needs a .tbi index" = all(file.exists(paste0(files, ".tbi")))) +samples <- sub("-.*$", "", sub("_modkit_pileup[.]bed[.]gz$", "", basename(files))) +stopifnot("sample ids are not unique after truncating at the first dash" = !anyDuplicated(samples)) +n_win <- (opt$length - 1L) %/% opt$window + 1L + +read_windows <- function(f, s) { + txt <- unlist(scanTabix(TabixFile(f), param = GRanges(opt$contig, IRanges(1, opt$length))), use.names = FALSE) + if (!length(txt)) return(data.table(sample_id = s, bin = seq_len(n_win) - 1L, depth = 0, rows = 0L)) + d <- fread(text = txt, header = FALSE, sep = "\t", select = c(2, 4, 10), col.names = c("start", "mod", "N")) + d <- d[mod == opt$mod] + d[, bin := start %/% opt$window] + w <- d[, .(depth = mean(N), rows = .N), by = bin] + w <- merge(data.table(bin = seq_len(n_win) - 1L), w, by = "bin", all.x = TRUE) + w[is.na(depth), `:=`(depth = 0, rows = 0L)] + w[, sample_id := s] + w[] +} + +prof <- rbindlist(Map(read_windows, files, samples)) +prof[, start := bin * opt$window + 1L] +wide <- dcast(prof, start ~ sample_id, value.var = "depth") +fwrite(wide, file.path(opt$out, "mt_coverage_profile.tsv"), sep = "\t") + +low <- prof[, .(windows_below_min = sum(depth < opt$min), median_depth = median(depth), min_depth = min(depth)), by = sample_id] +all_low <- prof[, .(all_below = all(depth < opt$min)), by = start][all_below == TRUE, start] +gap <- if (length(all_low)) sprintf("%d-%d", min(all_low), max(all_low) + opt$window - 1L) else "none" +gaps <- rbind(low, data.table(sample_id = "ALL_SAMPLES", windows_below_min = length(all_low), median_depth = NA_real_, min_depth = NA_real_), fill = TRUE) +gaps[, all_sample_gap := gap] +fwrite(gaps, file.path(opt$out, "mt_coverage_gaps.tsv"), sep = "\t") + +groups <- NULL +if (!is.null(opt$samplesheet) && file.exists(opt$samplesheet)) { + ss <- fread(opt$samplesheet) + if (!is.null(opt$`run-id`) && "run_id" %in% names(ss)) ss <- ss[run_id == opt$`run-id`] + if (all(c(opt$`id-col`, opt$`group-col`) %in% names(ss))) groups <- unique(ss[, .(sample_id = get(opt$`id-col`), group = get(opt$`group-col`))]) +} +prof[, group := if (is.null(groups)) "all" else groups$group[match(sample_id, groups$sample_id)]] +prof[is.na(group), group := "unassigned"] + +p <- ggplot(prof, aes(start, pmax(depth, 0.5), colour = group, group = sample_id)) + + geom_step(linewidth = 0.45, alpha = 0.85) + + scale_y_log10(limits = c(0.25, NA), breaks = c(1, 3, 10, 30, 100, 300, 1000)) + + geom_hline(yintercept = opt$min, linetype = "dotted", colour = "grey40") + + labs(title = sprintf("%s coverage per %d-bp window, %d samples (mean Nvalid_cov of %s rows)", opt$contig, opt$window, length(samples), opt$mod), + subtitle = sprintf("%sdotted = %g reads; windows below it in every sample: %s", if (nzchar(opt$label)) paste0(opt$label, "; ") else "", opt$min, gap), + x = sprintf("%s position (bp)", opt$contig), y = "depth (log)", colour = NULL) + + theme_bw(base_size = 11) + theme(legend.position = "top", plot.subtitle = element_text(size = 8.5, colour = "grey30")) +if (length(all_low)) p <- p + annotate("rect", xmin = min(all_low), xmax = max(all_low) + opt$window - 1L, ymin = 0.5, ymax = Inf, alpha = 0.15, fill = "grey30") +if (!is.null(opt$features) && file.exists(opt$features)) { + ft <- fread(opt$features, header = FALSE, select = 1:3, col.names = c("name", "start", "end")) + p <- p + geom_segment(data = ft, aes(x = start, xend = end, y = 0.6, yend = 0.6), inherit.aes = FALSE, colour = "grey20", linewidth = 2) + + geom_text(data = ft, aes(x = (start + end) / 2, y = 0.42, label = name), inherit.aes = FALSE, size = 2.1, angle = 90, hjust = 1) +} +ggsave(file.path(opt$out, "mt_coverage_profile.png"), p, width = 12, height = 5.5, dpi = 180, bg = "white") + +cat(sprintf("[mt_coverage_profile] %d samples, %d windows of %d bp on %s\n", length(samples), n_win, opt$window, opt$contig)) +print(low) +cat(sprintf("[mt_coverage_profile] windows below %g reads in every sample: %d (%s)\n", opt$min, length(all_low), gap)) +if (length(all_low)) { cat("[mt_coverage_profile] shared coverage hole detected: check for a NUMT dropout before any test on this contig\n"); quit(status = 2) } diff --git a/src/bin/mt_deletion_scan.sh b/src/bin/mt_deletion_scan.sh new file mode 100755 index 0000000000000000000000000000000000000000..d27e76ce13bc64feef5cc8af26feb9437d7704fc --- /dev/null +++ b/src/bin/mt_deletion_scan.sh @@ -0,0 +1,85 @@ +#!/usr/bin/env bash +#' First-pass scan for mtDNA deletions in realigned mitochondrial BAMs. +#' +#' A molecule carrying a deletion aligns as two same-strand segments with a reference gap +#' between them, or as one alignment with a long D operation. For every read with at least two +#' alignments on the contig, consecutive segments in read order are paired; a pair whose +#' reference gap is between MIN_DEL and MAX_DEL bp and does not wrap the origin is a junction +#' candidate. D operations of at least MIN_DEL bp inside one alignment are candidates too. +#' Candidates are clustered per sample in BIN-bp windows on both breakpoints; a cluster needs +#' MIN_READS reads to be reported, because a single junction read is indistinguishable from a +#' chimera. Depth at each breakpoint gives a rough deletion fraction. Runs inside +#' images/debian-nanopore.sif. +#' +#' Environment: BAM_DIRS (colon-separated dirs holding *_realigned.bam) MT_CONTIG MT_LEN OUT_DEL +#' [MIN_DEL=50] [MAX_DEL=16000] [BIN=30] [MIN_READS=2] [MAX_QGAP=200] +set -euo pipefail +export LC_ALL=C +: "${BAM_DIRS:?}" "${MT_CONTIG:?}" "${MT_LEN:?}" "${OUT_DEL:?}" +MIN_DEL="${MIN_DEL:-50}"; MAX_DEL="${MAX_DEL:-16000}"; BIN="${BIN:-30}"; MIN_READS="${MIN_READS:-2}"; MAX_QGAP="${MAX_QGAP:-200}" +mkdir -p "$(dirname "$OUT_DEL")" +CAND="${OUT_DEL%.tsv}_candidates.tsv" +printf 'sample\tqname\tkind\tstrand\tbreak_a\tbreak_b\tdel_size\tquery_gap\tmapq_min\n' > "$CAND" +for d in $(echo "$BAM_DIRS" | tr ":" " "); do + for bam in "$d"/*_realigned.bam; do + [ -s "$bam" ] || continue + id=$(basename "$bam"); id=${id%%-*} + samtools view -F 0x104 "$bam" "$MT_CONTIG" | awk -F"\t" -v OFS="\t" -v s="$id" -v L="$MT_LEN" -v mind="$MIN_DEL" -v maxd="$MAX_DEL" ' + function parse(cig, ln, op, lead, tail, q, r, seen) { + lead=0; tail=0; q=0; r=0; seen=0 + while (match(cig, /^[0-9]+[MIDNSHP=X]/)) { + ln=substr(cig,1,RLENGTH-1)+0; op=substr(cig,RLENGTH,1); cig=substr(cig,RLENGTH+1) + if (op=="S"||op=="H") { if (!seen) lead+=ln; else tail+=ln } + else { seen=1 + if (op=="M"||op=="="||op=="X") { q+=ln; r+=ln } + else if (op=="I") q+=ln + else if (op=="D"||op=="N") { if (ln>=mind) print s, qn, "cigar_D", strand, pos+r-1, pos+r+ln, ln, 0, mq; r+=ln } } } + LEAD=lead; TAIL=tail; QLEN=q; RLEN=r } + { qn=$1; flag=$2; pos=$4; mq=$5; strand=(int(flag/16)%2 ? "-" : "+") + parse($6) + readlen=LEAD+QLEN+TAIL + if (strand=="+") { qs=LEAD; qe=LEAD+QLEN } else { qs=TAIL; qe=TAIL+QLEN } + k=qn; cnt[k]++; i=cnt[k] + S[k,i]=strand; QS[k,i]=qs; QE[k,i]=qe; RB[k,i]=pos; REN[k,i]=pos+RLEN-1; MQ[k,i]=mq } + END { + for (k in cnt) { if (cnt[k]<2) continue + m=cnt[k]; for (i=1;i<=m;i++) idx[i]=i + for (i=2;i<=m;i++) { t=idx[i]; j=i-1; while (j>0 && QS[k,idx[j]]>QS[k,t]) { idx[j+1]=idx[j]; j-- } idx[j+1]=t } + for (i=1;iL-300 && bb<300) continue + if (del>=mind && del<=maxd) { mqm=(MQ[k,a]> "$CAND" + echo "[$id] scanned" + done +done +printf 'sample\tbreak_a_bin\tbreak_b_bin\tdel_size_median\tn_reads\tn_split\tn_cigar\tquery_gap_median\tdepth_a\tdepth_b\tdeletion_fraction\n' > "$OUT_DEL" +awk -F"\t" -v OFS="\t" -v bin="$BIN" -v minr="$MIN_READS" -v maxq="$MAX_QGAP" 'NR>1 && ($8<=maxq) { + ka=int($5/bin)*bin; kb=int($6/bin)*bin; key=$1"\t"ka"\t"kb + c[key]++; if ($3=="split") sp[key]++; else cg[key]++ + d[key, c[key]]=$7; g[key, c[key]]=$8 } + function med(key, m, i, j, t) { split("", arr); for (i=1;i<=m;i++) arr[i]=d[key,i]; for (i=2;i<=m;i++) { t=arr[i]; j=i-1; while (j>0 && arr[j]>t) { arr[j+1]=arr[j]; j-- } arr[j+1]=t } return arr[int((m+1)/2)] } + function medg(key, m, i, j, t) { split("", arr2); for (i=1;i<=m;i++) arr2[i]=g[key,i]; for (i=2;i<=m;i++) { t=arr2[i]; j=i-1; while (j>0 && arr2[j]>t) { arr2[j+1]=arr2[j]; j-- } arr2[j+1]=t } return arr2[int((m+1)/2)] } + END { for (key in c) if (c[key]>=minr) print key, med(key, c[key]), c[key], sp[key]+0, cg[key]+0, medg(key, c[key]), "NA", "NA", "NA" }' "$CAND" | sort -t "$(printf '\t')" -k1,1 -k5,5nr > "$OUT_DEL.tmp" +while IFS=$'\t' read -r s ka kb dmed nreads nsplit ncig gmed da db frac; do + bam=$(ls $(echo "$BAM_DIRS" | tr ":" " ") 2>/dev/null | grep -E "^${s}-.*_realigned.bam$" | head -n 1) + bamp=""; for d in $(echo "$BAM_DIRS" | tr ":" " "); do [ -s "$d/$bam" ] && bamp="$d/$bam"; done + if [ -n "$bamp" ]; then + da=$(samtools depth -a -r "$MT_CONTIG:$((ka+1))-$((ka+1))" "$bamp" | awk '{print $3+0}'); db=$(samtools depth -a -r "$MT_CONTIG:$((kb+1))-$((kb+1))" "$bamp" | awk '{print $3+0}') + frac=$(awk -v n="$nreads" -v a="${da:-0}" -v b="${db:-0}" 'BEGIN{m=(a+b)/2; printf "%.3f", (m>0 ? n/(m+n) : 0)}') + fi + printf '%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n' "$s" "$ka" "$kb" "$dmed" "$nreads" "$nsplit" "$ncig" "$gmed" "${da:-NA}" "${db:-NA}" "$frac" >> "$OUT_DEL" +done < "$OUT_DEL.tmp" +rm -f "$OUT_DEL.tmp" +echo +echo "[candidates] junction reads by sample (all, before clustering):" +awk -F"\t" 'NR>1{c[$1]++; if($3=="split") s[$1]++} END{for(k in c) printf "%s: %d candidate reads (%d split, %d cigar)\n", k, c[k], s[k]+0, c[k]-s[k]}' "$CAND" | sort +echo +echo "[clusters] junctions with at least $MIN_READS reads (breakpoints in $BIN-bp bins, query gap <= $MAX_QGAP):" +cat "$OUT_DEL" +echo +echo "[recurrent] breakpoint pairs seen in more than one sample:" +awk -F"\t" 'NR>1{k=$2"-"$3; n[k]++; r[k]+=$5; sz[k]=$4} END{for(k in n) if(n[k]>1) printf "%s (del %s bp): %d samples, %d reads\n", k, sz[k], n[k], r[k]}' "$OUT_DEL" | sort -t" " -k4,4nr +echo "[ok] wrote $OUT_DEL and $CAND" diff --git a/src/bin/mt_depth_compare.sh b/src/bin/mt_depth_compare.sh new file mode 100755 index 0000000000000000000000000000000000000000..6d0c279fd95ba9e74511f417fa33ddbd519fb825 --- /dev/null +++ b/src/bin/mt_depth_compare.sh @@ -0,0 +1,48 @@ +#!/usr/bin/env bash +#' Depth of the mitochondrial contig before and after the rescue, per window and per site, +#' computed from BAMs so that it can run before any pileup exists. +#' +#' Before is what step 2 keeps: records with MAPQ >= MAPQ on the mitochondrial contig of the +#' unfiltered whole-genome BAM (secondary records excluded, supplementary kept, as in +#' filter_bam). After is the realigned mitochondrial BAM written by rescue_mt.sh. +#' Writes /window_depth_before_after.tsv and /site_depth_before_after.tsv. +#' Runs inside images/debian-nanopore.sif. +#' +#' Environment: SRC_DIR (-Unfiltered.bam) AFTER_BAM_DIR (-Filtered__realigned.bam) +#' MT_CONTIG OUT_DEPTH [MAPQ=10] [WINDOW=250] [SITES colon-separated 1-based positions] +set -euo pipefail +export LC_ALL=C +: "${SRC_DIR:?}" "${AFTER_BAM_DIR:?}" "${MT_CONTIG:?}" "${OUT_DEPTH:?}" +MAPQ="${MAPQ:-10}" +WINDOW="${WINDOW:-250}" +SITES="${SITES:-}" +mkdir -p "$OUT_DEPTH" +WIN="$OUT_DEPTH/window_depth_before_after.tsv" +SITE="$OUT_DEPTH/site_depth_before_after.tsv" +printf 'sample\twindow_start\tbefore_depth\tafter_depth\n' > "$WIN" +printf 'sample\tposition\tbefore_depth\tafter_depth\n' > "$SITE" +windows() { + samtools depth -a -r "$MT_CONTIG" -Q "$1" "$2" \ + | awk -v w="$WINDOW" 'BEGIN{OFS="\t"}{b=int(($2-1)/w); s[b]+=$3; n[b]++} END{for(b in s) print b*w+1, s[b]/n[b]}' | sort -k1,1n +} +site() { samtools depth -a -r "$MT_CONTIG:$3-$3" -Q "$1" "$2" | awk '{d=$3+0} END{print d+0}'; } +n=0 +for ubam in "$SRC_DIR"/*-Unfiltered.bam; do + id=$(basename "$ubam" -Unfiltered.bam) + abam="$AFTER_BAM_DIR/$id-Filtered_${MT_CONTIG}_realigned.bam" + if [ ! -s "$abam" ]; then echo "[warn] no realigned BAM for $id"; continue; fi + n=$((n + 1)) + windows "$MAPQ" "$ubam" > "$OUT_DEPTH/$id.before.tmp" + windows 0 "$abam" > "$OUT_DEPTH/$id.after.tmp" + awk -v s="$id" 'BEGIN{OFS="\t"} NR==FNR{b[$1]=$2; next} {print s,$1,b[$1]+0,$2}' "$OUT_DEPTH/$id.before.tmp" "$OUT_DEPTH/$id.after.tmp" >> "$WIN" + for p in $(echo "$SITES" | tr ":" " "); do + printf '%s\t%s\t%s\t%s\n' "$id" "$p" "$(site "$MAPQ" "$ubam" "$p")" "$(site 0 "$abam" "$p")" >> "$SITE" + done + rm -f "$OUT_DEPTH/$id.before.tmp" "$OUT_DEPTH/$id.after.tmp" +done +echo "[info] $n samples compared (before = $MT_CONTIG records at MAPQ >= $MAPQ in the unfiltered BAM; after = realigned BAM)" +echo "[sites] depth summed over samples, before -> after:" +awk -F"\t" 'NR>1{b[$2]+=$3; a[$2]+=$4} END{for(p in b) print p, b[p], "->", a[p]}' "$SITE" | sort -k1,1n +echo "[windows] mean depth over samples per window, before -> after:" +awk -F"\t" 'NR>1{b[$2]+=$3; a[$2]+=$4; n[$2]++} END{for(w in b) printf "%d\t%.1f\t%.1f\t%.2f\n", w, b[w]/n[w], a[w]/n[w], (b[w]>0 ? a[w]/b[w] : 0)}' "$WIN" | sort -k1,1n +echo "[ok] wrote $WIN and $SITE" diff --git a/src/bin/mt_pileup_check.sh b/src/bin/mt_pileup_check.sh new file mode 100755 index 0000000000000000000000000000000000000000..2e627afd171a3c6e9ac6032bab25710cd8a28f89 --- /dev/null +++ b/src/bin/mt_pileup_check.sh @@ -0,0 +1,55 @@ +#!/usr/bin/env bash +#' Checklist items on rescued pileups: row counts per modification code and strand, and the +#' modkit pass thresholds of the rescued run against the original run of the same samples. +#' +#' For every *_modkit_pileup.bed.gz in AFTER_DIR, counts rows on the mitochondrial contig by +#' modification code and strand (one row per covered base per strand, so the count of code +#' `a` rows approaches the number of adenines on both strands when coverage is complete), and +#' reads pass_threshold lines from the after summary and, when present, from the matching +#' before summary in BEFORE_DIR (sample matched by the token before the first dash). +#' Writes . Runs inside images/debian-nanopore.sif. +#' +#' Also, per sample, the mean per-site read accounting of the `a` rows (valid, fail, diff, +#' nocall, delete) before and after, which separates a threshold shift from lost reads, and, +#' when AFTER_BAM_DIR and SRC_DIR are given, the share of primary reads carrying an MM tag in +#' the realigned BAM and on the contig of the unfiltered BAM. +#' +#' Environment: AFTER_DIR MT_CONTIG OUT_CHECK [BEFORE_DIR] [REF_MT: counts A/C/G/T on the contig] +#' [AFTER_BAM_DIR (-Filtered__realigned.bam)] [SRC_DIR (-Unfiltered.bam)] +set -euo pipefail +export LC_ALL=C +: "${AFTER_DIR:?}" "${MT_CONTIG:?}" "${OUT_CHECK:?}" +BEFORE_DIR="${BEFORE_DIR:-}" +REF_MT="${REF_MT:-}" +AFTER_BAM_DIR="${AFTER_BAM_DIR:-}" +SRC_DIR="${SRC_DIR:-}" +mkdir -p "$(dirname "$OUT_CHECK")" +if [ -n "$REF_MT" ] && [ -s "$REF_MT" ]; then + echo "[reference] base counts on $MT_CONTIG (both strands: adenines = A + T, cytosines = C + G):" + grep -v "^>" "$REF_MT" | tr -d "\n\r" | tr "a-z" "A-Z" | fold -w1 | sort | uniq -c | tr "\n" " "; echo +fi +thr() { grep -E "pass_threshold" "$1" 2>/dev/null | tr -s " " | awk '{printf "%s=%s;", $2, $3}' || true; } +find_one() { local c; for c in "$@"; do if [ -f "$c" ]; then echo "$c"; return 0; fi; done; return 0; } +acct() { tabix "$1" "$MT_CONTIG" | awk -F"\t" '$4=="a"{n++; v+=$10; f+=$16; d+=$17; nc+=$18; del+=$15} END{if(n) printf "%.1f\t%.1f\t%.1f\t%.1f\t%.1f", v/n, f/n, d/n, nc/n, del/n; else printf "NA\tNA\tNA\tNA\tNA"}'; } +tagshare() { samtools view -F 0x904 "$@" | awk '{m=0; for(i=12;i<=NF;i++) if($i ~ /^MM:Z:/){m=1; break}; t++; if(m) k++} END{printf "%d/%d", k+0, t+0}'; } +printf 'sample\trows_total\ta_plus\ta_minus\tm_plus\tm_minus\th_plus\th_minus\tafter_thresholds\tbefore_thresholds\tafter_a_valid\tafter_a_fail\tafter_a_diff\tafter_a_nocall\tafter_a_del\tbefore_a_valid\tbefore_a_fail\tbefore_a_diff\tbefore_a_nocall\tbefore_a_del\tmm_tag_after\tmm_tag_before\n' > "$OUT_CHECK" +for f in "$AFTER_DIR"/*_modkit_pileup.bed.gz; do + id=$(basename "$f" _modkit_pileup.bed.gz); tok=${id%%-*} + counts=$(tabix "$f" "$MT_CONTIG" | awk -F"\t" '{n++; c[$4" "$6]++} END{printf "%d\t%d\t%d\t%d\t%d\t%d\t%d", n, c["a +"], c["a -"], c["m +"], c["m -"], c["h +"], c["h -"]}') + after=$(thr "$AFTER_DIR/${id}_modkit_summary.txt") + aacc=$(acct "$f") + before=""; bacc=$(printf 'NA\tNA\tNA\tNA\tNA') + if [ -n "$BEFORE_DIR" ]; then + bf=$(find_one "$BEFORE_DIR/$tok"_modkit_summary.txt "$BEFORE_DIR/$tok"-*_modkit_summary.txt) + if [ -n "$bf" ]; then before=$(thr "$bf"); fi + bp=$(find_one "$BEFORE_DIR/$tok"_modkit_pileup.bed.gz "$BEFORE_DIR/$tok"-*_modkit_pileup.bed.gz) + if [ -n "$bp" ]; then bacc=$(acct "$bp"); fi + fi + ta="NA"; tb="NA" + if [ -n "$AFTER_BAM_DIR" ] && [ -s "$AFTER_BAM_DIR/$tok-Filtered_${MT_CONTIG}_realigned.bam" ]; then ta=$(tagshare "$AFTER_BAM_DIR/$tok-Filtered_${MT_CONTIG}_realigned.bam"); fi + if [ -n "$SRC_DIR" ] && [ -s "$SRC_DIR/$tok-Unfiltered.bam" ]; then tb=$(tagshare "$SRC_DIR/$tok-Unfiltered.bam" "$MT_CONTIG"); fi + printf '%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n' "$tok" "$counts" "${after:-NA}" "${before:-NA}" "$aacc" "$bacc" "$ta" "$tb" >> "$OUT_CHECK" +done +echo "[pileup rows and thresholds]" +cat "$OUT_CHECK" +echo "[ok] wrote $OUT_CHECK" diff --git a/src/bin/mt_read_accounting.sh b/src/bin/mt_read_accounting.sh new file mode 100755 index 0000000000000000000000000000000000000000..401016221f88ffa7a257441192513a30aeba5d99 --- /dev/null +++ b/src/bin/mt_read_accounting.sh @@ -0,0 +1,64 @@ +#!/usr/bin/env bash +#' Read accounting for the mitochondrial rescue: what the unfiltered BAMs contain, which +#' flowcells they pool, and how the realigned reads distribute by aligned fraction. +#' +#' Per unfiltered BAM: read groups (one per flowcell or basecalling run), mapped records in +#' total and on the mitochondrial contig (idxstats), primary reads (flagstat). Per realigned +#' BAM: for every read, origin (read group set by rescue_mt.sh), read length, aligned fraction +#' summed over primary and supplementary records, and soft-clipped bases; tabulated as +#' origin x aligned-fraction bin and, for reads below MIN_ALN_FRAC, origin x clip-size bin. +#' Runs inside images/debian-nanopore.sif. +#' +#' Environment: SRC_DIR (-Unfiltered.bam) AFTER_BAM_DIR (-Filtered__realigned.bam) +#' MT_CONTIG OUT_ACCT [MIN_ALN_FRAC=0.8] [THREADS=2] +set -euo pipefail +export LC_ALL=C +: "${SRC_DIR:?}" "${AFTER_BAM_DIR:?}" "${MT_CONTIG:?}" "${OUT_ACCT:?}" +MIN_ALN_FRAC="${MIN_ALN_FRAC:-0.8}" +T="${THREADS:-2}" +mkdir -p "$OUT_ACCT" +ACC="$OUT_ACCT/read_accounting.tsv" +READS="$OUT_ACCT/realigned_reads.tsv" +printf 'sample\tread_groups\tflowcells\tprimary_reads\tmapped_records\tmt_records\tmt_pct\n' > "$ACC" +printf 'sample\torigin\tread_len\taligned_frac\tclip_bp\tn_records\n' > "$READS" +for ubam in "$SRC_DIR"/*-Unfiltered.bam; do + id=$(basename "$ubam" -Unfiltered.bam) + abam="$AFTER_BAM_DIR/$id-Filtered_${MT_CONTIG}_realigned.bam" + rg=$(samtools view -H "$ubam" | grep -c "^@RG" || true) + fc=$(samtools view -H "$ubam" | awk -F"\t" '$1=="@RG"{for(i=2;i<=NF;i++) if($i ~ /^(PU|ID):/) {print $i; break}}' | sort -u | tr "\n" "," | sed "s/,\$//") + pri=$(samtools flagstat -@ "$T" "$ubam" | awk '/ primary$/{print $1}') + read -r mapped mt <<< "$(samtools idxstats "$ubam" | awk -v c="$MT_CONTIG" '{m+=$3; if($1==c) x=$3} END{print m, x+0}')" + pct=$(awk -v a="$mt" -v b="$mapped" 'BEGIN{printf "%.4f", (b>0 ? 100*a/b : 0)}') + printf '%s\t%s\t%s\t%s\t%s\t%s\t%s\n' "$id" "$rg" "$fc" "$pri" "$mapped" "$mt" "$pct" >> "$ACC" + echo "[$id] read groups $rg ($fc); primary reads $pri; mapped records $mapped; on $MT_CONTIG $mt ($pct percent)" + if [ -s "$abam" ]; then + samtools view -F 0x104 "$abam" | awk -F"\t" -v s="$id" 'BEGIN{OFS="\t"} + { rg="NA"; for (i=12;i<=NF;i++) if ($i ~ /^RG:Z:/) { rg=substr($i,6); sub(/^.*\./,"",rg); break } + cig=$6; aln=0; clip=0 + while (match(cig, /^[0-9]+[MIDNSHP=X]/)) { n=substr(cig,1,RLENGTH-1)+0; op=substr(cig,RLENGTH,1); if (op=="M"||op=="I"||op=="="||op=="X") aln+=n; else if (op=="S"||op=="H") clip+=n; cig=substr(cig,RLENGTH+1) } + A[$1]+=aln; N[$1]++ + if (int($2/2048)%2==0) { L[$1]=aln+clip; R[$1]=rg } } + END { for (q in L) if (L[q]>0) print s, R[q], L[q], sprintf("%.3f", A[q]/L[q]), L[q]-A[q], N[q] }' >> "$READS" + fi +done +echo +echo "[accounting]" +cat "$ACC" +echo +echo "[aligned fraction by origin, all samples pooled] bins: <0.5 0.5-0.6 0.6-0.7 0.7-0.8 0.8-0.9 0.9-0.95 >=0.95" +awk -F"\t" 'NR>1{f=$4+0; b=(f<0.5)?1:(f<0.6)?2:(f<0.7)?3:(f<0.8)?4:(f<0.9)?5:(f<0.95)?6:7; c[$2,b]++; t[$2]++} + END{for(o in t){printf "%-5s n=%-6d", o, t[o]; for(b=1;b<=7;b++) printf " %6d", c[o,b]+0; printf " (fractions:"; for(b=1;b<=7;b++) printf " %.3f", (c[o,b]+0)/t[o]; print ")"}}' "$READS" +echo +echo "[reads below $MIN_ALN_FRAC: unaligned bases per read by origin] bins: <100 100-300 300-1000 >=1000 bp" +awk -F"\t" -v m="$MIN_ALN_FRAC" 'NR>1 && $4=5000 bp" +awk -F"\t" -v m="$MIN_ALN_FRAC" 'NR>1 && $4=5000 bp; median length" +awk -F"\t" 'NR>1{l=$3+0; b=(l<500)?1:(l<1000)?2:(l<2000)?3:(l<5000)?4:5; c[$2,b]++; t[$2]++; len[$2, t[$2]]=l} + END{for(o in t){printf "%-5s n=%-6d", o, t[o]; for(b=1;b<=5;b++) printf " %6d", c[o,b]+0; n=t[o]; for(i=1;i<=n;i++) a[i]=len[o,i]; asort_ok=0; print ""}}' "$READS" +awk -F"\t" 'NR>1{print $2"\t"$3}' "$READS" | sort -t "$(printf '\t')" -k1,1 -k2,2n | awk -F"\t" '{v[$1,++n[$1]]=$2} END{for(o in n){m=int((n[o]+1)/2); printf "%s median read length %d bp (n=%d)\n", o, v[o,m], n[o]}}' +echo "[ok] wrote $ACC and $READS" diff --git a/src/bin/mt_read_bases.sh b/src/bin/mt_read_bases.sh new file mode 100755 index 0000000000000000000000000000000000000000..63352c49b57f0feffab459c0fe0a66608e1f67c4 --- /dev/null +++ b/src/bin/mt_read_bases.sh @@ -0,0 +1,43 @@ +#!/usr/bin/env bash +#' Per-read bases at called heteroplasmic sites, for phasing. +#' +#' For every .sites.bed in SITES_DIR (written by mt_spectrum_phasing.R, stage call), runs +#' samtools mpileup restricted to those positions with --output-QNAME and writes one row per +#' read and site: pos qname base strand. Bases are strand-resolved to the reference strand +#' (upper case); deletions and reference skips are recorded as "-". Runs inside +#' images/debian-nanopore.sif. +#' +#' Environment: BAM_DIRS (colon-separated dirs holding *_realigned.bam) REF_MT MT_CONTIG +#' SITES_DIR OUT_BASES [MIN_BQ=10] +set -euo pipefail +export LC_ALL=C +: "${BAM_DIRS:?}" "${REF_MT:?}" "${MT_CONTIG:?}" "${SITES_DIR:?}" "${OUT_BASES:?}" +MIN_BQ="${MIN_BQ:-10}" +mkdir -p "$OUT_BASES" +n=0 +for bed in "$SITES_DIR"/*.sites.bed; do + [ -s "$bed" ] || continue + id=$(basename "$bed" .sites.bed) + bam="" + for d in $(echo "$BAM_DIRS" | tr ":" " "); do for f in "$d/$id"-*_realigned.bam; do if [ -s "$f" ]; then bam="$f"; fi; done; done + [ -n "$bam" ] || { echo "[warn] no BAM for $id"; continue; } + n=$((n + 1)) + samtools mpileup -f "$REF_MT" -r "$MT_CONTIG" -l "$bed" -Q "$MIN_BQ" -q 0 -B -d 1000000 --output-QNAME \ + --ff UNMAP,SECONDARY,QCFAIL,DUP "$bam" 2>/dev/null \ + | awk -F"\t" -v OFS="\t" 'BEGIN{print "pos","qname","base","strand"} + { s=$5; ref=toupper($3); nq=split($7, q, ","); k=0; i=1; L=length(s) + while (i<=L) { ch=substr(s,i,1) + if (ch=="^") { i+=2; continue } + if (ch=="$") { i++; continue } + if (ch=="+" || ch=="-") { i++; num=""; while (substr(s,i,1) ~ /[0-9]/) { num=num substr(s,i,1); i++ } i+=num+0; continue } + k++ + if (ch=="*" || ch=="#") { b="-"; st=(ch=="*"?"+":"-") } + else if (ch==".") { b=ref; st="+" } else if (ch==",") { b=ref; st="-" } + else if (ch ~ /[ACGTN]/) { b=ch; st="+" } else if (ch ~ /[acgtn]/) { b=toupper(ch); st="-" } + else if (ch==">" || ch=="<") { b="-"; st="." } else { b="?"; st="." } + if (k<=nq) print $2, q[k], b, st + i++ } }' \ + > "$OUT_BASES/$id.readbases.tsv" + echo "[$id] $(($(wc -l < "$OUT_BASES/$id.readbases.tsv") - 1)) read-site observations at $(wc -l < "$bed") sites" +done +echo "[ok] $n samples with per-read bases in $OUT_BASES" diff --git a/src/bin/mt_site_compare.sh b/src/bin/mt_site_compare.sh new file mode 100755 index 0000000000000000000000000000000000000000..860bebd4b281a84a11ca3a5a296dc7efced2fc69 --- /dev/null +++ b/src/bin/mt_site_compare.sh @@ -0,0 +1,47 @@ +#!/usr/bin/env bash +#' Compare modkit pileups before and after the mitochondrial rescue at named sites and regions. +#' +#' For every *_modkit_pileup.bed.gz in AFTER_DIR the matching pileup in BEFORE_DIR is found by +#' the sample token (the part of the file name before the first dash, as in the pipeline). +#' Writes : one row per sample, site and strand with valid coverage and percent modified +#' before and after; and _regions.tsv: one row per sample and region with the +#' number of covered sites, mean valid coverage and sites at 3 reads or more, before and after. +#' Runs inside images/debian-nanopore.sif (tabix). +#' +#' Environment: BEFORE_DIR AFTER_DIR MT_CONTIG OUT_SITES [SITES colon-separated 1-based positions] +#' [REGIONS colon-separated name=start-end, 1-based inclusive] [MOD=a] +set -euo pipefail +export LC_ALL=C +: "${BEFORE_DIR:?}" "${AFTER_DIR:?}" "${MT_CONTIG:?}" "${OUT_SITES:?}" +SITES="${SITES:-}" +REGIONS="${REGIONS:-}" +MOD="${MOD:-a}" +mkdir -p "$(dirname "$OUT_SITES")" +OUT_SITES_REG="${OUT_SITES%.tsv}_regions.tsv" +find_before() { local c; for c in "$BEFORE_DIR/$1"_modkit_pileup.bed.gz "$BEFORE_DIR/$1"-*_modkit_pileup.bed.gz; do if [ -f "$c" ]; then echo "$c"; return 0; fi; done; return 0; } +site_row() { tabix "$1" "$MT_CONTIG:$2-$2" | awk -F"\t" -v m="$MOD" -v s="$3" '$4==m && $6==s {print $10"\t"$11; f=1} END{if(!f) print "0\tNA"}'; } +region_row() { tabix "$1" "$MT_CONTIG:$2-$3" | awk -F"\t" -v m="$MOD" '$4==m {n++; c+=$10; if($10>=3) k++} END{printf "%d\t%.1f\t%d", n+0, (n ? c/n : 0), k+0}'; } +printf 'sample\tposition\tstrand\tbefore_cov\tbefore_pct\tafter_cov\tafter_pct\n' > "$OUT_SITES" +printf 'sample\tregion\tspan\tbefore_sites\tbefore_mean_cov\tbefore_sites_ge3\tafter_sites\tafter_mean_cov\tafter_sites_ge3\n' > "$OUT_SITES_REG" +n=0 +for f in "$AFTER_DIR"/*_modkit_pileup.bed.gz; do + id=$(basename "$f" _modkit_pileup.bed.gz); id=${id%%-*} + b=$(find_before "$id") + if [ -z "$b" ]; then echo "[warn] no before pileup for $id in $BEFORE_DIR"; continue; fi + n=$((n + 1)) + for s in $(echo "$SITES" | tr ":" " "); do + for strand in + -; do + printf '%s\t%s\t%s\t%s\t%s\n' "$id" "$s" "$strand" "$(site_row "$b" "$s" "$strand")" "$(site_row "$f" "$s" "$strand")" >> "$OUT_SITES" + done + done + for r in $(echo "$REGIONS" | tr ":" " "); do + name=${r%%=*}; span=${r#*=}; st=${span%-*}; en=${span#*-} + printf '%s\t%s\t%s\t%s\t%s\n' "$id" "$name" "$span" "$(region_row "$b" "$st" "$en")" "$(region_row "$f" "$st" "$en")" >> "$OUT_SITES_REG" + done +done +echo "[info] $n samples compared" +echo "[sites] pooled valid coverage per site and strand, before -> after (sum over samples):" +awk -F"\t" 'NR>1{b[$2" "$3]+=$4; a[$2" "$3]+=$6} END{for(k in b) print k, b[k], "->", a[k]}' "$OUT_SITES" | sort -k1,1n -k2,2 +echo "[regions] mean valid coverage per region, before -> after (mean over samples):" +awk -F"\t" 'NR>1{b[$2]+=$5; a[$2]+=$8; n[$2]++} END{for(k in b) printf "%s %.1f -> %.1f\n", k, b[k]/n[k], a[k]/n[k]}' "$OUT_SITES_REG" | sort +echo "[ok] wrote $OUT_SITES and $OUT_SITES_REG" diff --git a/src/bin/mt_spectrum_phasing.R b/src/bin/mt_spectrum_phasing.R new file mode 100755 index 0000000000000000000000000000000000000000..4bb21a96507e9d3617b5eeeff81446b2a5e51dff --- /dev/null +++ b/src/bin/mt_spectrum_phasing.R @@ -0,0 +1,248 @@ +#!/usr/bin/env Rscript +#' Heteroplasmic sites, mutational spectrum and per-molecule phasing on mitochondrial BAMs. +#' +#' Stage `call` reads the per-position allele counts written by mt_allele_counts.sh, calls +#' heteroplasmic and homoplasmic sites per sample, flags systematic-error hotspots (a minor +#' allele recurring at low frequency across most samples), homopolymer context and barcode +#' cross-talk candidates (a minor allele that is homoplasmic in another sample of the same run), +#' tabulates the substitution spectrum per sample and per group, tests groups with depth as a +#' covariate, writes one BED of called sites per sample for mt_read_bases.sh, and draws the +#' figures. Stage `phase` reads the per-read bases at those sites and tests, for every pair of +#' heteroplasmic sites in a sample, whether the minor alleles co-occur on the same molecules +#' (Fisher's exact test on the 2 by 2 table of reads covering both sites). +#' Runs in the rocker container with the study library (data.table, ggplot2). +#' +#' Usage: +#' Rscript mt_spectrum_phasing.R --stage call --counts --ref --contig MT \ +#' --samplesheet --id-col real_barcode --group-col group --run-col run_id --out \ +#' [--min-depth 20] [--min-vaf 0.02] [--max-vaf 0.9] [--min-alt-strand 3] [--hotspot-frac 0.5] +#' [--crosstalk-vaf 0.1] [--mixture-sites 10] [--cohort label] +#' Rscript mt_spectrum_phasing.R --stage phase --out --bases [--min-cover 20] + +suppressPackageStartupMessages({ library(data.table); library(ggplot2) }) + +parse_args <- function(args) { + d <- list(stage = "call", counts = NULL, ref = NULL, contig = "MT", samplesheet = NULL, `id-col` = "real_barcode", + `group-col` = "group", `run-col` = "run_id", out = NULL, bases = NULL, `min-depth` = 20, `min-vaf` = 0.02, + `min-alt-strand` = 3, `hotspot-frac` = 0.5, `min-cover` = 20, `mixture-sites` = 10, `max-vaf` = 0.9, `crosstalk-vaf` = 0.1, cohort = "cohort") + i <- 1L + while (i <= length(args)) { + key <- sub("^--", "", args[i]); stopifnot("unknown argument" = key %in% names(d), "missing value" = i + 1L <= length(args)) + d[[key]] <- args[i + 1L]; i <- i + 2L + } + for (k in c("min-depth", "min-vaf", "min-alt-strand", "hotspot-frac", "min-cover", "mixture-sites", "max-vaf", "crosstalk-vaf")) d[[k]] <- as.numeric(d[[k]]) + stopifnot("--out is required" = !is.null(d$out)) + d +} +opt <- parse_args(commandArgs(trailingOnly = TRUE)) +dir.create(opt$out, showWarnings = FALSE, recursive = TRUE) +classes6 <- c("C>T/G>A", "T>C/A>G", "C>A/G>T", "C>G/G>C", "T>A/A>T", "T>G/A>C") +collapse6 <- function(ref, alt) { + key <- paste0(ref, ">", alt) + map <- c("C>T" = 1, "G>A" = 1, "T>C" = 2, "A>G" = 2, "C>A" = 3, "G>T" = 3, "C>G" = 4, "G>C" = 4, "T>A" = 5, "A>T" = 5, "T>G" = 6, "A>C" = 6) + factor(classes6[map[key]], levels = classes6) +} + +if (opt$stage == "call") { + stopifnot(!is.null(opt$counts), !is.null(opt$ref)) + files <- list.files(opt$counts, pattern = "[.]counts[.]tsv$", full.names = TRUE) + stopifnot("no counts files" = length(files) > 0) + cnt <- rbindlist(lapply(files, function(f) { d <- fread(f); d[, sample_id := sub("[.]counts[.]tsv$", "", basename(f))]; d })) + ref_seq <- toupper(paste(sub("\r$", "", grep("^>", readLines(opt$ref), value = TRUE, invert = TRUE)), collapse = "")) + L <- nchar(ref_seq) + bases <- strsplit(ref_seq, "")[[1]] + run_id <- rle(bases); ends <- cumsum(run_id$lengths); starts <- ends - run_id$lengths + 1L + hp <- rep(FALSE, L) + for (k in which(run_id$lengths >= 4)) hp[max(1L, starts[k] - 1L):min(L, ends[k] + 1L)] <- TRUE + hp_dt <- data.table(pos = seq_len(L), homopolymer = hp) + + groups <- NULL + if (!is.null(opt$samplesheet) && file.exists(opt$samplesheet)) { + ss <- fread(opt$samplesheet) + cols <- c(opt$`id-col`, opt$`group-col`) + if (all(cols %in% names(ss))) { + groups <- unique(ss[, .(sample_id = as.character(get(opt$`id-col`)), group = as.character(get(opt$`group-col`)), + run = if (opt$`run-col` %in% names(ss)) as.character(get(opt$`run-col`)) else "run")]) + groups <- groups[!duplicated(sample_id)] + } + } + if (is.null(groups)) groups <- data.table(sample_id = unique(cnt$sample_id), group = "all", run = "run") + cnt <- merge(cnt, groups, by = "sample_id", all.x = TRUE) + cnt[is.na(group), group := "unassigned"]; cnt[is.na(run), run := "run"] + cnt <- merge(cnt, hp_dt, by = "pos", all.x = TRUE) + + long <- melt(cnt, id.vars = c("sample_id", "group", "run", "pos", "ref", "depth", "homopolymer"), + measure.vars = list(fwd = c("Af", "Cf", "Gf", "Tf"), rev = c("Ar", "Cr", "Gr", "Tr")), variable.name = "allele_idx") + long[, allele := c("A", "C", "G", "T")[as.integer(allele_idx)]] + long <- long[allele != ref] + long[, alt := fwd + rev] + long[, vaf := ifelse(depth > 0, alt / depth, 0)] + setorder(long, sample_id, pos, -alt) + top <- long[, .SD[1], by = .(sample_id, pos)] + top[, callable := depth >= opt$`min-depth`] + top[, het := callable & vaf >= opt$`min-vaf` & vaf < opt$`max-vaf` & fwd >= opt$`min-alt-strand` & rev >= opt$`min-alt-strand`] + top[, hom := callable & vaf >= opt$`max-vaf`] + n_samples <- uniqueN(top$sample_id) + hot <- top[vaf >= 0.005 & vaf < opt$`max-vaf` & alt >= 2, .(n_low = .N), by = .(pos, allele)] + hot[, hotspot := n_low / n_samples >= opt$`hotspot-frac`] + top <- merge(top, hot[, .(pos, allele, n_low_samples = n_low, hotspot)], by = c("pos", "allele"), all.x = TRUE) + top[is.na(hotspot), `:=`(hotspot = FALSE, n_low_samples = 0L)] + homs <- unique(top[hom == TRUE, .(pos, allele, run, hom_sample = sample_id)]) + ct <- merge(top[het == TRUE & vaf < opt$`crosstalk-vaf`, .(sample_id, pos, allele, run)], homs, by = c("pos", "allele", "run"), allow.cartesian = TRUE) + ct <- ct[hom_sample != sample_id, .(crosstalk_candidate = TRUE), by = .(sample_id, pos, allele)] + top <- merge(top, ct, by = c("sample_id", "pos", "allele"), all.x = TRUE) + top[is.na(crosstalk_candidate), crosstalk_candidate := FALSE] + top[, class6 := collapse6(ref, allele)] + top[, class12 := paste0(ref, ">", allele)] + top[, clean := het & !hotspot & !crosstalk_candidate] + + sites <- top[het | hom, .(sample_id, group, run, pos, ref, alt_allele = allele, depth, alt_fwd = fwd, alt_rev = rev, vaf = round(vaf, 4), + het, hom, homopolymer, hotspot, n_low_samples, crosstalk_candidate, clean, class12, class6)] + fwrite(sites, file.path(opt$out, "sites_called.tsv"), sep = "\t") + fwrite(hot[hotspot == TRUE][order(pos)], file.path(opt$out, "hotspots.tsv"), sep = "\t") + + callable <- top[, .(callable_pos = sum(callable), mean_depth = mean(depth)), by = .(sample_id, group, run)] + spec <- top[clean == TRUE, .(n_het = .N, n_het_nohp = sum(!homopolymer), vaf_median = median(vaf), + transitions = sum(class6 %in% classes6[1:2]), oxidative = sum(class6 == classes6[3]), + oxidative_nohp = sum(class6 == classes6[3] & !homopolymer)), by = .(sample_id, group, run)] + per_sample <- merge(callable, spec, by = c("sample_id", "group", "run"), all.x = TRUE) + for (k in c("n_het", "n_het_nohp", "transitions", "oxidative", "oxidative_nohp")) per_sample[is.na(get(k)), (k) := 0L] + per_sample[, burden_per_kb := n_het / (callable_pos / 1000)] + per_sample[, transition_frac := ifelse(n_het > 0, transitions / n_het, NA_real_)] + per_sample[, oxidative_frac := ifelse(n_het > 0, oxidative / n_het, NA_real_)] + per_sample[, n_hom := top[hom == TRUE, .N, by = sample_id][match(per_sample$sample_id, sample_id), N]] + per_sample[is.na(n_hom), n_hom := 0L] + mid <- top[clean == TRUE & vaf >= 0.2 & vaf <= 0.8, .(n_het_mid = .N), by = sample_id] + per_sample <- merge(per_sample, mid, by = "sample_id", all.x = TRUE) + per_sample[is.na(n_het_mid), n_het_mid := 0L] + per_sample[, mixture_flag := n_het_mid >= opt$`mixture-sites`] + mixed <- per_sample[mixture_flag == TRUE, sample_id] + if (length(mixed)) cat(sprintf("[call] mixture flag (>= %d clean sites at VAF 0.2-0.8, the signature of DNA from two individuals): %s; excluded from group statistics and spectra\n", opt$`mixture-sites`, paste(mixed, collapse = ", "))) + top[, mixture := sample_id %in% mixed] + fwrite(per_sample[order(group, sample_id)], file.path(opt$out, "spectrum_per_sample.tsv"), sep = "\t") + sites[, mixture_sample := sample_id %in% mixed] + fwrite(sites, file.path(opt$out, "sites_called.tsv"), sep = "\t") + per_sample <- per_sample[mixture_flag == FALSE] + + cls <- top[clean == TRUE & mixture == FALSE, .N, by = .(group, class12)] + cls6 <- top[clean == TRUE & mixture == FALSE, .N, by = .(group, class6)] + cls6w <- dcast(cls6, group ~ class6, value.var = "N", fill = 0) + fwrite(cls, file.path(opt$out, "spectrum12_per_group.tsv"), sep = "\t") + fwrite(cls6w, file.path(opt$out, "spectrum6_per_group.tsv"), sep = "\t") + + tests <- list() + grp <- sort(unique(per_sample$group)); grp <- grp[grp != "unassigned"] + if (length(grp) >= 2) { + ps <- per_sample[group %in% grp]; ps[, group := factor(group, levels = grp)]; ps[, log_depth := log(mean_depth)] + for (v in c("burden_per_kb", "transition_frac", "oxidative_frac", "vaf_median")) { + x <- ps[!is.na(get(v))] + if (nrow(x) < 4 || uniqueN(x$group) < 2) next + p_np <- if (length(grp) == 2) wilcox.test(get(v) ~ group, data = x, exact = FALSE)$p.value else kruskal.test(get(v) ~ group, data = x)$p.value + full <- lm(as.formula(paste(v, "~ group + log_depth")), data = x); red <- lm(as.formula(paste(v, "~ log_depth")), data = x) + p_adj <- anova(red, full)$`Pr(>F)`[2] + meds <- x[, .(median = median(get(v)), n = .N), by = group] + tests[[v]] <- data.table(variable = v, groups = paste(sprintf("%s: median %.3g (n=%d)", meds$group, meds$median, meds$n), collapse = "; "), + p_nonparametric = p_np, p_depth_adjusted = p_adj) + } + m <- as.matrix(cls6w[group %in% grp, -1]); rownames(m) <- cls6w[group %in% grp, group] + if (nrow(m) >= 2 && sum(m) > 0) { + ct_p <- tryCatch(suppressWarnings(chisq.test(m[, colSums(m) > 0, drop = FALSE], simulate.p.value = TRUE, B = 20000)$p.value), error = function(e) NA_real_) + tests[["spectrum6_pooled"]] <- data.table(variable = "spectrum6_pooled", groups = paste(apply(m, 1, function(r) paste(r, collapse = "/")), collapse = " vs "), + p_nonparametric = ct_p, p_depth_adjusted = NA_real_) + } + } + tests_dt <- rbindlist(tests, fill = TRUE) + fwrite(tests_dt, file.path(opt$out, "group_tests.tsv"), sep = "\t") + + sdir <- file.path(opt$out, "sites"); dir.create(sdir, showWarnings = FALSE) + for (s in unique(top$sample_id)) { + b <- top[sample_id == s & het == TRUE & hotspot == FALSE][order(pos), .(chrom = opt$contig, start = pos - 1L, end = pos, name = paste0(ref, ">", allele))] + if (nrow(b) >= 2) fwrite(b, file.path(sdir, paste0(s, ".sites.bed")), sep = "\t", col.names = FALSE) + } + + if (nrow(cls6) > 0) { + cls6[, frac := N / sum(N), by = group] + p1 <- ggplot(cls6, aes(class6, frac, fill = group)) + geom_col(position = position_dodge(width = 0.8)) + + geom_text(aes(label = N, y = frac), position = position_dodge(width = 0.8), vjust = -0.3, size = 2.6) + + labs(title = sprintf("%s: substitution spectrum of clean heteroplasmic sites (VAF >= %g, both strands, no hotspot or cross-talk)", opt$cohort, opt$`min-vaf`), + x = NULL, y = "fraction of sites in group", fill = NULL) + theme_bw(base_size = 11) + theme(legend.position = "top") + ggsave(file.path(opt$out, "spectrum6_by_group.png"), p1, width = 9, height = 5, dpi = 160, bg = "white") + } + het_sites <- top[clean == TRUE & mixture == FALSE] + if (nrow(het_sites) > 0) { + p2 <- ggplot(het_sites, aes(vaf, fill = group)) + geom_histogram(bins = 40, position = "identity", alpha = 0.55) + scale_x_log10() + + labs(title = sprintf("%s: minor allele fraction of clean heteroplasmic sites", opt$cohort), x = "VAF (log)", y = "sites", fill = NULL) + + theme_bw(base_size = 11) + theme(legend.position = "top") + ggsave(file.path(opt$out, "vaf_distribution.png"), p2, width = 8, height = 4.5, dpi = 160, bg = "white") + } + p3 <- ggplot(merge(callable, spec, by = c("sample_id", "group", "run"), all.x = TRUE)[is.na(n_het), n_het := 0L][, burden_per_kb := n_het / (callable_pos / 1000)][, mixture := sample_id %in% mixed], aes(mean_depth, burden_per_kb, colour = group, shape = mixture, label = sample_id)) + geom_point(size = 2.2) + geom_text(vjust = -0.7, size = 2.3) + + scale_x_log10() + labs(title = sprintf("%s: heteroplasmic sites per callable kb against mean depth", opt$cohort), x = "mean depth (log)", y = "clean heteroplasmic sites per kb", colour = NULL) + + theme_bw(base_size = 11) + theme(legend.position = "top") + ggsave(file.path(opt$out, "burden_vs_depth.png"), p3, width = 8, height = 5, dpi = 160, bg = "white") + pos_plot <- top[het == TRUE][, .(n = .N), by = .(pos, hotspot)] + p4 <- ggplot(pos_plot, aes(pos, n, colour = hotspot)) + geom_segment(aes(xend = pos, yend = 0)) + + labs(title = sprintf("%s: heteroplasmic calls per position across samples (hotspots = systematic minor alleles)", opt$cohort), x = sprintf("%s position", opt$contig), y = "samples with a call", colour = "hotspot") + + theme_bw(base_size = 11) + theme(legend.position = "top") + ggsave(file.path(opt$out, "calls_along_genome.png"), p4, width = 12, height = 4, dpi = 160, bg = "white") + + cat(sprintf("[call] %d samples, %d callable positions per sample (median), %d heteroplasmic calls, %d clean (%d hotspot positions, %d cross-talk candidates), %d homoplasmic\n", + n_samples, as.integer(median(callable$callable_pos)), sum(top$het), sum(top$clean), sum(hot$hotspot), sum(top$het & top$crosstalk_candidate), sum(top$hom))) + print(per_sample[, .(sample_id, group, mean_depth = round(mean_depth, 1), callable_pos, n_het, burden_per_kb = round(burden_per_kb, 2), transition_frac = round(transition_frac, 2), oxidative_frac = round(oxidative_frac, 2), n_hom)], nrows = 100) + print(cls6w) + if (nrow(tests_dt)) print(tests_dt) +} + +if (opt$stage == "phase") { + stopifnot(!is.null(opt$bases)) + sites <- fread(file.path(opt$out, "sites_called.tsv")) + files <- list.files(opt$bases, pattern = "[.]readbases[.]tsv$", full.names = TRUE) + pairs <- list(); summ <- list() + for (f in files) { + s <- sub("[.]readbases[.]tsv$", "", basename(f)) + rb <- fread(f); if (nrow(rb) == 0) next + st <- sites[sample_id == s & het == TRUE & hotspot == FALSE] + if (nrow(st) < 2) next + rb <- merge(rb, st[, .(pos, alt_allele, ref)], by = "pos") + rb <- rb[base %in% c("A", "C", "G", "T")] + rb[, minor := base == alt_allele] + wide <- dcast(rb[, .(qname, pos, minor)], qname ~ pos, value.var = "minor") + ps <- sort(st$pos) + res <- list() + for (i in seq_len(length(ps) - 1L)) for (j in (i + 1L):length(ps)) { + a <- as.character(ps[i]); b <- as.character(ps[j]) + if (!(a %in% names(wide)) || !(b %in% names(wide))) next + w <- wide[!is.na(get(a)) & !is.na(get(b))] + n <- nrow(w); if (n < opt$`min-cover`) next + nab <- sum(w[[a]] & w[[b]]); na <- sum(w[[a]] & !w[[b]]); nb <- sum(!w[[a]] & w[[b]]); n0 <- sum(!w[[a]] & !w[[b]]) + ft <- fisher.test(matrix(c(nab, na, nb, n0), 2)) + res[[length(res) + 1L]] <- data.table(sample_id = s, group = st$group[1], pos_a = ps[i], pos_b = ps[j], distance = ps[j] - ps[i], + n_cover_both = n, n_minor_both = nab, n_minor_a_only = na, n_minor_b_only = nb, n_neither = n0, + expected_both = round((nab + na) * (nab + nb) / n, 2), odds_ratio = unname(ft$estimate), p_fisher = ft$p.value, + frac_a_carrying_b = ifelse(nab + na > 0, nab / (nab + na), NA_real_), frac_b_carrying_a = ifelse(nab + nb > 0, nab / (nab + nb), NA_real_)) + } + if (length(res)) { + r <- rbindlist(res); pairs[[s]] <- r + summ[[s]] <- data.table(sample_id = s, group = st$group[1], n_sites = nrow(st), n_pairs_tested = nrow(r), + n_linked = sum(r$p_fisher < 0.01 & r$odds_ratio > 1 & r$n_minor_both >= 3), + n_exclusive = sum(r$p_fisher < 0.01 & r$odds_ratio < 1), max_frac_shared = max(c(r$frac_a_carrying_b, r$frac_b_carrying_a), na.rm = TRUE)) + } + } + pairs_dt <- rbindlist(pairs); summ_dt <- rbindlist(summ) + fwrite(pairs_dt, file.path(opt$out, "phasing_pairs.tsv"), sep = "\t") + fwrite(summ_dt, file.path(opt$out, "phasing_summary.tsv"), sep = "\t") + if (nrow(pairs_dt)) { + pairs_dt[, log_or := log10(pmin(pmax(odds_ratio, 0.01), 100))] + top_s <- summ_dt[order(-n_pairs_tested)][seq_len(min(12, nrow(summ_dt))), sample_id] + pd <- pairs_dt[sample_id %in% top_s] + p <- ggplot(pd, aes(factor(pos_a), factor(pos_b), fill = log_or)) + geom_tile() + + geom_text(aes(label = n_minor_both), size = 2.2) + facet_wrap(~ sample_id, scales = "free") + + scale_fill_gradient2(low = "#2166ac", mid = "grey92", high = "#b2182b", midpoint = 0, name = "log10 odds ratio") + + labs(title = "Co-occurrence of minor alleles on the same molecules (number = reads carrying both)", x = "site A", y = "site B") + + theme_bw(base_size = 9) + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, size = 6), axis.text.y = element_text(size = 6)) + ggsave(file.path(opt$out, "phasing_pairs.png"), p, width = 13, height = 9, dpi = 150, bg = "white") + } + cat(sprintf("[phase] %d samples with >= 2 sites, %d pairs tested, %d linked pairs (p < 0.01, OR > 1, >= 3 shared reads), %d mutually exclusive pairs\n", + nrow(summ_dt), nrow(pairs_dt), sum(pairs_dt$p_fisher < 0.01 & pairs_dt$odds_ratio > 1 & pairs_dt$n_minor_both >= 3), sum(pairs_dt$p_fisher < 0.01 & pairs_dt$odds_ratio < 1))) + if (nrow(summ_dt)) print(summ_dt[order(-n_linked)], nrows = 100) + if (nrow(pairs_dt)) print(pairs_dt[p_fisher < 0.01][order(p_fisher)][seq_len(min(30, sum(pairs_dt$p_fisher < 0.01)))]) +} diff --git a/src/bin/mt_supplementary_check.sh b/src/bin/mt_supplementary_check.sh new file mode 100755 index 0000000000000000000000000000000000000000..605a1cfc79116b71284e31c098e75aece4bcdad4 --- /dev/null +++ b/src/bin/mt_supplementary_check.sh @@ -0,0 +1,52 @@ +#!/usr/bin/env bash +#' Account for supplementary records on the mitochondrial contig before and after the rescue. +#' +#' A read that crosses the origin of the circular contig (or is split elsewhere) aligns as one +#' primary record plus one or more supplementary records. The rescue extracts primary records +#' only, because a primary record carries the whole read (soft-clipped) and its modification +#' tags; realigning it with minimap2 -Y regenerates the supplementary records on the contig. +#' The only reads the rescue cannot carry are those whose primary record lies outside the +#' candidate regions while a supplementary record sits on the contig (chimeric reads). This +#' tool counts both classes per sample. Runs inside images/debian-nanopore.sif. +#' +#' Environment: SRC_DIR (-Unfiltered.bam) AFTER_BAM_DIR (-Filtered__realigned.bam) +#' CAND_DIR (.candidates.bam) MT_CONTIG OUT_SUPPL [MAPQ=10] +set -euo pipefail +export LC_ALL=C +: "${SRC_DIR:?}" "${AFTER_BAM_DIR:?}" "${CAND_DIR:?}" "${MT_CONTIG:?}" "${OUT_SUPPL:?}" +MAPQ="${MAPQ:-10}" +mkdir -p "$(dirname "$OUT_SUPPL")" +TMP=$(mktemp -d) +trap 'rm -rf "$TMP"' EXIT +printf 'sample\tbefore_primary_all\tbefore_primary_q%s\tbefore_suppl_all\tbefore_suppl_q%s\tbefore_split_reads_q%s\tsuppl_reads_lost_all\tsuppl_reads_lost_q%s\tafter_primary\tafter_suppl\tafter_split_reads\tafter_origin_spanners\n' "$MAPQ" "$MAPQ" "$MAPQ" "$MAPQ" > "$OUT_SUPPL" +for ubam in "$SRC_DIR"/*-Unfiltered.bam; do + id=$(basename "$ubam" -Unfiltered.bam) + abam="$AFTER_BAM_DIR/$id-Filtered_${MT_CONTIG}_realigned.bam" + cand="$CAND_DIR/$id.candidates.bam" + if [ ! -s "$abam" ] || [ ! -s "$cand" ]; then echo "[warn] missing realigned or candidate BAM for $id"; continue; fi + mt_len=$(samtools view -H "$abam" | awk -v c="$MT_CONTIG" -F"\t" '$1=="@SQ"{for(i=2;i<=NF;i++){if($i=="SN:"c)f=1; if($i ~ /^LN:/)l=substr($i,4)} if(f){print l; exit}}') + b_pri_all=$(samtools view -c -F 0x904 "$ubam" "$MT_CONTIG") + b_pri_q=$(samtools view -c -F 0x904 -q "$MAPQ" "$ubam" "$MT_CONTIG") + b_sup_all=$(samtools view -c -f 0x800 -F 0x104 "$ubam" "$MT_CONTIG") + b_sup_q=$(samtools view -c -f 0x800 -F 0x104 -q "$MAPQ" "$ubam" "$MT_CONTIG") + b_split_q=$(samtools view -F 0x104 -q "$MAPQ" "$ubam" "$MT_CONTIG" | cut -f1 | sort | uniq -d | wc -l) + samtools view "$cand" | cut -f1 | sort -u > "$TMP/cand.names" + lost_all=$(samtools view -f 0x800 -F 0x104 "$ubam" "$MT_CONTIG" | cut -f1 | sort -u | comm -23 - "$TMP/cand.names" | wc -l) + lost_q=$(samtools view -f 0x800 -F 0x104 -q "$MAPQ" "$ubam" "$MT_CONTIG" | cut -f1 | sort -u | comm -23 - "$TMP/cand.names" | wc -l) + a_pri=$(samtools view -c -F 0x904 "$abam") + a_sup=$(samtools view -c -f 0x800 -F 0x104 "$abam") + a_split=$(samtools view -F 0x104 "$abam" | cut -f1 | sort | uniq -d | wc -l) + a_origin=$(samtools view -F 0x104 "$abam" | awk -v L="$mt_len" '{ + cig=$6; ref=0 + while (match(cig, /^[0-9]+[MIDNSHP=X]/)) { n=substr(cig,1,RLENGTH-1)+0; op=substr(cig,RLENGTH,1); if (op=="M"||op=="D"||op=="N"||op=="="||op=="X") ref+=n; cig=substr(cig,RLENGTH+1) } + s=$4; e=$4+ref-1 + if (s<=50) head[$1]=1 + if (e>=L-50) tail[$1]=1 } + END { for (q in head) if (q in tail) k++; print k+0 }') + printf '%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n' "$id" "$b_pri_all" "$b_pri_q" "$b_sup_all" "$b_sup_q" "$b_split_q" "$lost_all" "$lost_q" "$a_pri" "$a_sup" "$a_split" "$a_origin" >> "$OUT_SUPPL" + echo "[$id] before on $MT_CONTIG: primary $b_pri_all (MAPQ>=$MAPQ $b_pri_q), supplementary $b_sup_all (MAPQ>=$MAPQ $b_sup_q), split reads kept by step 2 $b_split_q; supplementary-only reads the rescue cannot carry: $lost_all (MAPQ>=$MAPQ $lost_q); after: primary $a_pri, supplementary $a_sup, split reads $a_split, of which origin spanners $a_origin" +done +echo +echo "[totals]" +awk -F"\t" 'NR>1{for(i=2;i<=NF;i++) t[i]+=$i} END{printf "before primary all %d, MAPQ-filtered %d; before supplementary all %d, MAPQ-filtered %d; split reads kept by step 2 %d\n", t[2], t[3], t[4], t[5], t[6]; printf "supplementary-only reads lost by the rescue: %d (all), %d (MAPQ-filtered)\n", t[7], t[8]; printf "after primary %d, supplementary %d, split reads %d, origin spanners %d\n", t[9], t[10], t[11], t[12]}' "$OUT_SUPPL" +echo "[ok] wrote $OUT_SUPPL" diff --git a/src/bin/mt_tail_detail.sh b/src/bin/mt_tail_detail.sh new file mode 100755 index 0000000000000000000000000000000000000000..3fa914eb26f2bd0f163712b98ad7fa41af24417d --- /dev/null +++ b/src/bin/mt_tail_detail.sh @@ -0,0 +1,24 @@ +#!/usr/bin/env bash +#' Print the whole-genome record details of dropped MT-origin reads whose tails aligned to +#' recurrent nuclear loci: primary MAPQ, MT segment length and identity, and every SA entry +#' (target, position, strand, CIGAR, MAPQ, NM) so that a spurious low-complexity hit can be told +#' from a long, confident nuclear alignment. Runs inside images/debian-nanopore.sif. +#' +#' Environment: READS (human_chimera_check_reads.tsv) RUNS_ROOT MT_CONTIG PATTERN (regex on sa_targets) +set -euo pipefail +export LC_ALL=C +: "${READS:?}" "${RUNS_ROOT:?}" "${MT_CONTIG:?}" "${PATTERN:?}" +printf 'run\tsample\tqname\tread_len\tprimary_mapq\tmt_aligned\tmt_identity\tsa_entries(target,pos,strand,cigar_summary,mapq,nm)\n' +awk -F"\t" -v p="$PATTERN" 'NR>1 && $8 ~ p {print $1"\t"$2"\t"$3}' "$READS" | while IFS=$'\t' read -r run s q; do + cand="$RUNS_ROOT/${run}_chrM/qc/candidates/$s.candidates.bam" + samtools view -F 0x904 "$cand" "$MT_CONTIG" | awk -F"\t" -v OFS="\t" -v run="$run" -v s="$s" -v q="$q" '$1==q { + cig=$6; aln=0; nm=0; sa="" + while (match(cig, /^[0-9]+[MIDNSHP=X]/)) { n=substr(cig,1,RLENGTH-1)+0; op=substr(cig,RLENGTH,1); if (op=="M"||op=="I"||op=="D"||op=="="||op=="X") aln+=n; cig=substr(cig,RLENGTH+1) } + for (i=12;i<=NF;i++) { if ($i ~ /^NM:i:/) nm=substr($i,6)+0; if ($i ~ /^SA:Z:/) sa=substr($i,6) } + out=""; n2=split(sa, parts, ";") + for (k=1;k<=n2;k++) { if (parts[k]=="") continue; split(parts[k], f, ",") + c2=f[4]; a2=0; s2=0 + while (match(c2, /^[0-9]+[MIDNSHP=X]/)) { m=substr(c2,1,RLENGTH-1)+0; o=substr(c2,RLENGTH,1); if (o=="M"||o=="I"||o=="D"||o=="="||o=="X") a2+=m; if (o=="S"||o=="H") s2+=m; c2=substr(c2,RLENGTH+1) } + out=out sprintf("%s:%s%s aligned=%d clipped=%d mapq=%s nm=%s; ", f[1], f[2], f[3], a2, s2, f[5], f[6]) } + print run, s, substr(q,1,8), length($10), $5, aln, sprintf("%.4f", 1-nm/aln), out }' +done diff --git a/src/bin/rescue_mt.sh b/src/bin/rescue_mt.sh new file mode 100755 index 0000000000000000000000000000000000000000..2195ccf89922474201060f4e77b36ad3076b55fd --- /dev/null +++ b/src/bin/rescue_mt.sh @@ -0,0 +1,222 @@ +#!/usr/bin/env bash +#' Rebuild the mitochondrial BAM of every sample from its unfiltered whole-genome BAM. +#' +#' Stage 1 extracts the candidates: primary records placed on the mitochondrial contig or on +#' any NUMT locus, whatever their MAPQ. Stage 2 realigns them to the mitochondrial contig alone +#' with their modification tags, keeping the origin of each read in its read group (.mt or +#' .numt), drops every read that leaves MAX_CLIP or more bases unaligned to the contig +#' (summed over its primary and supplementary records: the signature of a nuclear read whose +#' flank has no mitochondrial counterpart, while adapters and barcodes stay below it), and +#' profiles depth in 250-bp windows. The gate then applies the decision rule: when the +#' candidates taken from NUMT loci exceed MAX_NUMT_FRAC times those taken from the +#' mitochondrial contig in any sample, the script stops with exit status 3 unless FORCE=1, so +#' that the realignment diagnostics can be judged before any pileup exists. Stage 3 runs modkit +#' with the mt scope. Every stage reuses its outputs when they are newer than its inputs. +#' Runs inside images/debian-nanopore.sif. +#' Rationale and the decision rule: docs/pipeline_notes.md#numt-rescue +#' +#' Environment: SRC_DIR (holds -Unfiltered.bam and .bai) OUT_DIR REF_MT NUMT_BED MT_CONTIG +#' [NUMT_CORE_BED] [SAMPLES colon-separated, default every *-Unfiltered.bam] [THREADS=8] +#' [NUM_READS=50000] [FILTER_MODE=percentile|fixed] [FILTER_PCT=0.1] [FILTER_THR=0.75] +#' [MAX_NUMT_FRAC=0.01] [MIN_ALN_FRAC=0.8] [MAX_CLIP=300, 0 disables] +#' [STAGE=full|candidates|realign] [FORCE=0] [DEST_R] [CODE_REV] [LABEL] [SORT_MEM=512M] +set -euo pipefail +export LC_ALL=C +: "${SRC_DIR:?}" "${OUT_DIR:?}" "${REF_MT:?}" "${NUMT_BED:?}" "${MT_CONTIG:?}" +NUMT_CORE_BED="${NUMT_CORE_BED:-}" +T="${THREADS:-8}" +NUM_READS="${NUM_READS:-50000}" +FILTER_PCT="${FILTER_PCT:-0.1}" +FILTER_MODE="${FILTER_MODE:-percentile}" +FILTER_THR="${FILTER_THR:-0.75}" +if [ "$FILTER_MODE" = fixed ]; then + PILEUP_FILTER="--filter-threshold $FILTER_THR"; SUMMARY_FILTER="--filter-threshold $FILTER_THR"; FILTER_DESC="fixed pass threshold $FILTER_THR" +else + PILEUP_FILTER="--filter-percentile $FILTER_PCT"; SUMMARY_FILTER="--filter-quantile $FILTER_PCT"; FILTER_DESC="pass threshold at percentile $FILTER_PCT of the read probabilities" +fi +MAX_NUMT_FRAC="${MAX_NUMT_FRAC:-0.01}" +MIN_ALN_FRAC="${MIN_ALN_FRAC:-0.8}" +MAX_CLIP="${MAX_CLIP:-300}" +STAGE="${STAGE:-full}" +FORCE="${FORCE:-0}" +DEST_R="${DEST_R:-}" +CODE_REV="${CODE_REV:-unknown}" +LABEL="${LABEL:-$(basename "$OUT_DIR")}" +SORT_MEM="${SORT_MEM:-512M}" + +[ -s "$NUMT_BED" ] || { echo "[error] NUMT bed missing: $NUMT_BED"; exit 1; } +[ -s "$REF_MT" ] || { echo "[error] mitochondrial reference missing: $REF_MT"; exit 1; } +[ -s "$REF_MT.fai" ] || samtools faidx "$REF_MT" +MT_LEN=$(awk -v c="$MT_CONTIG" '$1==c{print $2}' "$REF_MT.fai") +[ -n "$MT_LEN" ] || { echo "[error] $MT_CONTIG not in $REF_MT"; exit 1; } +if [ -z "${SAMPLES:-}" ]; then + SAMPLES=$(ls "$SRC_DIR"/*-Unfiltered.bam | sed "s#.*/##; s#-Unfiltered.bam\$##" | sort -V | tr "\n" " ") +else + SAMPLES=$(echo "$SAMPLES" | tr ":" " ") +fi +N_SAMPLES=$(echo $SAMPLES | wc -w) +[ "$N_SAMPLES" -gt 0 ] || { echo "[error] no *-Unfiltered.bam in $SRC_DIR"; exit 1; } + +CAND="$OUT_DIR/qc/candidates"; DEPTH="$OUT_DIR/qc/mt_depth"; BAMS="$OUT_DIR/qc/bam_filtering" +MODK="$OUT_DIR/modkit/modkit"; PROB="$OUT_DIR/modkit/modkit_prob"; TMP="$OUT_DIR/tmp" +mkdir -p "$CAND" "$DEPTH" "$BAMS" "$MODK" "$PROB" "$TMP" +REG="$CAND/candidate_regions.bed" +{ printf '%s\t0\t%s\t%s\n' "$MT_CONTIG" "$MT_LEN" "$MT_CONTIG"; awk 'BEGIN{OFS="\t"}{print $1,$2,$3,"numt_"NR}' "$NUMT_BED"; } > "$REG" +echo "[info] host=$(hostname) threads=$T samples($N_SAMPLES)=$SAMPLES" +echo "[info] source=$SRC_DIR output=$OUT_DIR mt=$MT_CONTIG ($MT_LEN bp) numt loci=$(wc -l < "$NUMT_BED") stage=$STAGE force=$FORCE" + +PROV="$OUT_DIR/PROVENANCE.txt" +{ + echo "rescue_mt $LABEL, $(date -u +%FT%TZ) on $(hostname), Slurm job ${SLURM_JOB_ID:-none}, pipeline code $CODE_REV" + samtools --version | head -n1; echo "minimap2 $(minimap2 --version)"; modkit --version + echo "mitochondrial reference: $REF_MT ($MT_CONTIG, $MT_LEN bp)" + echo "NUMT loci: $NUMT_BED" + echo "source BAMs: $SRC_DIR" + echo "candidates: primary records, any MAPQ, on $MT_CONTIG or on any NUMT locus (samtools view -F 0x904 -L)" + echo "realignment: samtools fastq -T MM,ML,MN | minimap2 -y -Y -ax lr:hq --secondary=no , origin kept in RG ID .mt or .numt" + echo "read filter: reads with $MAX_CLIP or more bases unaligned to $MT_CONTIG over primary plus supplementary records are dropped (0 = disabled); dropped reads listed per sample in qc/candidates/.clip_dropped.tsv" + echo "modkit: sample-probs, pileup, summary with --region $MT_CONTIG --num-reads $NUM_READS --allow-non-primary --modified-bases 6mA 5mC 5hmC; $FILTER_DESC (FILTER_MODE=$FILTER_MODE)" + echo "decision rule: cand_from_numt / cand_from_mt > $MAX_NUMT_FRAC in any sample stops before modkit unless FORCE=1 (FORCE=$FORCE)" +} > "$PROV" + +CS="$OUT_DIR/qc/candidate_summary.tsv" +printf 'sample\tcand_total\tcand_from_mt\tcand_from_numt\tcand_numt_core\tnumt_frac\tcand_mapq0\n' > "$CS" +for id in $SAMPLES; do + ubam="$SRC_DIR/$id-Unfiltered.bam" + [ -s "$ubam" ] && [ -s "$ubam.bai" ] || { echo "[error] missing $ubam or its index"; exit 1; } + cand="$CAND/$id.candidates.bam" + if [ -s "$cand" ] && [ -s "$cand.bai" ] && [ "$cand" -nt "$ubam" ]; then + echo "[$id] candidates already extracted, reusing $cand" + else + echo "[$id] extracting candidate reads (primary records, any MAPQ) from $MT_CONTIG and NUMT loci" + samtools view -@ "$T" -b -F 0x904 -L "$REG" "$ubam" -o "$cand" + samtools index "$cand" + fi + n_tot=$(samtools view -c "$cand") + n_mt=$(samtools view -c "$cand" "$MT_CONTIG") + n_numt=$((n_tot - n_mt)) + n_core=0 + if [ -n "$NUMT_CORE_BED" ]; then n_core=$(samtools view -c -L "$NUMT_CORE_BED" "$cand"); fi + n_q0=$((n_tot - $(samtools view -c -q 1 "$cand"))) + frac=$(awk -v a="$n_numt" -v b="$n_mt" 'BEGIN{printf "%.5f", (b>0) ? a/b : 0}') + echo "[$id] candidates $n_tot: on $MT_CONTIG $n_mt, on NUMT loci $n_numt (overlapping a NUMT core $n_core), numt/mt $frac, MAPQ 0: $n_q0" + printf '%s\t%s\t%s\t%s\t%s\t%s\t%s\n' "$id" "$n_tot" "$n_mt" "$n_numt" "$n_core" "$frac" "$n_q0" >> "$CS" +done +if [ "$STAGE" = candidates ]; then echo "[ok] candidates stage complete"; cat "$CS"; exit 0; fi + +RS="$OUT_DIR/qc/realign_summary.tsv" +printf 'sample\tcand_total\tcand_from_mt\tcand_from_numt\tcand_numt_core\tnumt_frac\tcand_mapq0\trealigned_mapped\trealigned_primary\tsupplementary\tfrom_mt\tfrom_numt\tfrom_mt_alnfrac_low\tfrom_numt_alnfrac_low\tclip_dropped_mt\tclip_dropped_numt\tdepth_windows_lt3\tmedian_depth\n' > "$RS" +for id in $SAMPLES; do + cand="$CAND/$id.candidates.bam" + all_bam="$CAND/$id.realigned_all.bam" + out_bam="$BAMS/$id-Filtered_${MT_CONTIG}_realigned.bam" + reads="$CAND/$id.reads.tsv" + st="$CAND/$id.origin_stats.tsv" + if [ -s "$all_bam" ] && [ -s "$all_bam.bai" ] && [ "$all_bam" -nt "$cand" ]; then + echo "[$id] realigned BAM already present, reusing $all_bam" + else + echo "[$id] realigning candidates to $MT_CONTIG alone, by origin" + for origin in mt numt; do + sub="$TMP/$id.$origin.cand.bam" + if [ "$origin" = mt ]; then + samtools view -@ "$T" -b "$cand" "$MT_CONTIG" -o "$sub" + else + samtools view -@ "$T" -b -L "$NUMT_BED" "$cand" -o "$sub" + fi + samtools fastq -@ "$T" -T MM,ML,MN "$sub" 2> "$CAND/$id.$origin.fastq.log" \ + | minimap2 -t "$T" -y -Y -ax lr:hq --secondary=no -R "@RG\tID:$id.$origin\tSM:$id" "$REF_MT" - 2> "$CAND/$id.$origin.minimap2.log" \ + | samtools view -@ "$T" -b -F 4 - \ + | samtools sort -@ "$T" -m "$SORT_MEM" -o "$TMP/$id.$origin.bam" - + done + samtools merge -@ "$T" -f -o "$all_bam" "$TMP/$id.mt.bam" "$TMP/$id.numt.bam" + samtools index "$all_bam" + rm -f "$TMP/$id."* + fi + samtools view -F 0x104 "$all_bam" | awk -F"\t" 'BEGIN{OFS="\t"} + { rg="NA"; for (i=12;i<=NF;i++) if ($i ~ /^RG:Z:/) { rg=substr($i,6); sub(/^.*\./,"",rg); break } + cig=$6; aln=0; clip=0 + while (match(cig, /^[0-9]+[MIDNSHP=X]/)) { + n=substr(cig,1,RLENGTH-1)+0; op=substr(cig,RLENGTH,1) + if (op=="M"||op=="I"||op=="="||op=="X") aln+=n; else if (op=="S"||op=="H") clip+=n + cig=substr(cig,RLENGTH+1) } + A[$1]+=aln + if (int($2/2048)%2==0) { L[$1]=aln+clip; R[$1]=rg } } + END { for (q in L) if (L[q]>0) print q, R[q], L[q], sprintf("%.4f", A[q]/L[q]), L[q]-A[q] }' > "$reads" + awk -F"\t" -v minf="$MIN_ALN_FRAC" -v mc="$MAX_CLIP" 'BEGIN{OFS="\t"} { tot[$2]++; if ($4 < minf) low[$2]++; if (mc > 0 && $5 >= mc) drop[$2]++ } + END { for (r in tot) print r, tot[r], low[r]+0, drop[r]+0 }' "$reads" > "$st" + if [ "$MAX_CLIP" -gt 0 ]; then + awk -F"\t" -v mc="$MAX_CLIP" '$5 >= mc' "$reads" > "$CAND/$id.clip_dropped.tsv" + awk -F"\t" -v mc="$MAX_CLIP" '$5 < mc {print $1}' "$reads" > "$TMP/$id.keep.txt" + samtools view -@ "$T" -b -N "$TMP/$id.keep.txt" -o "$out_bam" "$all_bam" + rm -f "$TMP/$id.keep.txt" + else + : > "$CAND/$id.clip_dropped.tsv" + cp "$all_bam" "$out_bam" + fi + samtools index "$out_bam" + samtools flagstat -@ "$T" "$out_bam" > "$BAMS/$id-Filtered_${MT_CONTIG}_realigned.flagstat" + n_map=$(samtools view -c "$out_bam") + n_pri=$(samtools view -c -F 0x900 "$out_bam") + n_sup=$(samtools view -c -f 0x800 "$out_bam") + from_mt=$(awk -F"\t" '$1=="mt"{print $2}' "$st"); from_mt=${from_mt:-0} + from_numt=$(awk -F"\t" '$1=="numt"{print $2}' "$st"); from_numt=${from_numt:-0} + low_mt=$(awk -F"\t" '$1=="mt"{print $3}' "$st"); low_mt=${low_mt:-0} + low_numt=$(awk -F"\t" '$1=="numt"{print $3}' "$st"); low_numt=${low_numt:-0} + drop_mt=$(awk -F"\t" '$1=="mt"{print $4}' "$st"); drop_mt=${drop_mt:-0} + drop_numt=$(awk -F"\t" '$1=="numt"{print $4}' "$st"); drop_numt=${drop_numt:-0} + echo "[$id] realigned: primaries from $MT_CONTIG $from_mt (aligned fraction < $MIN_ALN_FRAC: $low_mt, dropped for $MAX_CLIP+ unaligned bases: $drop_mt), from NUMT loci $from_numt (aligned fraction < $MIN_ALN_FRAC: $low_numt, dropped: $drop_numt); kept $n_map records ($n_pri primary, $n_sup supplementary)" + dep="$DEPTH/$id.depth250.tsv" + samtools depth -a -r "$MT_CONTIG" "$out_bam" \ + | awk 'BEGIN{OFS="\t"}{b=int(($2-1)/250); s[b]+=$3; n[b]++} END{for(b in s) print b*250+1, s[b]/n[b]}' | sort -k1,1n > "$dep" + lt3=$(awk '$2<3' "$dep" | wc -l) + med=$(sort -k2,2n "$dep" | awk '{a[NR]=$2} END{print a[int((NR+1)/2)]}') + echo "[$id] 250-bp windows below 3 reads: $lt3 of $(wc -l < "$dep"); median window depth $med" + row=$(awk -F"\t" -v s="$id" '$1==s' "$CS" | cut -f2-) + printf '%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n' "$id" "$row" "$n_map" "$n_pri" "$n_sup" "$from_mt" "$from_numt" "$low_mt" "$low_numt" "$drop_mt" "$drop_numt" "$lt3" "$med" >> "$RS" +done +echo +echo "[realignment summary]" +cat "$RS" + +maxfrac=$(awk -F"\t" 'NR>1{if($6>m)m=$6} END{printf "%.5f", m+0}' "$CS") +over=$(awk -F"\t" -v t="$MAX_NUMT_FRAC" 'NR>1 && $6>t{n++} END{print n+0}' "$CS") +echo "[gate] highest numt/mt candidate fraction $maxfrac; $over of $N_SAMPLES samples above the limit $MAX_NUMT_FRAC" +if [ "$over" -gt 0 ] && [ "$FORCE" != 1 ]; then + echo "[decision] pileups are not produced automatically for this project: $over sample(s) take more than $MAX_NUMT_FRAC of their candidates from NUMT loci." + echo "[decision] judge $RS (from_numt_alnfrac_low counts reads from NUMT loci that align only partially to $MT_CONTIG, the nuclear signature), choose per docs/pipeline_notes.md#numt-rescue, then rerun with FORCE=1 or with the competitive reference." + exit 3 +fi +if [ "$STAGE" = realign ]; then echo "[ok] realignment stage complete"; exit 0; fi + +SUMMARY="$OUT_DIR/qc/rescue_summary.tsv" +{ head -1 "$RS" | tr -d "\n"; printf '\tthresholds\n'; } > "$SUMMARY" +for id in $SAMPLES; do + out_bam="$BAMS/$id-Filtered_${MT_CONTIG}_realigned.bam" + echo "[$id] modkit on $MT_CONTIG (thresholds from every read of the contig)" + pfx="$MODK/${id}_modkit" + modkit sample-probs --force --hist --threads "$T" --interval-size 10000 --region "$MT_CONTIG" --num-reads "$NUM_READS" \ + --out-dir "$PROB/$id" "$out_bam" > "$PROB/$id.sample_probs.log" 2>&1 + modkit pileup --threads "$T" --interval-size 10000 $PILEUP_FILTER --modified-bases 6mA 5mC 5hmC \ + --reference "$REF_MT" --region "$MT_CONTIG" --num-reads "$NUM_READS" --allow-non-primary --bgzf \ + "$out_bam" "${pfx}_pileup.bed.gz" --log-filepath "${pfx}_pileup.log" + tabix -f -p bed "${pfx}_pileup.bed.gz" + modkit summary --threads "$T" --interval-size 10000 $SUMMARY_FILTER --region "$MT_CONTIG" --num-reads "$NUM_READS" \ + "$out_bam" > "${pfx}_summary.txt" 2> "${pfx}_summary.log" + thr=$(grep -E "pass_threshold" "${pfx}_summary.txt" | tr -s " " | tr "\n" ";") + echo "[$id] modkit done; $thr" + { awk -F"\t" -v s="$id" '$1==s' "$RS" | tr -d "\n"; printf '\t%s\n' "$thr"; } >> "$SUMMARY" +done + +echo +echo "[summary]" +cat "$SUMMARY" +N=$(ls "$MODK"/*_modkit_pileup.bed.gz | wc -l) +[ "$N" -eq "$N_SAMPLES" ] || { echo "[error] expected $N_SAMPLES pileups in $MODK, found $N"; exit 1; } +if [ -n "$DEST_R" ]; then + mkdir -p "$DEST_R" + cp "$MODK"/*_modkit_pileup.bed.gz "$MODK"/*_modkit_pileup.bed.gz.tbi "$MODK"/*_modkit_pileup.log "$MODK"/*_modkit_summary.txt "$DEST_R/" + cp "$SUMMARY" "$DEST_R/rescue_summary.tsv"; cp "$PROV" "$DEST_R/PROVENANCE.txt" +fi +rm -rf "$TMP" +du -sh "$OUT_DIR" +echo "[ok] mitochondrial rescue complete for $N_SAMPLES samples in $OUT_DIR" diff --git a/src/main.nf b/src/main.nf index 21d49a69b0216f62e20f23231acd47f9ca7fb978..fb514d73eaf58b485c3c8d6ca49cdc542797f695 100644 --- a/src/main.nf +++ b/src/main.nf @@ -51,10 +51,14 @@ workflow { min mapped reads per sample/barcode : ${params.min_mapped_reads_thresh} BAMs are barcoded : ${params.is_barcoded} samtools threads : ${params.samtools_threads} + mitochondrial rescue (rescue_mt) : ${params.rescue_mt} (NUMT loci: ${params.numt_bed}) expected output dirs : bam_filtering, intermediate_qc_reports, multiqc_input QC output directory : ${params.qc_out_dir} ====================================== """ + if (params.rescue_mt && !params.numt_bed) { + error "rescue_mt requires numt_bed (references/numt//numt_loci_.bed); see docs/pipeline_notes.md#numt_rescue" + } } else if (params.step == 3) { log.info """ =============================================== @@ -69,6 +73,7 @@ workflow { reference fasta : ${params.reference_file} samtools threads : ${params.samtools_threads} expected inputs : bam_filtering/*-Filtered*.bam + intermediate_qc_reports + MT-only BAM source : ${params.rescue_mt ? 'bam_filtering_mt/*.mt.bam (step-2 rescue)' : 'EXTRACT_MT from the filtered BAM'} modkit output directory : ${params.modkit_out_dir} mtDNA contig : ${params.mt_contig} =============================================== @@ -119,10 +124,14 @@ workflow { .flatten() mapq_ch = channel.value(params.mapq) qscore_thresh_ch = channel.value(params.qscore_thresh) + mt_reference_ch = file("${file(params.reference_file).parent}/${params.mt_contig}.fa") + numt_bed_ch = file(params.numt_bed ?: "${params.qc_out_dir}/numt_bed_not_set") } else if (params.step == "2_from_minknow") { input_dir = channel.fromPath("${params.steps_2_input_directory}/") mapq_ch = channel.value(params.mapq) qscore_thresh_ch = channel.value(params.qscore_thresh) + mt_reference_ch = file("${file(params.reference_file).parent}/${params.mt_contig}.fa") + numt_bed_ch = file(params.numt_bed ?: "${params.qc_out_dir}/numt_bed_not_set") } else if (params.step == 3) { filtered_bams = channel.fromPath("${params.steps_3_input_directory}/bam_filtering/*-Filtered*.bam") .map { file -> tuple(file.baseName, file) } @@ -133,7 +142,7 @@ workflow { if (id.endsWith('.bam')) id = id[0..-5] tuple(id, bai) } - // join BAM+BAI by id, not emission order, or pairs misalign (see docs/pipeline_notes.md#main) + // see docs/pipeline_notes.md#main paired_bam_bai = filtered_bams.join(filtered_bais, by: 0) filtered_bams = paired_bam_bai.map { id, bam, bai -> tuple(id, bam) } .toSortedList { a, b -> a[0] <=> b[0] } @@ -154,6 +163,11 @@ workflow { modkit_qsize_ch = channel.value(params.modkit_extract_qsize) samplesheet_ch = file(params.samplesheet) mtdna_metrics_files = channel.fromPath("${params.steps_3_input_directory}/intermediate_qc_reports/mtdna_metrics/*") + // see docs/pipeline_notes.md#numt_rescue + rescued_mt_bams = params.rescue_mt + ? channel.fromPath("${params.steps_3_input_directory}/bam_filtering_mt/*.mt.bam").map { f -> tuple(f.simpleName, f) } + .join(channel.fromPath("${params.steps_3_input_directory}/bam_filtering_mt/*.mt.bam.bai").map { f -> tuple(f.simpleName, f) }, by: 0) + : channel.empty() } if (params.step == 1) { INDEX_REFERENCE(reference_file) @@ -178,14 +192,18 @@ workflow { txt_files, mapq_ch, qscore_thresh_ch, - samtools_threads_ch + samtools_threads_ch, + mt_reference_ch, + numt_bed_ch ) } else if (params.step == "2_from_minknow") { FILTERING_AND_QC_FROM_MINKNOW( input_dir, mapq_ch, qscore_thresh_ch, - samtools_threads_ch + samtools_threads_ch, + mt_reference_ch, + numt_bed_ch ) } else if (params.step == 3) { MODKIT_AND_MULTIQC( @@ -202,7 +220,8 @@ workflow { modkit_isize_ch, modkit_qsize_ch, samplesheet_ch, - mtdna_metrics_files + mtdna_metrics_files, + rescued_mt_bams ) } } \ No newline at end of file diff --git a/src/modules/basecall.nf b/src/modules/basecall.nf index 24fcb90a6a5efc1ce8b646a4e2452235fae6eb56..6de99b8d87c800df974a7b8d6e999eaf63e0b26e 100644 --- a/src/modules/basecall.nf +++ b/src/modules/basecall.nf @@ -40,7 +40,7 @@ process FAST5_to_POD5 { } process BASECALL { - // No publishDir: adding one re-duplicates every demux BAM in results/ (see docs/pipeline_notes.md#basecall) + // see docs/pipeline_notes.md#basecall label 'gpu' input: @@ -104,11 +104,11 @@ process BASECALL { echo "Processing \$file..." - # trim to temp first so header is fully written before aligner reads it (see docs/pipeline_notes.md#basecall) + # see docs/pipeline_notes.md#basecall dorado trim "\$file" --sequencing-kit "${barcoding_kit}" > "trimmed_temp.bam" if [ -s "trimmed_temp.bam" ]; then - # -Y soft-clips supplementary so they keep seq+MM/ML/MN for modkit --allow-non-primary; hard-clip breaks D-loop (see docs/pipeline_notes.md#basecall) + # see docs/pipeline_notes.md#basecall dorado aligner "\$REF_PATH" "trimmed_temp.bam" --mm2-opts "-Y" \ | samtools sort -@ ${samtools_threads} -o "../\$file" diff --git a/src/modules/extract_mt.nf b/src/modules/extract_mt.nf new file mode 100644 index 0000000000000000000000000000000000000000..d6199e8aa14b35c76cbeb7a68ff1263ddccdfac4 --- /dev/null +++ b/src/modules/extract_mt.nf @@ -0,0 +1,27 @@ +// see docs/pipeline_notes.md#extract_mt +process EXTRACT_MT { + publishDir "${params.modkit_out_dir}/mt_bams/", mode: "copy", overwrite: true + label 'cpu' + errorStrategy 'ignore' + + input: + tuple val(id), path(bam), path(bai) + + output: + tuple val(id), path("${id}.mt.bam"), path("${id}.mt.bam.bai"), emit: mt_bam, optional: true + + script: + // see docs/pipeline_notes.md#extract_mt + def supp = params.mt_keep_supplementary ? "" : "-F 0x800" + """ + set -euo pipefail + samtools view -b ${supp} "${bam}" "${params.mt_contig}" > mt_raw.bam + # see docs/pipeline_notes.md#extract_mt + n=\$(samtools view -c mt_raw.bam) + [ "\$n" -gt 0 ] || { echo "[EXTRACT_MT] 0 reads on ${params.mt_contig} for ${id} -- check params.mt_contig vs the BAM contig name" >&2; exit 1; } + samtools addreplacerg -r 'ID:${id}' -r 'SM:${id}' mt_raw.bam -o mt_rg.bam + samtools sort -@ ${params.samtools_threads} mt_rg.bam -o "${id}.mt.bam" + samtools index "${id}.mt.bam" + rm -f mt_raw.bam mt_rg.bam + """ +} diff --git a/src/modules/filter_bam.nf b/src/modules/filter_bam.nf index 579813488d8d00a2c17dccc13ee74e17838f7066..6662867e86e6662834b55b59a72c288d8928e4c8 100644 --- a/src/modules/filter_bam.nf +++ b/src/modules/filter_bam.nf @@ -27,7 +27,7 @@ process FILTER_BAM { samtools index -@ ${samtools_threads} "${id}-Unfiltered.bam" samtools flagstat -@ ${samtools_threads} "${id}-Unfiltered.bam" > "${id}-Unfiltered.flagstat" samtools idxstats -@ ${samtools_threads} "${id}-Unfiltered.bam" > "${id}-Unfiltered.idxstat" - # -F 0x100 keeps supplementary (drops only secondary) for F7-A recovery via modkit --allow-non-primary; see docs/pipeline_notes.md#filter_bam + # see docs/pipeline_notes.md#filter_bam samtools view -@ ${samtools_threads} -b -q ${mapq} -F 0x100 ${bam} > "intermediate.bam" samtools sort -@ ${samtools_threads} "intermediate.bam" -o "${id}-Filtered_primary_mapq_${mapq}.bam" samtools index -@ ${samtools_threads} "${id}-Filtered_primary_mapq_${mapq}.bam" diff --git a/src/modules/haplogroup.nf b/src/modules/haplogroup.nf index ef5a2c5230654f31ac468369093215dab571ef63..29b4a525e25f74a96d7376d86042e0b530f56ce6 100644 --- a/src/modules/haplogroup.nf +++ b/src/modules/haplogroup.nf @@ -2,7 +2,8 @@ process CLAIR3 { container "${params.clair3_sif}" publishDir "${params.modkit_out_dir}/haplogroup/clair3/", mode: "copy", overwrite: true label 'cpu' - errorStrategy 'ignore' // low coverage/no variants must not sink step 3 (see docs/pipeline_notes.md#haplogroup) + // see docs/pipeline_notes.md#haplogroup + errorStrategy 'ignore' input: tuple val(id), path(filtered_bam), path(filtered_bai) @@ -17,7 +18,7 @@ process CLAIR3 { """ set -euo pipefail samtools faidx ${reference} - # mtDNA is HAPLOID (--haploid_precise); D-loop sits at the rCRS head/tail (--enable_..._head_and_tail), which Clair3 skips by default (see docs/pipeline_notes.md#haplogroup) + # see docs/pipeline_notes.md#haplogroup run_clair3.sh \\ --bam_fn="${filtered_bam}" \\ --ref_fn="${reference}" \\ @@ -50,13 +51,12 @@ process HAPLOGROUP { script: """ set -uo pipefail - # always exit 0 and write a PUBLISHED log: a failed task skips publishDir, which is how F6 originally vanished silently; the .txt exists only on success (see docs/pipeline_notes.md#haplogroup) + # see docs/pipeline_notes.md#haplogroup { echo "[HAPLOGROUP] id=${id} tree=${tree} vcf=${vcf}" if ! command -v haplogrep3 >/dev/null 2>&1; then echo "[FATAL] haplogrep3 not on PATH -- rebuild containers/debian-nanopore.def" else - # set -uo (NOT -e) + no --version probe: haplogrep3 has no --version (errors exit 2); set -e aborted before classify -> run01 F6 failure (see docs/pipeline_notes.md#haplogroup) if haplogrep3 classify --tree ${tree} --in "${vcf}" --out "${id}_haplogroup.txt"; then echo "[HAPLOGROUP] OK -> ${id}_haplogroup.txt" else @@ -84,7 +84,7 @@ process MERGE_HAPLOGROUPS { script: """ set -euo pipefail - # join rank-1 rows (\$3==1) to the samplesheet by real_barcode; strip Haplogrep3's quoted fields first; SampleID==real_barcode; samplesheet cols 1=run_id 3=real_barcode 4=subject_id 7=group (see docs/pipeline_notes.md#haplogroup) + # see docs/pipeline_notes.md#haplogroup awk 'BEGIN{ FS="\\t"; OFS="\\t" } FNR==NR { if (FNR==1) next diff --git a/src/modules/heteroplasmy.nf b/src/modules/heteroplasmy.nf new file mode 100644 index 0000000000000000000000000000000000000000..90d88a335cc07aa575b2ada9db2678c71df97d9e --- /dev/null +++ b/src/modules/heteroplasmy.nf @@ -0,0 +1,164 @@ +// see docs/pipeline_notes.md#heteroplasmy + + +process DOWNSAMPLE_MATCH { + publishDir "${params.modkit_out_dir}/heteroplasmy/matched_bams/", mode: "copy", pattern: "*.matched.bam*", overwrite: true + publishDir "${params.modkit_out_dir}/heteroplasmy/logs/", mode: "copy", pattern: "*.downsample.log", overwrite: true + label 'cpu' + errorStrategy 'ignore' + + input: + tuple val(id), path(mt_bam), path(mt_bai) + + output: + tuple val(id), path("${id}.matched.bam"), path("${id}.matched.bam.bai"), emit: matched, optional: true + path "${id}.downsample.log", emit: log + + script: + """ + set -euo pipefail + T=${params.het_downsample_target} + SEED=${params.het_downsample_seed} + # see docs/pipeline_notes.md#heteroplasmy + mean=\$(samtools coverage -r "${params.mt_contig}" "${mt_bam}" | awk 'NR==2 {print \$7}') + mean=\${mean:-0} + if awk -v m="\$mean" -v t="\$T" 'BEGIN{ exit !(m+0 >= t+0 && m+0 > 0) }'; then + frac=\$(awk -v c="\$T" -v m="\$mean" 'BEGIN{printf "%.4f", c/m}') + samtools view --subsample-seed \$SEED --subsample \$frac -b -o "${id}.matched.bam" "${mt_bam}" + samtools index "${id}.matched.bam" + printf '%s\\tmean=%s\\tT=%s\\tfraction=%s\\tstatus=downsampled\\n' "${id}" "\$mean" "\$T" "\$frac" > "${id}.downsample.log" + else + printf '%s\\tmean=%s\\tT=%s\\tfraction=NA\\tstatus=excluded_below_target\\n' "${id}" "\$mean" "\$T" > "${id}.downsample.log" + fi + """ +} + +process MUTSERVE { + container "${params.mtdnaserver_sif}" + publishDir "${params.modkit_out_dir}/heteroplasmy/mutserve/${arm}/", mode: "copy", overwrite: true + label 'cpu' + errorStrategy 'ignore' + + input: + tuple val(id), path(bam), path(bai) + path reference + val arm + + output: + tuple val(id), val(arm), path("${id}_${arm}.txt"), emit: table, optional: true + tuple val(id), val(arm), path("${id}_${arm}.vcf.gz"), path("${id}_${arm}.vcf.gz.tbi"), emit: vcf, optional: true + path "*_raw.txt", emit: raw, optional: true + + script: + // see docs/pipeline_notes.md#heteroplasmy + def avail_mem = task.memory ? (task.memory.mega * 0.8).intValue() : 1024 + """ + set -euo pipefail + samtools faidx ${reference} + samtools index ${bam} + java -Xmx${avail_mem}M -jar /opt/mutserve/mutserve.jar call \\ + --level ${params.het_detection_level} \\ + --reference ${reference} \\ + --mapQ ${params.het_mapq} \\ + --baseQ ${params.het_baseq} \\ + --output ${id}_${arm}.vcf.gz \\ + --no-ansi \\ + --strand-bias ${params.het_strand_bias} \\ + --write-raw \\ + ${bam} + bcftools norm -m-any -f ${reference} -o ${id}_${arm}.norm.vcf.gz -Oz ${id}_${arm}.vcf.gz + mv ${id}_${arm}.norm.vcf.gz ${id}_${arm}.vcf.gz + tabix -f ${id}_${arm}.vcf.gz + """ +} + +process PARSE_HETEROPLASMY { + label 'cpu' + errorStrategy 'ignore' + + input: + tuple val(id), val(arm), path(table) + + output: + path "${id}_${arm}.variants.tsv", emit: tsv, optional: true + + script: + // see docs/pipeline_notes.md#heteroplasmy + """ + set -euo pipefail + awk -F'\\t' -v OFS='\\t' -v id='${id}' -v arm='${arm}' 'NR>1 && \$2=="PASS" { + print id, arm, \$3, \$4, \$5, \$6, \$11, \$12, \$13, \$14 + }' "${table}" > "${id}_${arm}.variants.tsv" + """ +} + +process MERGE_HETEROPLASMY { + publishDir "${params.modkit_out_dir}/heteroplasmy/", mode: "copy", overwrite: true + label 'cpu' + stageInMode 'copy' + + input: + path tsvs + path samplesheet + + output: + path "heteroplasmy_variants_master.tsv", emit: master + + script: + // see docs/pipeline_notes.md#heteroplasmy + """ + set -euo pipefail + printf 'sample_id\\tarm\\treal_barcode\\trun_id\\tgroup\\tpos\\tref\\tvariant\\tvariant_level\\tcoverage\\tcov_fwd\\tcov_rev\\ttype\\n' > heteroplasmy_variants_master.tsv + awk 'BEGIN{ FS="\\t"; OFS="\\t" } + FNR==NR { if (FNR==1) next; split(\$0,c,","); run[c[3]]=c[1]; grp[c[3]]=c[7]; next } + { split(\$1,p,"-"); rb=p[1]; r=(rb in run)?run[rb]:"NA"; g=(rb in grp)?grp[rb]:"NA"; + print \$1, \$2, rb, r, g, \$3, \$4, \$5, \$6, \$7, \$8, \$9, \$10 }' \\ + "${samplesheet}" ${tsvs} >> heteroplasmy_variants_master.tsv + """ +} + +process MERGE_VCFS_HET { + container "${params.mtdnaserver_sif}" + label 'cpu' + errorStrategy 'ignore' + + input: + path vcfs + path tbis + + output: + tuple path("heteroplasmy_native_merged.vcf.gz"), path("heteroplasmy_native_merged.vcf.gz.tbi"), emit: merged, optional: true + + script: + """ + set -euo pipefail + # see docs/pipeline_notes.md#heteroplasmy + set -- ${vcfs} + if [ "\$#" -ge 2 ]; then + bcftools merge -m none -Oz -o heteroplasmy_native_merged.vcf.gz ${vcfs} + else + cp "\$1" heteroplasmy_native_merged.vcf.gz + fi + tabix -f heteroplasmy_native_merged.vcf.gz + """ +} + +process HAPLOCHECK_HET { + container "${params.mtdnaserver_sif}" + publishDir "${params.modkit_out_dir}/heteroplasmy/", mode: "copy", overwrite: true + label 'cpu' + errorStrategy 'ignore' + + input: + tuple path(merged_vcf), path(merged_tbi) + + output: + path "haplocheck*.txt", emit: report, optional: true + + script: + // see docs/pipeline_notes.md#heteroplasmy + """ + set -uo pipefail + java -jar /opt/haplocheck/haplocheck.jar --out haplocheck.txt --raw "${merged_vcf}" + """ +} diff --git a/src/modules/modkit.nf b/src/modules/modkit.nf index e997166784d9fa29ad2334e2e05e23dfcef004b8..fd9998fde83966ac7a78bb5ab5c2e0ffe7b0fd67 100644 --- a/src/modules/modkit.nf +++ b/src/modules/modkit.nf @@ -3,8 +3,7 @@ process MODKIT { label 'cpu' input: - tuple val(id), path(bam) - path bai + tuple val(id), path(bam), path(bai) path reference_file val modkit_threads val modkit_isize @@ -15,16 +14,21 @@ process MODKIT { path "*", emit: allfiles script: + def fixed_thr = params.modkit_filter_mode == 'fixed' + def filter_arg = fixed_thr ? "--filter-threshold ${params.modkit_filter_threshold}" : "--filter-percentile ${params.modkit_filter_percentile}" + def summary_arg = fixed_thr ? "--filter-threshold ${params.modkit_filter_threshold}" : "--filter-quantile ${params.modkit_filter_percentile}" + // see docs/pipeline_notes.md#modkit """ samtools faidx "${reference_file}" echo "[$id] starting base modification probability sampling (sample-probs)" - # modkit 0.6.3: sample-probs has NO --mapped-only flag; do not add it (see docs/pipeline_notes.md#modkit) + # see docs/pipeline_notes.md#modkit modkit sample-probs \ --hist \ --threads ${modkit_threads} \ --interval-size ${modkit_isize} \ --region ${params.mt_contig} \ + --num-reads ${params.modkit_num_reads} \ --out-dir "./modkit_prob/" "${bam}" echo "[$id] sample-probs completed" @@ -35,7 +39,7 @@ process MODKIT { --threads ${modkit_threads} \ --interval-size ${modkit_isize} \ --reference "${reference_file}" \ - --filter-threshold 0.75 \ + ${filter_arg} \ --region ${params.mt_contig} \ --mapped-only \ --allow-non-primary \ @@ -43,14 +47,15 @@ process MODKIT { echo "[$id] extract calls completed" echo "[$id] starting modkit pileup (site-level bedMethyl, bgzf + tabix for dmr/localize/motif)" - # do NOT add --header to pileup: it space-delimits cols >10 and breaks tabix/dmr (see docs/pipeline_notes.md#modkit) + # see docs/pipeline_notes.md#modkit modkit pileup \ --threads ${modkit_threads} \ --interval-size ${modkit_isize} \ - --filter-threshold 0.75 \ + ${filter_arg} \ --modified-bases 6mA 5mC 5hmC \ --reference "${reference_file}" \ --region ${params.mt_contig} \ + --num-reads ${params.modkit_num_reads} \ --allow-non-primary \ --bgzf \ "${bam}" "${id}_modkit_pileup.bed.gz" \ @@ -62,8 +67,9 @@ process MODKIT { modkit summary \ --threads ${modkit_threads} \ --interval-size ${modkit_isize} \ - --filter-threshold 0.75 \ + ${summary_arg} \ --region ${params.mt_contig} \ + --num-reads ${params.modkit_num_reads} \ ${bam} > "${id}_modkit_summary.txt" echo "[$id] summary successful" """ diff --git a/src/modules/modkit_dmr.nf b/src/modules/modkit_dmr.nf index 26799b07b1021965f06320ddfe6994bbf593cdd9..4f92760bf5355940974c64e86b25b97207ab94a8 100644 --- a/src/modules/modkit_dmr.nf +++ b/src/modules/modkit_dmr.nf @@ -1,7 +1,7 @@ process MODKIT_DMR { publishDir "${params.modkit_out_dir}/dmr/", mode: "copy", overwrite: true label 'cpu' - // errorStrategy 'ignore': a low-coverage region/contrast can exit non-zero; failure stays in the run report (see docs/pipeline_notes.md#modkit_dmr) + // see docs/pipeline_notes.md#modkit_dmr errorStrategy 'ignore' input: @@ -42,7 +42,7 @@ process MODKIT_DMR { arg="\${spec#*:}" outdir="dmr_regions/\${bname}_\${tag}" mkdir -p "\$outdir" - # \$arg MUST stay unquoted so it splits into separate tokens (see docs/pipeline_notes.md#modkit_dmr) + # see docs/pipeline_notes.md#modkit_dmr modkit dmr multi "\${SARGS[@]}" \\ --regions-bed "\$bed" \\ --ref ${reference} \\ @@ -57,7 +57,7 @@ process MODKIT_DMR { done mkdir -p dmr_sites - # CRITICAL: array MUST be GROUP_LIST not GROUPS (GROUPS is a bash special var, assignments ignored) (see docs/pipeline_notes.md#modkit_dmr) + # see docs/pipeline_notes.md#modkit_dmr GROUP_LIST=(\$(cut -f2 ${manifest} | sort -u | grep . || true)) for ((x=0; x<\${#GROUP_LIST[@]}; x++)); do for ((y=x+1; y<\${#GROUP_LIST[@]}; y++)); do diff --git a/src/modules/modkit_entropy.nf b/src/modules/modkit_entropy.nf index 6d3b8bfa976dfba93f4a2641148b6ecca389f4de..20fd99fa40bb4e4e97e5bc0d5ba1a452548737d6 100644 --- a/src/modules/modkit_entropy.nf +++ b/src/modules/modkit_entropy.nf @@ -1,8 +1,8 @@ -// MODKIT_ENTROPY (see docs/pipeline_notes.md#modkit_entropy) +// see docs/pipeline_notes.md#modkit_entropy process MODKIT_ENTROPY { publishDir "${params.modkit_out_dir}/entropy/", mode: "copy", overwrite: true label 'cpu' - // errorStrategy 'ignore': a low-coverage group/region exits non-zero; let rest of step 3 finish (see docs/pipeline_notes.md#modkit_entropy) + // see docs/pipeline_notes.md#modkit_entropy errorStrategy 'ignore' input: @@ -23,7 +23,7 @@ process MODKIT_ENTROPY { samtools faidx ${reference} mkdir -p entropy - # Array MUST be GROUP_LIST not GROUPS: GROUPS is special, assignments are silently ignored (see docs/pipeline_notes.md#modkit_entropy) + # see docs/pipeline_notes.md#modkit_entropy GROUP_LIST=(\$(cut -f2 ${manifest} | sort -u | grep . || true)) for grp in "\${GROUP_LIST[@]}"; do SARGS=() @@ -34,19 +34,17 @@ process MODKIT_ENTROPY { for bed in ${region_beds}; do bname=\$(basename "\$bed" .bed) for base in A C; do - # NO --combine-strands: heavy/light strands must stay separate (see docs/pipeline_notes.md#modkit_entropy) - # --regions requires --out-bed to be a DIRECTORY, one per group/region/base (see docs/pipeline_notes.md#modkit_entropy) + # see docs/pipeline_notes.md#modkit_entropy outdir="entropy/\${grp}_\${bname}_\${base}" mkdir -p "\$outdir" - # --filter-threshold 0.75 MUST match the pileup step (see docs/pipeline_notes.md#modkit_entropy) modkit entropy "\${SARGS[@]}" \\ --ref ${reference} \\ --base "\$base" \\ --regions "\$bed" \\ --out-bed "\$outdir" \\ --prefix "\${grp}_\${bname}_\${base}" \\ - --filter-threshold 0.75 \\ - --min-coverage 3 \\ + --filter-percentile ${params.modkit_filter_percentile} \\ + --min-coverage ${params.modkit_min_coverage} \\ --threads ${threads} \\ --force --header --drop-zeros \\ --log-filepath "entropy_\${grp}_\${bname}_\${base}.log" \\ diff --git a/src/modules/modkit_localize.nf b/src/modules/modkit_localize.nf index 6685ee68172dd81d71a57159d8d09b6371c02524..aed77d1a81e67ba21b17e590689893843e2a3363 100644 --- a/src/modules/modkit_localize.nf +++ b/src/modules/modkit_localize.nf @@ -1,8 +1,8 @@ process MODKIT_LOCALIZE { - // .tbi taken as explicit input so Nextflow stages it next to .bed.gz in same workdir (see docs/pipeline_notes.md#modkit_localize) + // see docs/pipeline_notes.md#modkit_localize publishDir "${params.modkit_out_dir}/localize/", mode: "copy", overwrite: true label 'cpu' - // errorStrategy 'ignore': low-coverage sample may exit non-zero; don't sink whole run (see docs/pipeline_notes.md#modkit_localize) + // see docs/pipeline_notes.md#modkit_localize errorStrategy 'ignore' input: @@ -24,7 +24,7 @@ process MODKIT_LOCALIZE { --regions ${control_elements_bed} \\ --genome-sizes genome_sizes.tsv \\ --window 200 \\ - --min-coverage 3 \\ + --min-coverage ${params.modkit_min_coverage} \\ --threads ${threads} \\ --out-file "${id}_localize.tsv" \\ --log-filepath "${id}_localize.log" diff --git a/src/modules/modkit_motif.nf b/src/modules/modkit_motif.nf index 871f089bc797ec9a2688120b21f694fa543194ab..1372e0e9cc6529023446c1d24450c612b40a5da6 100644 --- a/src/modules/modkit_motif.nf +++ b/src/modules/modkit_motif.nf @@ -1,8 +1,7 @@ process MODKIT_MOTIF_SEARCH { - // .tbi must be an explicit input so Nextflow stages it next to .bed.gz (see docs/pipeline_notes.md#modkit_motif) + // see docs/pipeline_notes.md#modkit_motif publishDir "${params.modkit_out_dir}/motif/", mode: "copy", overwrite: true label 'cpu' - // errorStrategy 'ignore': exploratory per-sample step, low coverage yields no motif (see docs/pipeline_notes.md#modkit_motif) errorStrategy 'ignore' input: @@ -18,12 +17,12 @@ process MODKIT_MOTIF_SEARCH { """ set -euo pipefail samtools faidx ${reference} - # --min-sites lowered 300->10 for 16.5kb mtDNA; default finds nothing (see docs/pipeline_notes.md#modkit_motif) + # see docs/pipeline_notes.md#modkit_motif modkit motif search \\ --in-bedmethyl ${bed_gz} \\ --ref ${reference} \\ --contig ${params.mt_contig} \\ - --min-coverage 5 \\ + --min-coverage ${params.modkit_min_coverage} \\ --min-sites 10 \\ --known-motif GATC 1 a \\ --threads ${threads} \\ diff --git a/src/modules/mosdepth.nf b/src/modules/mosdepth.nf index aa0abc1fe63ca9b5044533d35448cacb48d3ba17..567be05f435bda357c0035231b8647ee52ca6d4e 100644 --- a/src/modules/mosdepth.nf +++ b/src/modules/mosdepth.nf @@ -15,7 +15,7 @@ process MOSDEPTH { """ set -euo pipefail mosdepth -t ${threads} -n -x "${id}" "${filtered_bam}" - # mosdepth default -F 1796 keeps supplementary, drops secondary; must match modkit --allow-non-primary (see docs/pipeline_notes.md#mosdepth) + # see docs/pipeline_notes.md#mosdepth mosdepth -t ${threads} --chrom ${params.mt_contig} "${id}_mt" "${filtered_bam}" """ } diff --git a/src/modules/mt_variant_scan.nf b/src/modules/mt_variant_scan.nf new file mode 100644 index 0000000000000000000000000000000000000000..38079abafcea85282e27e5caf2ed4017ae8f38cc --- /dev/null +++ b/src/modules/mt_variant_scan.nf @@ -0,0 +1,101 @@ +// see docs/pipeline_notes.md#mt_variant_scan + +process MT_ALLELE_COUNTS { + publishDir "${params.modkit_out_dir}/mt_variants/counts/", mode: "copy", overwrite: true + label 'cpu' + + input: + tuple val(id), path(bam), path(bai) + path mt_fasta + + output: + path "${id}.counts.tsv", emit: counts + + script: + """ + set -euo pipefail + mkdir -p in + ln -sf "\$(readlink -f "${bam}")" "in/${id}-mt_realigned.bam" + ln -sf "\$(readlink -f "${bai}")" "in/${id}-mt_realigned.bam.bai" + cp "${mt_fasta}" mt_contig.fa + BAM_DIRS=in REF_MT=mt_contig.fa MT_CONTIG="${params.mt_contig}" OUT_COUNTS=. \\ + bash "\$(command -v mt_allele_counts.sh)" + """ +} + +process MT_VARIANT_CALL { + container "${params.mt_variant_r_sif}" + publishDir "${params.modkit_out_dir}/mt_variants/", mode: "copy", overwrite: true + label 'cpu' + + input: + path counts, stageAs: "counts/*" + path mt_fasta + path samplesheet + + output: + path "sites_called.tsv", emit: sites_called + path "sites/*.sites.bed", emit: sites, optional: true + path "*.tsv", emit: tables + path "*.png", emit: figures, optional: true + + script: + def r_libs = params.mt_variant_r_libs ? "export R_LIBS_USER=\"${params.mt_variant_r_libs}\"" : "true" + """ + set -euo pipefail + ${r_libs} + Rscript "\$(command -v mt_spectrum_phasing.R)" --stage call --counts counts --ref "${mt_fasta}" \\ + --contig "${params.mt_contig}" --samplesheet "${samplesheet}" \\ + --id-col real_barcode --group-col group --run-col run_id \\ + --out . --cohort "${params.run_id ?: params.project_name}" \\ + --min-depth ${params.mt_variant_min_depth} --min-vaf ${params.mt_variant_min_vaf} \\ + --max-vaf ${params.mt_variant_max_vaf} --min-alt-strand ${params.mt_variant_min_alt_strand} + """ +} + +process MT_READ_BASES { + publishDir "${params.modkit_out_dir}/mt_variants/readbases/", mode: "copy", overwrite: true + label 'cpu' + + input: + tuple val(id), path(bam), path(bai), path(sites_bed) + path mt_fasta + + output: + path "${id}.readbases.tsv", emit: bases + + script: + """ + set -euo pipefail + mkdir -p in sites + ln -sf "\$(readlink -f "${bam}")" "in/${id}-mt_realigned.bam" + ln -sf "\$(readlink -f "${bai}")" "in/${id}-mt_realigned.bam.bai" + ln -sf "\$(readlink -f "${sites_bed}")" "sites/${id}.sites.bed" + cp "${mt_fasta}" mt_contig.fa + BAM_DIRS=in REF_MT=mt_contig.fa MT_CONTIG="${params.mt_contig}" SITES_DIR=sites OUT_BASES=. \\ + bash "\$(command -v mt_read_bases.sh)" + """ +} + +process MT_PHASING { + container "${params.mt_variant_r_sif}" + publishDir "${params.modkit_out_dir}/mt_variants/", mode: "copy", overwrite: true + label 'cpu' + + input: + path bases, stageAs: "readbases/*" + path sites_called + + output: + path "phasing_*", emit: phasing, optional: true + + script: + def r_libs = params.mt_variant_r_libs ? "export R_LIBS_USER=\"${params.mt_variant_r_libs}\"" : "true" + """ + set -euo pipefail + ${r_libs} + mkdir -p readbases + Rscript "\$(command -v mt_spectrum_phasing.R)" --stage phase --out . --bases readbases \\ + --min-cover ${params.mt_variant_min_cover} + """ +} diff --git a/src/modules/mtdna_metrics.nf b/src/modules/mtdna_metrics.nf index 038f68dec497fa8cef5769e454e60e3e646f76d9..c249c13d284d646be3ae2ba768f021cfb175c2f6 100644 --- a/src/modules/mtdna_metrics.nf +++ b/src/modules/mtdna_metrics.nf @@ -25,7 +25,7 @@ process MTDNA_METRICS { meanmapq=\${meanmapq:-0} numreads_cov=\${numreads_cov:-0} - # -Q (mapQ) omitted on purpose: step 2 already filtered MAPQ>=10 (see docs/pipeline_notes.md#mtdna_metrics) + # see docs/pipeline_notes.md#mtdna_metrics meandepth_q10=\$(samtools depth -a -q 10 -r "${params.mt_contig}" "${filtered_bam}" | awk '{s+=\$3; n++} END {if (n>0) printf "%.4f", s/n; else print 0}') meandepth_q10=\${meandepth_q10:-0} @@ -33,13 +33,13 @@ process MTDNA_METRICS { reads_MT=\$(echo "\${idx}" | awk -v mt="${params.mt_contig}" '\$1==mt {print \$3}') reads_MT=\${reads_MT:-0} reads_nuclear=\$(echo "\${idx}" | awk -v mt="${params.mt_contig}" '\$1!=mt && \$1!="*" {s+=\$3} END {print s+0}') - # ratio is capture-enrichment efficiency, NOT mtDNA copy number (see docs/pipeline_notes.md#mtdna_metrics) + # see docs/pipeline_notes.md#mtdna_metrics ratio_MT_nuclear=\$(awk -v a="\${reads_MT}" -v b="\${reads_nuclear}" 'BEGIN {d=b; if (d<1) d=1; printf "%.6f", a/d}') readlen_mean=\$(samtools stats "${filtered_bam}" | grep '^SN' | grep 'average length:' | awk -F'\\t' '{print \$3}') readlen_mean=\${readlen_mean:-0} - # 0x800=supplementary, 0x900=supp+secondary; high supp/primary is expected on circular MT, not an error (see docs/pipeline_notes.md#mtdna_metrics) + # see docs/pipeline_notes.md#mtdna_metrics supp_MT=\$(samtools view -c -f 0x800 "${total_bam}" "${params.mt_contig}") primary_MT=\$(samtools view -c -F 0x900 "${total_bam}" "${params.mt_contig}") supp_MT=\${supp_MT:-0} @@ -92,7 +92,7 @@ process MERGE_MTDNA_METRICS { > mtdna_metrics_master.tsv cat mtdna_metrics_body.tsv >> mtdna_metrics_master.tsv - # custom-content '# ...' lines MUST precede the column header or MultiQC ignores the table (see docs/pipeline_notes.md#mtdna_metrics) + # see docs/pipeline_notes.md#mtdna_metrics echo "# plot_type: 'table'" > "mtDNA_Metrics_mqc.tsv" echo "# id: 'mtdna metrics custom'" >> "mtDNA_Metrics_mqc.tsv" echo "# section_name: 'mtDNA per-sample metrics'" >> "mtDNA_Metrics_mqc.tsv" diff --git a/src/modules/nanocomp.nf b/src/modules/nanocomp.nf index 81d3bd859619d4b1dcea5bbc080e4839745be2d7..4a92d72f68f43d76f6ae10a171ce4acd647e0e3d 100644 --- a/src/modules/nanocomp.nf +++ b/src/modules/nanocomp.nf @@ -15,10 +15,10 @@ process NANOCOMP { """ set -euo pipefail mkdir -p nanocomp - # build --bam and --names in ONE pass; NanoComp matches Nth bam to Nth name (see docs/pipeline_notes.md#nanocomp) + # see docs/pipeline_notes.md#nanocomp BAMS=(); NAMES=() while IFS=\$'\\t' read -r f n; do - # skip blank/trailing-newline lines or an empty bam path gets pushed (see docs/pipeline_notes.md#nanocomp) + # see docs/pipeline_notes.md#nanocomp [ -z "\$f" ] && continue BAMS+=( "\$f" ); NAMES+=( "\$n" ) done < ${manifest} diff --git a/src/modules/pycoqc.nf b/src/modules/pycoqc.nf index 0865ec1a75fddd8838073e096838539f624c5a3e..bc74dfa03a9ef3e9cdbf7008e9fcacf2c2842c4d 100644 --- a/src/modules/pycoqc.nf +++ b/src/modules/pycoqc.nf @@ -1,7 +1,7 @@ process PYCOQC_NO_FILTER { publishDir "${params.qc_out_dir}/pycoqc_no_filter/", mode: 'copy', overwrite: true, pattern: "*-Unfiltered_pycoqc*" label 'cpu' - // errorStrategy 'ignore': empty barcode -> empty seq_summary -> pycoQC EmptyDataError; drop only that sample's plot, not the step (see docs/pipeline_notes.md#pycoqc) + // see docs/pipeline_notes.md#pycoqc errorStrategy 'ignore' input: @@ -40,7 +40,8 @@ process PYCOQC_NO_FILTER { process PYCOQC_FILTER { publishDir "${params.qc_out_dir}/pycoqc_filtered/", mode: 'copy', overwrite: true, pattern: "*-Filtered_pycoqc*" label 'cpu' - errorStrategy 'ignore' // like PYCOQC_NO_FILTER: bad/empty sample skips its plot, not the step (see docs/pipeline_notes.md#pycoqc) + // see docs/pipeline_notes.md#pycoqc + errorStrategy 'ignore' input: val id diff --git a/src/modules/relabel_samples.nf b/src/modules/relabel_samples.nf index 94d40c8ce2619c80c405271b444aeab8d04336db..57cfd5ebdb3a89d5a39531ea76693b030acecf71 100644 --- a/src/modules/relabel_samples.nf +++ b/src/modules/relabel_samples.nf @@ -16,7 +16,8 @@ process RELABEL_SAMPLES { script: """ set -euo pipefail - ulimit -c 0 || true # dorado 2.0.0 summary SIGABRTs -> avoid multi-GB core dumps (see docs/pipeline_notes.md#relabel_samples) + # see docs/pipeline_notes.md#relabel_samples + ulimit -c 0 || true if [ ! -s "${samplesheet}" ]; then echo "[RELABEL] ERROR: samplesheet '${samplesheet}' is empty/missing." >&2 @@ -30,7 +31,7 @@ process RELABEL_SAMPLES { inbams="" for bc in \$bcs; do - # match by barcode token: dorado emits ..._barcode07.bam (no trailing token); _barcode07_*.bam never matches (see docs/pipeline_notes.md#relabel_samples) + # see docs/pipeline_notes.md#relabel_samples for f in *\${bc}*.bam; do [ -e "\$f" ] && inbams="\$inbams \$f" done @@ -43,16 +44,16 @@ process RELABEL_SAMPLES { n=\$(echo \$inbams | wc -w) if [ "\$n" -eq 1 ]; then - cp \$inbams "\${real}.bam" # demux BAM already coordinate-sorted from BASECALL; no re-sort (see docs/pipeline_notes.md#relabel_samples) + cp \$inbams "\${real}.bam" else - # merging coordinate-sorted inputs stays coordinate-sorted; no re-sort (see docs/pipeline_notes.md#relabel_samples) + # see docs/pipeline_notes.md#relabel_samples echo "[RELABEL] merging \$n barcodes into \$real (\$inbams)" samtools merge -f -@ ${samtools_threads} "\${real}.bam" \$inbams fi samtools index -@ ${samtools_threads} "\${real}.bam" - # dorado 2.0.0 summary SIGABRTs on the MERGED bam (@PG collisions); summarize each INPUT bam and concat (see docs/pipeline_notes.md#relabel_samples) + # see docs/pipeline_notes.md#relabel_samples : > "\${real}.txt" hdr=1 for f in \$inbams; do diff --git a/src/modules/rescue_mt.nf b/src/modules/rescue_mt.nf new file mode 100644 index 0000000000000000000000000000000000000000..55319408e2f975eb20f7166d387477e407ad128b --- /dev/null +++ b/src/modules/rescue_mt.nf @@ -0,0 +1,61 @@ +process RESCUE_MT { + publishDir "${params.qc_out_dir}/bam_filtering_mt/", mode: "copy", pattern: "*.mt.ba*", overwrite: true + publishDir "${params.qc_out_dir}/mt_rescue/", mode: "copy", pattern: "${id}.*.tsv", overwrite: true + label 'cpu' + + input: + tuple val(id), path(bam), path(bai) + path mt_fasta + path numt_bed + val samtools_threads + + output: + tuple val(id), path("${id}.mt.bam"), path("${id}.mt.bam.bai"), emit: mt_bam + path "${id}.rescue_summary.tsv", emit: summary + path "${id}.depth250.tsv", emit: depth + path "${id}.clip_dropped.tsv", emit: dropped + + script: + """ + set -euo pipefail + mkdir -p in + ln -sf "\$(readlink -f "${bam}")" "in/${id}-Unfiltered.bam" + ln -sf "\$(readlink -f "${bai}")" "in/${id}-Unfiltered.bam.bai" + cp "${mt_fasta}" mt_contig.fa + SRC_DIR=in OUT_DIR=rescue REF_MT=mt_contig.fa NUMT_BED="${numt_bed}" MT_CONTIG="${params.mt_contig}" \\ + SAMPLES="${id}" THREADS=${samtools_threads} MAX_CLIP=${params.rescue_max_clip} STAGE=realign FORCE=1 \\ + LABEL="${params.run_id ?: params.project_name}" CODE_REV="${workflow.commitId ?: 'uncommitted'}" \\ + bash "\$(command -v rescue_mt.sh)" + mv "rescue/qc/bam_filtering/${id}-Filtered_${params.mt_contig}_realigned.bam" "${id}.mt.bam" + mv "rescue/qc/bam_filtering/${id}-Filtered_${params.mt_contig}_realigned.bam.bai" "${id}.mt.bam.bai" + cp rescue/qc/realign_summary.tsv "${id}.rescue_summary.tsv" + cp "rescue/qc/mt_depth/${id}.depth250.tsv" "${id}.depth250.tsv" + cp "rescue/qc/candidates/${id}.clip_dropped.tsv" "${id}.clip_dropped.tsv" + """ +} + +process RESCUE_MT_SUMMARY { + publishDir "${params.qc_out_dir}/mt_rescue/", mode: "copy", overwrite: true + label 'cpu' + + input: + path summaries + + output: + path "rescue_summary.tsv", emit: table + + script: + def force = params.rescue_force ? "1" : "0" + """ + set -euo pipefail + awk 'FNR==1 && NR!=1 {next} {print}' ${summaries} > rescue_summary.tsv + awk -F'\\t' -v limit=${params.rescue_max_numt_frac} -v force=${force} ' + NR==1 { for (i=1;i<=NF;i++) col[\$i]=i; next } + { p=\$col["realigned_primary"]+0; d=\$col["clip_dropped_numt"]+0; f=(p>0 ? d/p : 0) + printf "[RESCUE_MT] %s: primaries %d, NUMT-origin %d, dropped as nuclear %d (%.4f of primaries), windows below 3 reads %d\\n", \$1, p, \$col["from_numt"], d, f, \$col["depth_windows_lt3"] + if (f > limit) over++ } + END { if (over > 0 && force != "1") { + printf "[RESCUE_MT] %d sample(s) exceed rescue_max_numt_frac=%s: nuclear reads reach the mt contig through the NUMT loci; judge rescue_summary.tsv, see docs/pipeline_notes.md#numt_rescue, then rerun with --rescue_force true or disable --rescue_mt\\n", over, limit > "/dev/stderr" + exit 1 } }' rescue_summary.tsv + """ +} diff --git a/src/nextflow.config b/src/nextflow.config index e4b0faf01d7ce1d6658735fe6b585ec01f13bfa5..ad788c3b0e4787f30abda2bd767afa8ea825d0a5 100644 --- a/src/nextflow.config +++ b/src/nextflow.config @@ -1,5 +1,6 @@ cleanup = true params { + // parameters are documented in README.md (Pipeline parameters) and docs/pipeline_notes.md#nextflow project_name = "default" run_id = null step = null @@ -28,10 +29,10 @@ params { steps_2_input_directory = params.run_id ? params.basecalling_out_dir : null steps_3_input_directory = params.run_id ? params.qc_out_dir : null - // mt_contig must match the reference AND the dorado-aligned BAM (e.g. MT vs chrM); see docs/pipeline_notes.md#nextflow + // see docs/pipeline_notes.md#nextflow mt_contig = "MT" - // samplesheet holds participant clinical data; keep OUTSIDE version control (gitignored); see docs/pipeline_notes.md#nextflow + // see docs/pipeline_notes.md#nextflow samplesheet = "/home/joaochrusciel/bio_projects/nanopore/references/samplesheet.csv" mt_macro_bed = "/home/joaochrusciel/bio_projects/nanopore/references/mt_macro.bed" @@ -44,12 +45,24 @@ params { enable_haplogroup = true clair3_sif = "./images/clair3.sif" clair3_model = "r1041_e82_400bps_hac_v500" - haplogrep_tree = "phylotree-rcrs@17.2" // MUST match HAPLOGREP_TREE baked in containers/versions.txt; see docs/pipeline_notes.md#nextflow - - mapq = "10" + haplogrep_tree = "phylotree-rcrs@17.2" + + // see docs/pipeline_notes.md#heteroplasmy + enable_heteroplasmy = false + enable_haplocheck = true + mtdnaserver_sif = "./images/mtdna-server-2.sif" + het_detection_level = "0.01" + het_downsample_target = 40 + het_downsample_seed = 42 + mt_keep_supplementary = true + het_mapq = 20 + het_baseq = 20 + het_strand_bias = "1.6" + + mapq = "10" qscore_thresh = "9" basecall_speed = "hac@v6.0.0" - basecall_mods = "5mC_5hmC,6mA" // can't use more than one modification per nucleotide; see docs/pipeline_notes.md#nextflow + basecall_mods = "5mC_5hmC,6mA" barcoding_kit = "SQK-RBK114-24" min_mapped_reads_thresh = 500 basecall_config = null @@ -65,6 +78,29 @@ params { modkit_threads = 8 modkit_extract_qsize = 1000 modkit_isize = 10000 + // see docs/pipeline_notes.md#modkit + modkit_filter_percentile = "0.1" + modkit_filter_mode = "percentile" + modkit_filter_threshold = "0.75" + modkit_num_reads = 50000 + modkit_min_coverage = 10 + + // see docs/pipeline_notes.md#numt_rescue + rescue_mt = false + numt_bed = null + rescue_max_clip = 300 + rescue_max_numt_frac = 0.01 + rescue_force = false + + // see docs/pipeline_notes.md#mt_variant_scan + enable_mt_variant_scan = false + mt_variant_r_sif = "./images/rstudio_tidyverse_4.5.3.sif" + mt_variant_r_libs = null + mt_variant_min_depth = 20 + mt_variant_min_vaf = "0.02" + mt_variant_max_vaf = "0.9" + mt_variant_min_alt_strand = 3 + mt_variant_min_cover = 20 } includeConfig ({ diff --git a/src/sub_workflows/FILTERING_AND_QC_FROM_MINKNOW.nf b/src/sub_workflows/FILTERING_AND_QC_FROM_MINKNOW.nf index dca6c2fab723cc3e9edb71fad74d283110e6bc67..4932ebde316df319bdd248f993706ad983fd2316 100644 --- a/src/sub_workflows/FILTERING_AND_QC_FROM_MINKNOW.nf +++ b/src/sub_workflows/FILTERING_AND_QC_FROM_MINKNOW.nf @@ -2,6 +2,7 @@ include { PYCOQC_NO_FILTER ; PYCOQC_FILTER } from '../modules/pycoqc.nf' include { FILTER_BAM } from '../modules/filter_bam.nf' include { MAKE_QC_REPORT } from '../modules/num_reads_report.nf' include { CONVERT_INPUT_FROM_MINKNOW_BARCODED ; CONVERT_INPUT_FROM_MINKNOW_NOT_BARCODED } from '../modules/convert_input_from_minknow.nf' +include { RESCUE_MT ; RESCUE_MT_SUMMARY } from '../modules/rescue_mt.nf' workflow FILTERING_AND_QC_FROM_MINKNOW { take: @@ -9,6 +10,8 @@ workflow FILTERING_AND_QC_FROM_MINKNOW { mapq qscore_thresh samtools_threads + mt_reference + numt_bed main: if (params.is_barcoded) { @@ -24,6 +27,13 @@ workflow FILTERING_AND_QC_FROM_MINKNOW { mapq, samtools_threads ) + // see docs/pipeline_notes.md#numt_rescue + if (params.rescue_mt) { + tbam_keyed = FILTER_BAM.out.total_bam.map { f -> tuple(f.baseName.tokenize('-')[0], f) } + tbai_keyed = FILTER_BAM.out.total_bai.map { f -> tuple(f.name.tokenize('-')[0], f) } + RESCUE_MT(tbam_keyed.join(tbai_keyed), mt_reference, numt_bed, samtools_threads) + RESCUE_MT_SUMMARY(RESCUE_MT.out.summary.collect()) + } PYCOQC_NO_FILTER( FILTER_BAM.out.id, FILTER_BAM.out.total_bam, diff --git a/src/sub_workflows/FILTERING_AND_QC_FROM_STEP_1.nf b/src/sub_workflows/FILTERING_AND_QC_FROM_STEP_1.nf index 278bb321e1c3b6724f0c82ec906f772c7ff291ae..6477e3e67c84a42984704b925dd3d786016bb52e 100644 --- a/src/sub_workflows/FILTERING_AND_QC_FROM_STEP_1.nf +++ b/src/sub_workflows/FILTERING_AND_QC_FROM_STEP_1.nf @@ -6,6 +6,7 @@ include { SAMTOOLS_STATS } from '../modules/samtools_qc.nf' include { MOSDEPTH } from '../modules/mosdepth.nf' include { NANOPLOT } from '../modules/nanoplot.nf' include { MTDNA_METRICS } from '../modules/mtdna_metrics.nf' +include { RESCUE_MT ; RESCUE_MT_SUMMARY } from '../modules/rescue_mt.nf' workflow FILTERING_AND_QC_FROM_STEP_1 { take: @@ -14,6 +15,8 @@ workflow FILTERING_AND_QC_FROM_STEP_1 { mapq qscore_thresh samtools_threads + mt_reference + numt_bed main: FILTER_BAM(bam_files, txt_files, mapq, samtools_threads) @@ -23,7 +26,7 @@ workflow FILTERING_AND_QC_FROM_STEP_1 { FILTER_BAM.out.filtered_bai, samtools_threads ) - // join BAM/BAI by real_barcode token (baseName before '-'), NOT emission order, or dropped samples misalign id<->BAM (see docs/pipeline_notes.md#filtering_and_qc_from_step_1) + // see docs/pipeline_notes.md#filtering_and_qc_from_step_1 fbam_keyed = FILTER_BAM.out.filtered_bam.map { f -> tuple(f.baseName.tokenize('-')[0], f) } fbai_keyed = FILTER_BAM.out.filtered_bai.map { f -> tuple(f.name.tokenize('-')[0], f) } tbam_keyed = FILTER_BAM.out.total_bam.map { f -> tuple(f.baseName.tokenize('-')[0], f) } @@ -38,6 +41,11 @@ workflow FILTERING_AND_QC_FROM_STEP_1 { if (params.enable_nanopack) { NANOPLOT(qc_in, samtools_threads) } + // see docs/pipeline_notes.md#numt_rescue + if (params.rescue_mt) { + RESCUE_MT(tbam_keyed.join(tbai_keyed), mt_reference, numt_bed, samtools_threads) + RESCUE_MT_SUMMARY(RESCUE_MT.out.summary.collect()) + } PYCOQC_NO_FILTER( FILTER_BAM.out.id, FILTER_BAM.out.total_bam, diff --git a/src/sub_workflows/MODKIT_AND_MULTIQC.nf b/src/sub_workflows/MODKIT_AND_MULTIQC.nf index 7a9072ecda6b078b0382584de1a836d48f9c6f68..82d4f87b25fdb0da38bf4073f91114323c47df42 100644 --- a/src/sub_workflows/MODKIT_AND_MULTIQC.nf +++ b/src/sub_workflows/MODKIT_AND_MULTIQC.nf @@ -9,6 +9,9 @@ include { MODKIT_LOCALIZE } from '../modules/modkit_localize.nf' include { MODKIT_MOTIF_SEARCH } from '../modules/modkit_motif.nf' include { NANOCOMP } from '../modules/nanocomp.nf' include { CLAIR3 ; HAPLOGROUP ; MERGE_HAPLOGROUPS } from '../modules/haplogroup.nf' +include { EXTRACT_MT } from '../modules/extract_mt.nf' +include { MT_ALLELE_COUNTS ; MT_VARIANT_CALL ; MT_READ_BASES ; MT_PHASING } from '../modules/mt_variant_scan.nf' +include { DOWNSAMPLE_MATCH ; MUTSERVE as MUTSERVE_NATIVE ; MUTSERVE as MUTSERVE_MATCHED ; PARSE_HETEROPLASMY ; MERGE_HETEROPLASMY ; MERGE_VCFS_HET ; HAPLOCHECK_HET } from '../modules/heteroplasmy.nf' workflow MODKIT_AND_MULTIQC { take: @@ -26,6 +29,7 @@ workflow MODKIT_AND_MULTIQC { modkit_extract_qsize samplesheet mtdna_metrics_files + rescued_mt_bams main: CALCULATE_COVERAGE(filtered_bams, filtered_bais, samtools_threads) @@ -46,19 +50,26 @@ workflow MODKIT_AND_MULTIQC { MERGE_COVERAGE.out.mqc, multiqc_config ) - // do NOT re-add strand collapsing: mtDNA 6mA/5mC is strand-asymmetric (see docs/pipeline_notes.md#modkit_and_multiqc) - modkit_bams = filtered_bams - modkit_bais = filtered_bais + // see docs/pipeline_notes.md#modkit_and_multiqc + bam_bai = filtered_bams.map { id, bam -> tuple(id.tokenize('-')[0], bam) } + .join(filtered_bais.map { bai -> tuple(bai.name.tokenize('-')[0], bai) }) + // see docs/pipeline_notes.md#numt_rescue + if (params.rescue_mt) { + mt_bam = rescued_mt_bams + } else { + EXTRACT_MT(bam_bai) + mt_bam = EXTRACT_MT.out.mt_bam + } + MODKIT( - modkit_bams, - modkit_bais, + mt_bam, reference_file, modkit_threads, modkit_isize, modkit_extract_qsize ) - // keep .unique(): a real_barcode appears on >1 row for dual-barcode samples (see docs/pipeline_notes.md#modkit_and_multiqc) + // see docs/pipeline_notes.md#modkit_and_multiqc sample_groups = channel.fromPath(samplesheet) .splitCsv(header: true) .map { row -> tuple(row.real_barcode, row.group) } @@ -72,11 +83,11 @@ workflow MODKIT_AND_MULTIQC { ]) control_elements_ch = channel.value(file(params.mt_control_elements_bed)) - // join pileup to group by real_barcode (id token), then drop user_unknown (see docs/pipeline_notes.md#modkit_and_multiqc) + // see docs/pipeline_notes.md#modkit_and_multiqc pileup_labeled = MODKIT.out.pileup .map { id, bed, tbi -> tuple(id.tokenize('-')[0], bed, tbi) } .join(sample_groups) - dmr_eligible = pileup_labeled.filter { it[3] != 'user_unknown' } + dmr_eligible = pileup_labeled dmr_beds = dmr_eligible.map { rb, bed, tbi, grp -> bed }.collect() dmr_tbis = dmr_eligible.map { rb, bed, tbi, grp -> tbi }.collect() dmr_manifest = dmr_eligible @@ -84,9 +95,8 @@ workflow MODKIT_AND_MULTIQC { .collectFile(name: 'dmr_manifest.tsv', newLine: true) MODKIT_DMR(dmr_beds, dmr_tbis, dmr_manifest, reference_file, region_beds_ch, modkit_threads) - bam_bai = filtered_bams.map { id, bam -> tuple(id.tokenize('-')[0], bam) } - .join(filtered_bais.map { bai -> tuple(bai.name.tokenize('-')[0], bai) }) - ent_labeled = bam_bai.join(sample_groups) + // see docs/pipeline_notes.md#modkit_entropy + ent_labeled = mt_bam.join(sample_groups) ent_bams = ent_labeled.map { rb, bam, bai, grp -> bam }.collect() ent_bais = ent_labeled.map { rb, bam, bai, grp -> bai }.collect() ent_manifest = ent_labeled @@ -106,13 +116,39 @@ workflow MODKIT_AND_MULTIQC { NANOCOMP(nanocomp_bams, nanocomp_bais, nanocomp_manifest, modkit_threads) } - // keep .ifEmpty([]) so MERGE runs even when all Haplogrep3 are empty (see docs/pipeline_notes.md#modkit_and_multiqc) + // see docs/pipeline_notes.md#modkit_and_multiqc if (params.enable_haplogroup) { CLAIR3(bam_bai, reference_file, modkit_threads) HAPLOGROUP(CLAIR3.out.vcf, params.haplogrep_tree) MERGE_HAPLOGROUPS(HAPLOGROUP.out.haplogroup.collect().ifEmpty([]), samplesheet) } + // see docs/pipeline_notes.md#heteroplasmy + if (params.enable_heteroplasmy) { + MUTSERVE_NATIVE(mt_bam, reference_file, 'native') + DOWNSAMPLE_MATCH(mt_bam) + MUTSERVE_MATCHED(DOWNSAMPLE_MATCH.out.matched, reference_file, 'matched') + het_tables = MUTSERVE_NATIVE.out.table.mix(MUTSERVE_MATCHED.out.table) + PARSE_HETEROPLASMY(het_tables) + MERGE_HETEROPLASMY(PARSE_HETEROPLASMY.out.tsv.collect().ifEmpty([]), samplesheet) + if (params.enable_haplocheck) { + // see docs/pipeline_notes.md#heteroplasmy + nat_vcfs = MUTSERVE_NATIVE.out.vcf.map { sid, arm, vcf, tbi -> vcf }.collect() + nat_tbis = MUTSERVE_NATIVE.out.vcf.map { sid, arm, vcf, tbi -> tbi }.collect() + MERGE_VCFS_HET(nat_vcfs, nat_tbis) + HAPLOCHECK_HET(MERGE_VCFS_HET.out.merged) + } + } + + // see docs/pipeline_notes.md#mt_variant_scan + if (params.enable_mt_variant_scan) { + MT_ALLELE_COUNTS(mt_bam, reference_file) + MT_VARIANT_CALL(MT_ALLELE_COUNTS.out.counts.collect(), reference_file, samplesheet) + site_beds = MT_VARIANT_CALL.out.sites.flatten().map { bed -> tuple(bed.name.replaceFirst(/[.]sites[.]bed$/, ''), bed) } + MT_READ_BASES(mt_bam.join(site_beds), reference_file) + MT_PHASING(MT_READ_BASES.out.bases.collect().ifEmpty([]), MT_VARIANT_CALL.out.sites_called) + } + emit: modkit_pileup = MODKIT.out.pileup modkit_all = MODKIT.out.allfiles