diff --git a/docs/candidate_studies.md b/docs/candidate_studies.md index 613aff2..fcaf0f5 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 `--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. | diff --git a/modules/rnaseq/differential_expression_deseq2.nf b/modules/rnaseq/differential_expression_deseq2.nf index e7cfbae..8665435 100644 --- a/modules/rnaseq/differential_expression_deseq2.nf +++ b/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), "}" diff --git a/tests/test_rnaseq_tag_seq.py b/tests/test_rnaseq_tag_seq.py new file mode 100644 index 0000000..a643d0a --- /dev/null +++ b/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") diff --git a/workflows/rnaseq/README.md b/workflows/rnaseq/README.md index cc683cb..927af8f 100644 --- a/workflows/rnaseq/README.md +++ b/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