From ee12f5de5f32567c103ecaab8bc79cd3da6f8200 Mon Sep 17 00:00:00 2001 From: Steven Roberts Date: Mon, 12 Oct 2026 16:42:35 -0700 Subject: [PATCH] Add a 3' Tag-seq mode to the RNA-seq workflow (--library_protocol three_prime_tag) 3' Tag-Seq libraries (Lexogen QuantSeq and similar) yield reads in proportion to a transcript's abundance, not its length, so length correction is wrong for them. The new --library_protocol parameter (default full_length, so existing runs are unchanged) selects: - Salmon --noLengthCorrection, the option Salmon documents for QuantSeq; - DESeq2 on raw counts (tximport countsFromAbundance "no" into DESeqDataSetFromMatrix), with no average-transcript-length offsets; with --replicate_unit, raw counts are summed per unit; - fastp --trim_poly_x for reads that run into the poly(A) tail. The protocol and the length model used are recorded in the manifest and the DESeq2 design JSON; an unknown protocol stops the run. Verified natively on simulated QuantSeq-like reads from the full M. galloprovincialis transcriptome: both modes recover 59/60 simulated 4-fold changes, but for 40 genes whose expressed isoform switched to a ~2.7x longer one with no change in 3' reads, full_length calls 26 differentially expressed (|log2FC| up to 3.0) and three_prime_tag calls none. Rerunning an earlier full_length smoke run with this code gives byte-identical DESeq2 results. CI stub-runs the three_prime_tag DAG and rejects an unknown protocol; tests/test_rnaseq_protocols.py guards the wiring. Needed by TUCCI2025_HEAT_PRIMING (PRJNA1132140, draft PR #41). Co-Authored-By: Claude Opus 5.5 --- .github/workflows/ci.yml | 21 +++++++ docs/candidate_studies.md | 2 +- modules/rnaseq/align_quant_salmon.nf | 7 +++ .../rnaseq/differential_expression_deseq2.nf | 27 ++++++++- modules/rnaseq/emit_manifest.nf | 1 + modules/rnaseq/trim_fastp.nf | 4 ++ tests/test_rnaseq_protocols.py | 59 +++++++++++++++++++ workflows/rnaseq/README.md | 30 ++++++++++ workflows/rnaseq/main.nf | 4 ++ workflows/rnaseq/nextflow.config | 6 ++ 10 files changed, 158 insertions(+), 3 deletions(-) create mode 100644 tests/test_rnaseq_protocols.py diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 747dd67..caf17e7 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -198,6 +198,27 @@ jobs: --viral_reference data/reference/pathogens/OsHV-1/NC_005881.2.fasta \ --outdir /tmp/aree-smoke/rnaseq-raw-se \ -work-dir /tmp/aree-smoke/work-rnaseq-raw-se + # 3' Tag-seq mode wires through the same DAG, and an unknown protocol stops the run. + nextflow run workflows/rnaseq -profile local,ci_stub -stub-run \ + --mode raw_reanalysis \ + --single_end \ + --library_protocol three_prime_tag \ + --study_id CI_RNASEQ_TAG \ + --comparison_id treatment_vs_control \ + --reads 'tests/fixtures/rnaseq_raw_se/ERR*.fastq' \ + --transcript_fasta tests/fixtures/rnaseq_raw/reference.fa \ + --gtf tests/fixtures/rnaseq_raw/annotation.gtf \ + --tx2gene tests/fixtures/rnaseq_raw/tx2gene.tsv \ + --sample_sheet tests/fixtures/rnaseq_raw_se/samples.tsv \ + --control_level control \ + --treatment_level treatment \ + --outdir /tmp/aree-smoke/rnaseq-raw-tag \ + -work-dir /tmp/aree-smoke/work-rnaseq-raw-tag + ! nextflow run workflows/rnaseq -profile local,ci_stub -stub-run \ + --mode raw_reanalysis --library_protocol quantseq \ + --study_id X --comparison_id y \ + --outdir /tmp/aree-smoke/rnaseq-bad-protocol \ + -work-dir /tmp/aree-smoke/work-rnaseq-bad-protocol # One trimmed file per library, none for a second mate. test "$(ls /tmp/aree-smoke/rnaseq-raw-se/rnaseq/trimmed/*_trimmed_R1.fastq.gz | wc -l)" -eq 4 ! ls /tmp/aree-smoke/rnaseq-raw-se/rnaseq/trimmed/*_trimmed_R2.fastq.gz 2>/dev/null diff --git a/docs/candidate_studies.md b/docs/candidate_studies.md index 613aff2..89e307f 100644 --- a/docs/candidate_studies.md +++ b/docs/candidate_studies.md @@ -94,7 +94,7 @@ pooling would run on anything but simulated data. |---|---|---|---|---|---| | 4 | `PRJNA678408` | 10 WGBS | 5 diploid + 5 triploid | Desiccation + acute heat | Would be AREE's **first real methylation study**. But the intake converter does not handle region tables yet, and coordinate-based evidence is exposed to the assembly change described in [handling_genome_versions.md](handling_genome_versions.md). | | 5 | `PRJNA762441` | 18 RNA-seq | diploid/triploid × 3 timepoints, n=3 | Thermal stress, ploidy contrast | Clean design; ploidy is a useful resilience covariate. No linked publication confirmed. | -| 6 | `PRJNA913164` | 72 | diploid/triploid, marine heatwave | Thermal tolerance | **Tag-seq (3′ counts), not standard RNA-seq** — quantification differs from the salmon-based workflow. Largest n on the list. No publication found; treat as unpublished data. | +| 6 | `PRJNA913164` | 72 | diploid/triploid, marine heatwave | Thermal tolerance | **Tag-seq (3′ counts), not standard RNA-seq.** Run it with `--library_protocol three_prime_tag` (no length correction in Salmon or DESeq2; see `workflows/rnaseq/README.md`). Largest n on the list. No publication found; treat as unpublished data. | | 7 ✅ | `PRJNA735889` | 76 RNA-seq (31 mislabelled AMPLICON) | 15 tanks, one per pH, × 5 oysters | OA "tipping point" | Registered as `LUTIER2022_OA_TIPPING`. **Do not filter on ENA library_strategy**: the AMPLICON runs are RNA-seq. One tank per pH, so fit on tanks, not oysters. | | 8 | `PRJNA1196326` | 6 RNA-seq | 2 groups × n=3 | Transgenerational OA | Very small, but transgenerational designs are directly relevant to breeding and rare in the public record. | | 9 | `PRJNA877226` | 6 RNA-seq | run labels blank in ENA | Vibrio × high temperature | Multi-stressor and cheap, but the deposited metadata does not describe groups — the design must be recovered from the paper before it is worth registering. | diff --git a/modules/rnaseq/align_quant_salmon.nf b/modules/rnaseq/align_quant_salmon.nf index 13d0be6..6ff3c10 100644 --- a/modules/rnaseq/align_quant_salmon.nf +++ b/modules/rnaseq/align_quant_salmon.nf @@ -64,13 +64,20 @@ process SALMON_QUANT { // Fragment GC-bias correction is applied to paired-end libraries only: // single-end reads do not observe fragment length, which the model needs, // and Salmon documents it as beta for single-end data. + // + // 3' Tag-Seq (params.library_protocol = 'three_prime_tag'): a transcript + // yields reads in proportion to its abundance, not its length, so the + // effective-length model would mis-weight them. Salmon documents + // --noLengthCorrection for exactly this case (QuantSeq). def fq = reads instanceof Collection ? reads as List : [ reads ] def input = fq.size() == 1 ? "-r ${fq[0]}" : "-1 ${fq[0]} -2 ${fq[1]} --gcBias" + def length_model = params.library_protocol == 'three_prime_tag' ? '--noLengthCorrection' : '' """ salmon quant \\ -i ${salmon_index} \\ -l A \\ ${input} \\ + ${length_model} \\ -p ${task.cpus} \\ --geneMap ${gtf} \\ --validateMappings \\ diff --git a/modules/rnaseq/differential_expression_deseq2.nf b/modules/rnaseq/differential_expression_deseq2.nf index 11c3ee7..5757dd5 100644 --- a/modules/rnaseq/differential_expression_deseq2.nf +++ b/modules/rnaseq/differential_expression_deseq2.nf @@ -61,6 +61,7 @@ process DIFFERENTIAL_EXPRESSION_DESEQ2 { exclude_samples = "${exclude_samples}", replicate_unit = "${replicate_unit}", covariates = "${covariates}", + library_protocol = "${params.library_protocol}", out_tsv = "${study_id}_${comparison_id}_deseq2_raw.tsv", out_rdata = "${study_id}_${comparison_id}_deseq2.RData", out_design = "${study_id}_${comparison_id}_deseq2_design.json" @@ -254,11 +255,29 @@ process DIFFERENTIAL_EXPRESSION_DESEQ2 { tx2gene\$transcript_id <- sub("\\\\..*", "", tx2gene\$transcript_id) tx2gene <- unique(tx2gene) + # 3' Tag-Seq counts do not scale with transcript length, so no length + # correction of any kind is applied: raw estimated counts go into DESeq2 + # without tximport's average-transcript-length offsets. Full-length RNA-seq + # keeps the offsets (or length-scaled counts when summing to units). + if (!args\$library_protocol %in% c("full_length", "three_prime_tag")) { + stop("library_protocol must be full_length or three_prime_tag, got: ", args\$library_protocol) + } + three_prime_tag <- identical(args\$library_protocol, "three_prime_tag") + length_model <- if (three_prime_tag) { + "none: raw counts, no length offsets (three_prime_tag)" + } else if (pseudo_bulk) { + "lengthScaledTPM counts summed per replicate unit" + } else { + "tximport average transcript length offsets" + } + message("Length model: ", length_model) + if (pseudo_bulk) { # Length-scaled counts can be summed across animals; tximport's - # per-sample length offsets cannot, so they are folded in here. + # per-sample length offsets cannot, so they are folded in here. For + # 3' Tag-Seq the raw counts are summed instead. txi <- tximport(quant_files, type = "salmon", tx2gene = tx2gene, ignoreTxVersion = TRUE, - countsFromAbundance = "lengthScaledTPM") + countsFromAbundance = if (three_prime_tag) "no" else "lengthScaledTPM") unit_counts <- t(rowsum(t(txi\$counts), group = samples\$replicate_unit)) unit_info <- unique(samples[, c("replicate_unit", "condition", covariates_used)]) unit_info\$n_libraries <- as.integer(table(samples\$replicate_unit)[unit_info\$replicate_unit]) @@ -291,6 +310,8 @@ process DIFFERENTIAL_EXPRESSION_DESEQ2 { dds <- if (pseudo_bulk) { DESeqDataSetFromMatrix(countData = round(unit_counts), colData = col_data, design = design_formula) + } else if (three_prime_tag) { + DESeqDataSetFromMatrix(countData = round(txi\$counts), colData = col_data, design = design_formula) } else { DESeqDataSetFromTximport(txi, colData = col_data, design = design_formula) } @@ -325,6 +346,8 @@ process DIFFERENTIAL_EXPRESSION_DESEQ2 { paste0(' "covariates_requested": ', json_arr(covariates_requested), ","), paste0(' "covariates_used": ', json_arr(covariates_used), ","), paste0(' "covariates_dropped_constant": ', json_arr(covariates_dropped), ","), + paste0(' "library_protocol": ', json_str(args\$library_protocol), ","), + paste0(' "length_model": ', json_str(length_model), ","), paste0(' "n_libraries_or_units": ', nrow(mm), ","), paste0(' "residual_df": ', residual_df), "}" diff --git a/modules/rnaseq/emit_manifest.nf b/modules/rnaseq/emit_manifest.nf index df4b367..4f2bc21 100644 --- a/modules/rnaseq/emit_manifest.nf +++ b/modules/rnaseq/emit_manifest.nf @@ -176,6 +176,7 @@ process EMIT_MANIFEST { "replicate_unit": replicate_unit, "covariates": covariates, "single_end": ${params.single_end ? 'True' : 'False'}, + "library_protocol": "${params.library_protocol}", "se_adapter": "${params.se_adapter}" if ${params.single_end ? 'True' : 'False'} else None, "viral_reference": "${params.viral_reference ?: ''}" or None, }, diff --git a/modules/rnaseq/trim_fastp.nf b/modules/rnaseq/trim_fastp.nf index 2bdeb79..713ac67 100644 --- a/modules/rnaseq/trim_fastp.nf +++ b/modules/rnaseq/trim_fastp.nf @@ -36,12 +36,16 @@ process TRIM_FASTP { // and is reliable, so it is kept. def fq = reads instanceof Collection ? reads as List : [ reads ] def se_adapter = params.se_adapter == 'auto' ? '' : "--adapter_sequence ${params.se_adapter}" + // 3' Tag-Seq reads run into the poly(A) tail; trim homopolymer tails so + // they neither fail to map nor map to A-rich sequence. + def poly_x = params.library_protocol == 'three_prime_tag' ? '--trim_poly_x' : '' def io = fq.size() == 1 \ ? "-i ${fq[0]} -o ${sample_id}_trimmed_R1.fastq.gz ${se_adapter}" \ : "-i ${fq[0]} -I ${fq[1]} -o ${sample_id}_trimmed_R1.fastq.gz -O ${sample_id}_trimmed_R2.fastq.gz --detect_adapter_for_pe" """ fastp \\ ${io} \\ + ${poly_x} \\ --thread ${task.cpus} \\ --json ${sample_id}.fastp.json \\ --html ${sample_id}.fastp.html \\ diff --git a/tests/test_rnaseq_protocols.py b/tests/test_rnaseq_protocols.py new file mode 100644 index 0000000..c65d169 --- /dev/null +++ b/tests/test_rnaseq_protocols.py @@ -0,0 +1,59 @@ +"""Guards on the RNA-seq workflow's library protocols (full_length, three_prime_tag). + +3' Tag-Seq (Lexogen QuantSeq and similar) needs no length correction anywhere: +not in Salmon, not in DESeq2. These tests read the workflow source, so they +run in the Python CI job without Nextflow, Salmon or R; they guard the wiring +against being undone. The behaviour itself was checked by running the workflow +on simulated Tag-seq reads in both modes (workflows/rnaseq/README.md), and CI +stub-runs the three_prime_tag DAG and rejects an unknown protocol. +""" +from __future__ import annotations + +import re + +from common import REPO_ROOT + +WORKFLOW = REPO_ROOT / "workflows" / "rnaseq" +MODULES = REPO_ROOT / "modules" / "rnaseq" + + +def _read(path): + return path.read_text(encoding="utf-8") + + +def test_full_length_is_the_default_so_existing_runs_are_unchanged(): + assert re.search(r"^\s*library_protocol\s*=\s*'full_length'", _read(WORKFLOW / "nextflow.config"), re.M) + + +def test_workflow_rejects_an_unknown_protocol(): + main = _read(WORKFLOW / "main.nf") + assert "params.library_protocol in ['full_length', 'three_prime_tag']" in main + assert "error \"--library_protocol must be" in main + + +def test_salmon_drops_length_correction_only_for_three_prime_tag(): + salmon = _read(MODULES / "align_quant_salmon.nf") + assert "params.library_protocol == 'three_prime_tag' ? '--noLengthCorrection' : ''" in salmon + assert "${length_model}" in salmon + + +def test_fastp_trims_poly_a_only_for_three_prime_tag(): + fastp = _read(MODULES / "trim_fastp.nf") + assert "params.library_protocol == 'three_prime_tag' ? '--trim_poly_x' : ''" in fastp + assert "${poly_x}" in fastp + + +def test_deseq2_uses_raw_counts_without_length_offsets_for_three_prime_tag(): + deseq2 = _read(MODULES / "differential_expression_deseq2.nf") + assert 'library_protocol = "${params.library_protocol}"' in deseq2 + # Per-library fit: raw counts into a plain count matrix, not FromTximport. + assert re.search( + r"else if \(three_prime_tag\) \{\s*DESeqDataSetFromMatrix\(countData = round\(txi\\\$counts\)", deseq2) + # Per-unit fit: raw counts summed, not length-scaled TPM. + assert 'countsFromAbundance = if (three_prime_tag) "no" else "lengthScaledTPM"' in deseq2 + # The length model actually used is written for the manifest. + assert '"length_model": ' in deseq2 and '"library_protocol": ' in deseq2 + + +def test_manifest_records_the_protocol(): + assert '"library_protocol": "${params.library_protocol}"' in _read(MODULES / "emit_manifest.nf") diff --git a/workflows/rnaseq/README.md b/workflows/rnaseq/README.md index f82647c..49d949a 100644 --- a/workflows/rnaseq/README.md +++ b/workflows/rnaseq/README.md @@ -77,6 +77,30 @@ sheet's `run_accession` column. The single-end path changes three steps: `single_end` and `se_adapter` are recorded in the manifest's `parameters`. +**3′ Tag-Seq libraries** (Lexogen QuantSeq and similar). Add +`--library_protocol three_prime_tag` (default `full_length`). Tag-seq reads +come from the 3′ end of each transcript, so a transcript yields reads in +proportion to its abundance, not its length. The mode changes three steps: + +- fastp also trims homopolymer tails (`--trim_poly_x`), since reads run into + the poly(A) tail. +- Salmon quantifies with `--noLengthCorrection`, the option Salmon documents + for QuantSeq-style protocols. +- DESeq2 gets the raw estimated counts (`tximport(countsFromAbundance = + "no")` into `DESeqDataSetFromMatrix`), with no average-transcript-length + offsets. With `--replicate_unit`, raw counts are summed per unit instead of + length-scaled ones. + +Why it matters, measured on simulated QuantSeq-like reads from the *M. +galloprovincialis* transcriptome (6 libraries, 3 vs 3, 2026-10-12). Both modes +recovered 59 of 60 genes with a simulated 4-fold change. But for 40 genes whose +expressed isoform switched to one about 2.7 times longer, with no change in +3′ reads, full-length mode called 26 differentially expressed (median |log2FC| +0.84, up to 3.0), because the per-sample length offsets moved with the isoform +mix. `three_prime_tag` called none (median |log2FC| 0.09). The protocol and +the length model are recorded in the manifest's `parameters` and +`deseq2_design`; any value other than the two above stops the run. + ### `processed_results_harmonization` An existing DE-like TSV (e.g. already-published DESeq2/edgeR/limma output) @@ -255,6 +279,12 @@ regenerated from raw data), not the evidence records themselves. tasks completed, through DESeq2, the manifest and the report. That checks the command lines, not the biology; no full-depth single-end run has been made. +- **3′ Tag-seq mode**: executed natively on 2026-10-12 against six simulated + QuantSeq-like single-end libraries (85,000 to 101,000 reads) and the full *M. + galloprovincialis* RefSeq transcriptome, in both protocols; see the + single-end section above for the result. Salmon received + `--noLengthCorrection` and fastp `--trim_poly_x` in every library (about + 6,500 poly(A) reads trimmed per library). No real Tag-seq study has been run. - **CI coverage**: tiny paired FASTQ, reference, annotation, tx2gene, and sample fixtures live under `../../tests/fixtures/rnaseq_raw`, and single-end ones (ENA-style file names, `run_accession` in the sample sheet) under diff --git a/workflows/rnaseq/main.nf b/workflows/rnaseq/main.nf index 7f6102d..e998e9b 100644 --- a/workflows/rnaseq/main.nf +++ b/workflows/rnaseq/main.nf @@ -67,6 +67,10 @@ workflow { no_deseq2_design = file("${projectDir}/assets/NO_DESEQ2_DESIGN") collect_versions_py = file("${projectDir}/../../modules/shared/collect_versions.py", checkIfExists: true) + if (!(params.library_protocol in ['full_length', 'three_prime_tag'])) { + error "--library_protocol must be 'full_length' or 'three_prime_tag', got '${params.library_protocol}'" + } + if (params.mode == 'raw_reanalysis') { // ---- raw_reanalysis branch -------------------------------------- diff --git a/workflows/rnaseq/nextflow.config b/workflows/rnaseq/nextflow.config index 6bcbc70..36cc949 100644 --- a/workflows/rnaseq/nextflow.config +++ b/workflows/rnaseq/nextflow.config @@ -33,6 +33,12 @@ params { reads = null // e.g. 'data/raw/GIGAS_HEAT01/*_R{1,2}.fastq.gz' single_end = false // true: --reads matches one FASTQ per library (sample id = file name stem) se_adapter = 'AGATCGGAAGAGCACACGTCTGAACTCCAGTCA' // single-end 3' adapter for fastp (Illumina TruSeq / NEBNext read 1); 'auto' = fastp detection + // 'full_length' (standard RNA-seq) or 'three_prime_tag' (3' Tag-Seq such as + // Lexogen QuantSeq): reads come from the 3' end of each transcript, so counts + // do not scale with transcript length. three_prime_tag turns off length + // correction in Salmon (--noLengthCorrection) and in DESeq2 (raw counts, + // no tximport length offsets), and trims poly(A) tails with fastp. + library_protocol = 'full_length' salmon_index = null // path to a prebuilt Salmon index directory transcript_fasta = null // used to build a Salmon index if salmon_index not given genome_fasta = null // reserved for STAR-based alternative alignment