Record the DESeq2 length model, guard --tag_seq, and document its measured effect

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 `--tag_seq` (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. |

modules/rnaseq/differential_expression_deseq2.nf

@@ -260,6 +260,15 @@ process DIFFERENTIAL_EXPRESSION_DESEQ2 {
     # a gene's count does not scale with transcript length. Their counts are
     # used as they are, with no tximport length offset or length scaling.
     if (args\$tag_seq) message("tag_seq: DESeq2 on counts without transcript-length offsets")
+    # Which length model was fitted, for the manifest: the one property that
+    # decides whether these fold changes can be compared with other studies'.
+    length_model <- if (args\$tag_seq) {
+        "none: counts without length offsets (tag_seq)"
+    } else if (pseudo_bulk) {
+        "lengthScaledTPM counts summed per replicate unit"
+    } else {
+        "tximport average transcript length offsets"
+    }
 
     if (pseudo_bulk) {
         # Length-scaled counts can be summed across animals; tximport's
@@ -336,6 +345,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('  "tag_seq": ', if (args\$tag_seq) "true" else "false", ","),
+        paste0('  "length_model": ', json_str(length_model), ","),
         paste0('  "n_libraries_or_units": ', nrow(mm), ","),
         paste0('  "residual_df": ', residual_df),
         "}"

tests/test_rnaseq_tag_seq.py

@@ -0,0 +1,55 @@
+"""Guards on the RNA-seq workflow's 3' Tag-seq path (--tag_seq).
+
+3' Tag-Seq (Lexogen QuantSeq and similar) needs no length correction anywhere:
+not in Salmon, not in DESeq2. Applying it can manufacture fold changes from
+isoform switching (workflows/rnaseq/README.md, "3' Tag-seq libraries"). These
+tests read the workflow source, so they run in the Python CI job without
+Nextflow, Salmon or R, and guard the wiring against being undone. CI's
+Nextflow job stub-runs the --tag_seq DAG.
+"""
+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_tag_seq_is_off_by_default():
+    assert re.search(r"^\s*tag_seq\s*=\s*false\b", _read(WORKFLOW / "nextflow.config"), re.M)
+
+
+def test_salmon_drops_length_correction_and_gc_bias_for_tag_seq():
+    salmon = _read(MODULES / "align_quant_salmon.nf")
+    assert "params.tag_seq ? '--noLengthCorrection' : ''" in salmon
+    assert "(params.tag_seq ? '' : ' --gcBias')" in salmon
+
+
+def test_fastp_trims_poly_a_for_tag_seq():
+    assert "params.tag_seq ? '--trim_poly_x' : ''" in _read(MODULES / "trim_fastp.nf")
+
+
+def test_deseq2_uses_counts_without_length_offsets_for_tag_seq():
+    deseq2 = _read(MODULES / "differential_expression_deseq2.nf")
+    # Per-library fit: counts into a plain matrix, not DESeqDataSetFromTximport.
+    assert re.search(r"else if \(args\\\$tag_seq\) \{.*?DESeqDataSetFromMatrix\(countData = round\(txi\\\$counts\)",
+                     deseq2, re.S)
+    # Per-unit fit: raw counts summed, not length-scaled TPM.
+    assert 'countsFromAbundance = if (args\\$tag_seq) "no" else "lengthScaledTPM"' in deseq2
+
+
+def test_the_length_model_used_is_written_for_the_manifest():
+    deseq2 = _read(MODULES / "differential_expression_deseq2.nf")
+    assert '"none: counts without length offsets (tag_seq)"' in deseq2
+    assert "paste0('  \"length_model\": ', json_str(length_model)" in deseq2
+    assert "paste0('  \"tag_seq\": '" in deseq2
+
+
+def test_manifest_records_tag_seq():
+    assert '"tag_seq": ${params.tag_seq ?' in _read(MODULES / "emit_manifest.nf")

workflows/rnaseq/README.md

@@ -165,7 +165,18 @@ than they are. A quick way to tell a QuantSeq FWD library from full-length
 RNA-seq when the deposit does not say: the fraction of reads containing 15+
 consecutive A's is several percent for QuantSeq and about 0.1% for
 full-length libraries. `tag_seq` is recorded in the workflow manifest's
-`parameters`.
+`parameters`, and the length model DESeq2 actually used in `deseq2_design`
+(`length_model`).
+
+The length offsets do more than shift expression levels. They can create fold
+changes. Measured on simulated QuantSeq-like reads from the *M.
+galloprovincialis* RefSeq transcriptome (six single-end 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, the standard path called 26
+differentially expressed (median |log2FC| 0.84, up to 3.0). The per-sample
+length offset moves with the isoform mix even when the 3' read count does not.
+With `--tag_seq`, none were called (median |log2FC| 0.09).
 
 ## Covariates (`raw_reanalysis` only)
 
@@ -274,6 +285,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**: executed natively on 2026-10-12, with and without
+  `--tag_seq`, against six simulated QuantSeq-like single-end libraries
+  (82,000 to 101,000 reads) and the full *M. galloprovincialis* RefSeq
+  transcriptome; see the Tag-seq section for the result. Salmon received
+  `--noLengthCorrection` and fastp `--trim_poly_x` in every 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