# Reproducibility plan Rebuilding this repository's analysis as a Quarto project that anyone can clone and render end to end. Audit performed 2026-08-29 against commit `94718d8`, tagged `v0-preexisting`. Local toolchain at that point: R 4.3.2, Quarto 1.6.40, DESeq2 1.42.1. **Status:** phases 0-4 are complete and phase 5 is under way; the analysis has rendered and been checked on a clean machine (see the 2026-09-06 entry at the end). --- ## Why this exists The gonad and liver differential-expression results in this repository are real, but at the time of the audit **no notebook in it ran as written**. The upstream half was a shell transcript bound to one machine; the downstream half assumed a working directory it was never rendered from; three of the six notebooks were copied wholesale from an unrelated oyster experiment; and nothing anywhere pinned a package version. Seventeen findings are recorded below, grouped by what they cost. `B*` stops a re-run outright. `R*` lets a re-run finish and hand back something subtly different from what the manuscript reports. `E*` is the environment and repository hygiene that makes both likely. --- ## Decisions These were open questions in the first draft. They are settled and the phases below reflect them. | # | Question | Decision | |---|---|---| | 1 | Finding R2: the outputs claim an `lfc = 1.5` filter that never ran. Apply the filter, or correct the labels? | **`padj < 0.05` only.** Delete the unused `lfc.cutoff` / `svalue.cutoff` variables and correct the misleading labels. Published DEG counts stay as they are. | | 2 | Should `DESEQ_output/` stay in git? | **Yes, keep it tracked.** Results stay browsable on GitHub without cloning. Phase 2 must therefore reproduce them in place rather than into an ignored directory. | | 3 | Is functional enrichment in scope, given the existing results came from the DAVID web tool? | **In scope, and DAVID is dropped.** Replace it with a fully scripted R workflow (see below — the annotation source turned out to be much better than expected). | | 4 | What happened to the brain tissue that the upstream notebook processes but nothing consumes? | **Presume the data is missing.** Do not delete the brain branch; park it, clearly marked, and revisit if the FASTQs are located. | ### On decision 3 — the annotation source The first draft assumed real enrichment meant building a sockeye GO annotation from scratch via BLAST against SwissProt. It does not. NCBI publishes a GO annotation file for this exact assembly: ``` https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/006/149/115/ GCF_006149115.2_Oner_1.1/GCF_006149115.2_Oner_1.1_gene_ontology.gaf.gz ``` Verified 2026-08-29: - 1.6 MB, GAF 2.2, generated by NCBI 2023-12-02, tied to `GCF_006149115.2` — the assembly this project already used. - 155,362 annotations over 32,629 genes (42,524 C / 55,639 F / 57,199 P). - Column 3 is the `LOC…` symbol, which matches our gene IDs **exactly**: all 32,629 GAF genes are present in the count matrix, so 32,629 of 37,942 genes (86%) carry GO terms. 24 of the 31 significant liver DEGs are annotated. So enrichment becomes one `download.file()` plus `goseq`, with gene lengths derived from `sequences/GCF_006149115.2_Oner_1.1_mRNA.gff`, which is already in the repository. No BLAST, no web tool, no manual step. **KEGG — correction.** An earlier draft of this plan said there was no KEGG organism for *Oncorhynchus nerka*. That was wrong. The check queried the KEGG REST organism list, got an empty response, and misread "no data" as "no match". KEGG carries sockeye as organism **`one`** (genome T07952, 37,942 proteins — the same gene count as our matrix), and its gene identifiers are NCBI GeneIDs, which join to our `LOC` ids directly. The giveaway was in the repository all along: `DAVID_KEGG_pathway.txt` contains strings like `one04110:Cell cycle`. Both KEGG REST responses are vendored next to the GAF, with their retrieval date recorded in `tag-seq/genome/kegg_one_SOURCE.txt`. No ortholog transfer needed. --- ## Findings ### Blocking — nothing rendered - **B1 — Three different working-directory assumptions, none declared.** The notebooks live in `tag-seq/code/`. The DESeq2 notebooks read `data/…` and write `DESEQ_output/…`, paths relative to `tag-seq/`, one level up. `3_generate_gene_tables-salmon.Rmd` reads `../DESEQ_output/…`, relative to `code/`. Nothing sets `knitr::opts_knit$set(root.dir=)`. These were run chunk-by-chunk with the console parked at `tag-seq/`; pressing Knit fails on the first read. - **B2 — The gene-table notebook operates on objects that never exist.** It loads `liver` and `gonad`, then manipulates `control_ploidy`, `heat_ploidy`, `diploid_heat` and five more, all leftovers from an oyster ploidy x multistressor experiment. It reads the DESeq2 CSVs with `sep=" "` when notebook 2 wrote them with `sep=","`. Its final chunk does not parse: `geom_text_repel(` is never closed. - **B3 — Both enrichment notebooks belong to a different experiment.** They read `gene_tables/filtered/HISAT2_multiqc_biplot/*.txt` and `blast/geneid2go.tab`, neither of which exists here, and pin their shell chunks to `engine.path="C:/Windows/System32/bash.exe"`. - **B4 — The upstream notebook is a transcript of one session on one machine.** Absolute paths throughout (`/home/shared/8TB_HDD_02/mattgeorgephd/…`). Two steps say "I moved files" by hand. It carries an abandoned bowtie2 branch alongside the HISAT2 branch actually used, and a brain branch (decision 4). ### Silently wrong — runs, but may not reproduce - **R1 — Sample labels are assigned by position, and the check is tautological.** `colnames(cts) <- row.names(trt_list)` overwrites `01B_S43_R1.gtf…` with `B01…` purely by column order. The next line tests `all(colnames(cts) %in% rownames(coldata))`, which is now unconditionally `TRUE`. The order is correct in the committed matrices (verified), but any re-run of `prepDE.py` that reorders columns mislabels every sample silently. **Highest-value fix in this plan.** - **R2 — Reported thresholds do not match the filtering that ran.** See decision 1. - **R3 — The reference level is implicit.** `lfcShrink(coef=2)` relies on `factor()` sorting `social` before `territorial`. Correct today; a renamed level silently flips the sign of every `log2FoldChange`. - **R4 — Heatmap group colours are hardcoded counts.** `ColSideColors = c(rep("royalblue1", 15), rep("red3", 15))` in liver, `14`/`14` in gonad, both assuming column order rather than deriving it from `coldata`. - **R8 — The gene-identifier step destroyed the identity of 4,731 genes.** *(found during phase 2, fully characterised and **fixed** afterwards)* Row names are `gene-|`, where the identifier is a `LOC` number for genes with no assigned symbol and a real symbol for the rest. The `.Rmd` split each name on the literal string `LOC` and kept the third piece, which works for `LOC115144855` and turns `gene-arhgap8|arhgap8` into the string `LOCNA` — for 4,731 genes, 12.5% of the annotation, and disproportionately the well-characterised ones. My first assessment — that these were dropped before reaching any output — was wrong. `data.frame()` de-duplicates row names, so they reached the published tables as `LOCNA.1`, `LOCNA.2292` and so on: identifiers corresponding to nothing, joining to nothing. In the committed results that was 936 rows in `liver-ALL-DEG-apeglm.csv` and 3,141 in the gonad equivalent, including **3 of 31 significant liver genes and 309 of 1,630 significant gonad genes**, all of which the gene tables then discarded as "no description". See *Identifier repair* below. - **R5 — Sample exclusion is a comment toggle with no recorded criterion.** Gonad drops `C05` and `C17`; liver has the same two lines commented out. The only justification in the repo is "MultiQC report: Pheatmap: C05, C17", with no MultiQC output committed to support it. - **R6 — `cache = TRUE` paired with `rm(list=ls())`.** Cached chunks restore objects the clear was meant to remove, and the cache is not invalidated when `data/` changes. - **R7 — `setwd()` inside chunks.** knitr resets the working directory at every chunk boundary, so this works interactively and quietly writes to the wrong place when rendered. ### Environment and repository - **E1 — No versions pinned; missing packages failed open.** *(fixed in phase 1)* - **E2 — Twelve of twenty-six loaded packages were never called.** *(fixed in phase 1)* - **E3 — RStudio session state committed; no `.gitignore`.** *(fixed in phase 0)* - **E4 — Inputs and derived outputs are mixed.** `DESEQ_output/` is 37 MB of regenerable results; `data/` and `genome/` are 79 MB of inputs. Per decision 2 the results stay tracked, but source and derived should still be signposted. - **E5 — The two DESeq2 notebooks are a copy-paste pair.** Normalised for tissue name, 571 and 582 lines differ in five places. - **E6 — The README documents logistics, not the analysis.** - **E7 — `sequences/GCF_006149115.2_Oner_1.1_mRNA.gff` is tracked and zero bytes.** *(new, found during phase 3)* The `sortBed` step in the upstream notebook that should have written it produced nothing that was committed. Gene lengths for `goseq` come from the assembly feature table instead. Phase 4 should either regenerate the file or stop tracking it. --- ## Target layout ``` _quarto.yml # project + render order README.md # how to run it, in five lines renv.lock # pinned R + Bioconductor versions .Rprofile # pins CRAN to a dated snapshot REPRODUCIBILITY-PLAN.md # this file index.qmd # landing page: design, samples, results index code/ 01-upstream-alignment.qmd # documented, parameterised, eval: false 02-differential-expression.qmd # params: tissue -> renders twice 03-gene-tables.qmd 04-enrichment.qmd # goseq against the NCBI GAF _common.R # theme, paths, load_counts(), assertions data/ # SOURCE ONLY - never written by code treatments-{liver,gonad}.csv onerka_gene_count_matrix-{liver,gonad}.csv genome/Onerka_LOCID_gene_table.txt genome/GCF_006149115.2_Oner_1.1_gene_ontology.gaf.gz DESEQ_output/ # DERIVED but tracked (decision 2) liver/ gonad/ docs/ # rendered site, published via GitHub Pages ``` --- ## Phases ### Phase 0 — Stop the bleeding  `DONE` (commit `66ff8b6`) Tagged `v0-preexisting` at `94718d8` first, so the outputs the manuscript cites stay retrievable. Then removed 217 files of editor session state from the index (`.Rproj.user/`, `.RData`, two `.Rhistory`), added a `.gitignore`, and set `RestoreWorkspace` / `SaveWorkspace` / `AlwaysSaveHistory` to `No` in the `.Rproj` so `.RData` cannot come back. Encoding `ISO8859-1` -> `UTF-8`. Working copies untouched; history not rewritten, so existing clones and links keep working. 299 tracked files -> 84. (The first draft estimated "under 60"; the remainder is real data, code and committed results.) ### Phase 1 — Pin the environment  `DONE` (commit `de3c217`) `renv.lock` now records 148 packages against R 4.3.2 and Bioconductor 3.18, including DESeq2 1.42.1, apeglm 1.24.0, PCAtools 2.14.0, ashr 2.2-63. Pinning was not the formality it looked like: - `renv::init()` defaults CRAN to `packagemanager.posit.co/cran/latest`. Against latest, `Deriv` 4.3.0 — pulled in transitively by `ggpubr -> rstatix -> car -> pbkrtest -> doBy` — fails to compile on R 4.3.2, because its sources now target R >= 4.5. A lockfile that cannot restore is not a pin. Fixed with a dated snapshot, `cran/2024-04-24`, contemporary with Bioconductor 3.18. *(Superseded 2026-09-27: the lockfile was later snapshotted with packages that snapshot does not hold, and the pin is now `cran/2025-05-18`. See "The pinned snapshot now restores the lockfile" below.)* - The pin must be applied **after** `renv/activate.R`, which restores the lockfile's repositories and would otherwise overwrite it; and it must **amend** the `CRAN` entry rather than replace the vector, which silently drops all five Bioconductor repositories and makes DESeq2 unresolvable. Both are commented in `.Rprofile`. - renv itself is pinned to 1.0.7, the version in that snapshot, so a clean clone can bootstrap from the same place. - The library was rebuilt entirely from the pinned repositories rather than hydrated from the author's user library, which had produced a mixed set (a newer `ggpubr` demanding a `cowplot` the snapshot cannot supply). Both DESeq2 notebooks now use explicit `library()` calls. Nothing installs at render time, and a missing package stops the render — the old `sapply(load.lib, require)` returned `FALSE` and carried on. Fifteen of the twenty-six loaded packages were dropped after verifying none of their functions are called. `tag-seq/code/.renvignore` excludes the three oyster-project notebooks from dependency discovery; `3_generate_gene_tables` additionally requests `kmisc`, which is not a CRAN package and would fail every restore. Phase 3 removes that exclusion for the rewritten file. **Verified:** `renv::status()` reports a consistent project, all ten analysis packages load from the project library, and the rewritten `load_libraries` chunk executes without error. ### Phase 2 — Port the DESeq2 analysis to one parameterised Quarto document  `DONE` Two 570-line near-duplicates become one document rendered twice. - `_quarto.yml` with the project root as execution directory, so every path is written once and means one thing. Delete every `setwd()`. *(B1, R7)* - `02-differential-expression.qmd` takes `params: tissue`, rendered with `liver` and `gonad`. Axis limits, group sizes and excluded samples become parameters. *(E5)* - Replace the positional column rename with a real join on a sample-ID column, and `stopifnot()` that the match is complete and one-to-one. *(R1)* - `relevel(trt, ref = "social")` and `stopifnot(resultsNames(dds)[2] == "trt_territorial_vs_social")`. *(R3)* - Delete the unused `lfc.cutoff` / `svalue.cutoff` and correct the labels that claim them, per decision 1. *(R2)* - Derive all heatmap annotation from `coldata`. *(R4)* - Move sample exclusion into a parameter with a written criterion, and commit the MultiQC output that justifies it. *(R5)* - Replace `cache: true` with Quarto's `freeze: auto`, which invalidates on source change. Remove `rm(list=ls())`. *(R6)* - Shared theme and loader into `code/_common.R`. `_quarto.yml` at the project root sets `execute-dir: project` and `freeze: auto`. `tag-seq/code/_common.R` holds the paths, the per-tissue configuration, the loader and the theme; `tag-seq/code/02-differential-expression.qmd` is 386 lines and replaces the two notebooks' 1,153. Render both with `tag-seq/code/render-differential-expression.sh`. **Verification.** Both tissues were rendered and the results diffed against the committed baseline. Every table reproduces: | table | rows old / new | same gene set | max abs d(log2FC) | max abs d(padj) | |---|---|---|---|---| | liver ALL-DEG (unshrunken) | 4076 / 4076 | yes | 4.4e-09 | 1.3e-07 | | liver ALL-DEG normal | 5646 / 5646 | yes | 2.9e-08 | 1.4e-07 | | liver ALL-DEG apeglm | 5646 / 5646 | yes | **3.0e-02** | 1.4e-07 | | liver ALL-DEG ashr | 5646 / 5646 | yes | 5.9e-08 | 1.4e-07 | | liver SIG-DEG unshrunken | 66 / 66 | yes | 1.9e-14 | 6.6e-13 | | liver SIG-DEG apeglm | 31 / 31 | yes | **1.1e-03** | 5.6e-13 | | gonad ALL-DEG apeglm | 18935 / 18935 | yes | **4.2e-02** | 1.7e-07 | | gonad SIG-DEG unshrunken | 1653 / 1653 | yes | 3.8e-10 | 2.1e-11 | | gonad SIG-DEG apeglm | 1630 / 1630 | yes | **5.3e-03** | 2.1e-11 | Gene sets, row counts, descriptions and sort order are identical throughout, and `baseMean` agrees to 1e-15 relative — the counts, normalisation and filtering are reproduced exactly. `apeglm` is the only estimator whose fold changes move at all, in the third decimal; its shrinkage is an iterative optimisation and is sensitive to the package version. No gene changes significance status. Two intended differences: - `liver-gene-counts.csv` row labels change per decision 1 (`DEG_apeglm-s0.005_lfc1.5` -> `DEG_apeglm-p0.05`). **Every count is unchanged** (66 / 31 / 31 / 31). `GONAD-gene-counts.csv` is byte-identical, since its labels were already honest. - Image files differ byte-wise. They are re-rendered by a different ggplot2 and R graphics stack; the content is equivalent. `renv::status()` remains clean. `shiny` is marked ignored in `renv/settings.json`: renv flags it for any document with a `params:` header, but params are always passed on the command line, never prompted for interactively. ### Phase 3 — Rebuild the gene tables, and script the enrichment  `DONE` - Rewrite `03-gene-tables.qmd` against the two tissues that exist. The characterized / uncharacterized / `LOCNA` split is about fifteen lines of dplyr in a function called twice, not 165 lines of oyster contrasts. *(B2)* - Pass tibbles directly rather than round-tripping through CSV with mismatched delimiters. - Delete `X_GOterm_Annotation.Rmd` and `X_GOSeq_annotation.Rmd`. *(B3)* - New `04-enrichment.qmd`: download the assembly GAF, build the `gene2GO` mapping, compute gene lengths from `sequences/GCF_006149115.2_Oner_1.1_mRNA.gff`, and run `goseq` with the length-bias correction. Commit the GAF alongside the other inputs so the analysis does not depend on NCBI staying up. Report GO only; KEGG needs ortholog transfer and is out of scope. - Drop the rewritten gene-table notebook from `.renvignore` so its dependencies get pinned; add `goseq` to the lockfile. **Gene tables.** `03-gene-tables.qmd` replaces the 165-line oyster notebook. Whatever produced the committed `gene_tables/` was never in the repository, so the rule was reconstructed from the outputs and verified: drop genes with no description (the old code's `LOCNA` class), then split the rest on whether the description begins with `uncharacterized`. Reproduces exactly — liver 28 / 0 / 28, gonad 1280 / 41 / 1321 — and `uncharacterized_stats.txt` is byte-identical. One deliberate format change: the committed tables were written with `quote = FALSE, sep = ","` while descriptions contain commas (`uncharacterized LOC115129809, transcript variant X3`). In `gonad_characterized.txt` only 754 of 1,280 rows had the expected 8 fields; 526 had 9-11 and could not be read back. They are now properly quoted CSV. Rows and values are unchanged. `gonad_all_LOCID.txt` also reorders 71 of 1,321 lines — adjacent pairs with near-tied fold changes, swapped by apeglm's third-decimal shift from phase 2. Same gene set. **Enrichment.** `04-enrichment.qmd` runs `goseq` — with its length-bias correction — for GO and KEGG, both tissues. Universe is every gene with an adjusted p-value; DE is `padj < 0.05`, the same rule the results tables use. Gene lengths could not come from `sequences/GCF_006149115.2_Oner_1.1_mRNA.gff` as planned: **that file is tracked and zero bytes** (new finding E7). They come from the assembly feature table instead — median mRNA interval length per gene, keyed by `GeneID`, which is our `LOC` id without the prefix. 37,929 of 37,942 genes get a length. Results, written to `tag-seq/GO_output//`: | tissue | categories tested | FDR < 0.05 | |---|---|---| | liver GO | 3275 | 0 | | liver KEGG | 187 | 0 | | gonad GO | 5209 | 1 | | gonad KEGG | 194 | 0 | The single hit is `GO:0051131` *chaperone-mediated protein complex assembly* (BP, 6 of 6 genes DE, FDR 0.028). Liver has only 31 significant genes, so no enrichment there is the expected result rather than a failure. **On the DAVID files.** They are not enrichment results and never were. `DAVID_GOterms.txt` is a 53-row gene report (symbol, name, species) with no GO terms in it at all, and `DAVID_KEGG_pathway.txt` is a 13-row gene report listing each gene's pathways. There was no committed enrichment analysis to reproduce, so these numbers are the project's first. **Done:** `gene_tables/` and the enrichment results regenerate from a render, and every notebook in the repository either runs or is gone. ### Phase 4 — Make the upstream steps honest  `DONE` ### Identifier repair  `DONE` Fixes R8. Identifiers are now taken from after the pipe, which is unique across all 37,942 genes. Descriptions and lengths come from the assembly feature table instead of `Onerka_LOCID_gene_table.txt`, which is keyed entirely by `LOC` identifiers and so could never annotate a symbol-named gene. The feature table covers 37,929 genes and **agrees with the old table on all 33,210 they share** — verified, so this is a strict superset, not a different annotation. KEGG's GeneID keys are mapped through the same table rather than by pasting a `LOC` prefix, which would have been right for LOC-numbered genes and wrong for the rest. **No published statistic changed.** Row names never entered the model, so `baseMean`, `log2FoldChange` and `padj` are bit-identical to the previous render in both tissues — maximum difference 0.00e+00. What changed is which genes can be identified: | | before | after | |---|---|---| | fabricated `LOCNA.*` ids, liver / gonad ALL-DEG | 936 / 3141 | 0 / 0 | | rows with no description, liver / gonad | 936 / 3142 | 11 / 12 | | gene tables, liver (char / unchar / total) | 28 / 0 / 28 | **31 / 0 / 31** | | gene tables, gonad | 1280 / 41 / 1321 | **1587 / 41 / 1628** | | GO-annotated universe, liver / gonad | ~4710 / ~15794 | **5259 / 17086** | | GO categories tested, liver / gonad | 3275 / 5209 | 3749 / 6066 | | significant GO terms, gonad | 1 | **2** | | significant KEGG pathways, gonad | 0 | **1** | The recovered biology is not marginal. Gonad KEGG now returns **`one00190` Oxidative phosphorylation, 59 of 188 pathway genes differentially expressed, p = 2.3e-05, FDR = 0.0044** — invisible before precisely because oxidative phosphorylation genes are well characterised, carry real symbols, and were therefore the ones being discarded. GO adds *negative regulation of non-canonical Wnt signaling pathway* alongside the existing *chaperone-mediated protein complex assembly*. **Completed:** `01-upstream-alignment.qmd` replaces the 640-line shell-transcript `.Rmd` with: - All paths parameterized at the top and sourced from a single configuration block - Bowtie2 branch removed; HISAT2 only - Brain tissue clearly marked as unresolved with decision 4 rationale - All "I moved files" comments replaced with actual shell commands - Heavy chunks marked with `eval: false` and documented as Raven-only - Tool versions documented in comments (HISAT2 2.2.1, StringTie 2.2.1, samtools 1.12, cutadapt, FastQC) - Alignment rates and sample-exclusion justification committed in `tag-seq/QC/ALIGNMENT_SUMMARY.md` - Old notebook deprecated but retained for reference ### Phase 5 — Publish and keep it honest - Render to `docs/`, enable GitHub Pages. Every figure and table becomes a URL for the manuscript's data-availability statement. - Rewrite the README: what the experiment was, how to run it (`renv::restore()`, then `quarto render`), what each notebook does, where the raw data lives on Gannet. Keep the existing logistics links. *(E6)* - GitHub Action running `renv::restore()` and `quarto render` on push. It will not run the FASTQ steps, but it proves the R analysis still renders from a clean checkout — exactly the property that decayed here. - Archive to Zenodo for a citable DOI at submission. --- ## Post-phase-4 audit (2026-08-30) An independent audit re-checked the phase 0-4 claims by recomputing the committed result tables rather than re-rendering. Everything checkable held: zero fabricated `LOCNA.*` identifiers remain in any of the eight `ALL-DEG` tables; missing descriptions are 11 of 5,646 (liver) and 12 of 18,935 (gonad); the gene tables re-derive exactly from `*-SIG-DEG-apeglm.csv` (liver 31/0, gonad 1,587/41, 1,628 LOCIDs, identical gene sets); recomputed `padj < 0.05` counts match the committed `*-gene-counts.csv` exactly; and all 32,629 genes in the vendored GAF are present in the 37,942-gene count matrix. It also found twelve open items. Three have been closed: **Step 1 — one command means the whole analysis** *(B5, B6, E8)* - `_quarto.yml` now declares an explicit `project.render` list. Without it, Quarto rendered every `.qmd` and `.Rmd` in the project, including the deprecated upstream transcript — ~40 `{bash}` chunks with no `eval` guard on any of them, pointing at `/home/shared/` and Gannet *(B5)*. - That transcript moved to `archive/`, out of render scope, with a `knitr::opts_chunk$set(eval = FALSE)` guard added as a second line of defence. - `render-all.sh` at the project root is now the documented entry point. A bare `quarto render` renders the parameterised notebook 02 once, with its default tissue, so it regenerated liver and silently left gonad alone *(B6)*. The script renders 02 per tissue via the existing `tag-seq/code/render-differential-expression.sh`, then 03, then 04, and prints the `git diff --stat` command that constitutes the reproduction check. - `tag-seq.Rproj` moved to the project root as `project-sockeye-tagseq.Rproj`. In `tag-seq/` it set the interactive working directory one level below what `execute-dir: project` assumes, so every `source("tag-seq/code/_common.R")` failed when chunks were run by hand *(E8)*. - README rewritten: `./render-all.sh` not `quarto render`, a table of what each notebook reads and writes, and the expected-diff statement. **Step 4 — derived-file bookkeeping** *(E7, E11, E13)* - `sequences/GCF_006149115.2_Oner_1.1_mRNA.gff` untracked. It was a zero-byte file asserting that a derived artefact existed; it is regenerated by the `bedtools` chunk in 01 on Raven, and nothing downstream reads it — goseq's lengths come from the feature table. Recorded in a new `tag-seq/sequences/README.md`, with the prose in `_common.R` and 04 corrected *(E7)*. - The DAVID files moved from `tag-seq/DESEQ_output/gonad/` — where a render overwrites its neighbours — to `archive/david-2023/`, with `archive/README.md` stating why they are not maintained: they cannot be reproduced from this repository, and both are per-gene reports (columns: gene symbol, name, species, pathway) with no over-representation statistic and no FDR anywhere in either file *(E13)*. - Signpost READMEs in `tag-seq/data/`, `tag-seq/genome/`, `tag-seq/sequences/` and `figures/` classify every tracked file as input, provenance-only or superseded, and name the notebook that reads it. The 33 MB of provenance-only files (both transcript matrices, `onerka_merged-liver.gtf`) and the 3.8 MB superseded `Onerka_LOCID_gene_table.txt` are now labelled as such, with keep-or-move recorded as an open decision rather than resolved by deletion. `figures/figure_X.*` is labelled manuscript-only, produced by no script *(E11)*. **Step 5 — input integrity** *(E12, part of)* - `CHECKSUMS.sha256` records sha256 for the eight files the analysis reads: both count matrices, both treatment tables, the feature table, the GAF, and the two vendored KEGG REST responses. Verified against the current working tree. - `verify_inputs()` in `_common.R` checks the manifest whenever the file is sourced, i.e. at the top of 02, 03 and 04. A mismatch warns rather than stops — a deliberate data update should be loud, not blocking — and the check skips with a message if `digest` is unavailable, so no notebook gains a hard dependency. This is finding R1's failure mode one layer out: R1 was a count matrix whose columns had been silently reordered; this catches the input file changing at all. - Still open under E12: `LICENSE` and `CITATION.cff`. **Steps 6 and 7 — publish, and check it from outside** *(E6, P5)* - `index.qmd` at the project root is the landing page: design, sample table, results as committed, links to the five rendered notebooks and the tracked result directories, and a data-availability statement. It carries the C05/C17 caveat as a callout, because a reader who lands on the site should meet that before they read the gonad numbers. - `output-dir: docs`, and `/docs/` is gitignored. Nothing generated is committed, so a stale rendered page cannot contradict the tracked tables. The Action publishes `docs/` to Pages as a build artifact instead, which needs Pages set to "GitHub Actions" as its source. - `check-reproduction.R` replaces `git diff --exit-code` as the reproduction test. A bare diff cannot work here: re-rendered PNG bytes always differ, and apeglm's iterative shrinkage moves `log2FoldChange` in the third decimal across package versions — exactly the columns 02, 03 and 04 write. So images are reported and ignored, tables are compared on matched keys with a relative tolerance (1e-3, and 1e-2 for p-value columns, which span thirty orders of magnitude), and significance calls at 0.05 must be identical with no tolerance at all. Files it cannot parse count as failures, not passes. - `.github/workflows/render.yml` restores `renv.lock` on a clean Ubuntu runner, verifies `CHECKSUMS.sha256`, runs `./render-all.sh`, then `check-reproduction.R`; it uploads the rendered site and, on failure, the result diff. First run builds the Bioconductor stack from source (budget an hour); later runs restore from cache. `check-reproduction.R` was exercised against deliberately perturbed result tables before being committed: fold changes moved 0.05% pass, 5% fail, a single padj nudged across 0.05 fails on the significance call even though the numeric move is 0.6%, a dropped gene fails on dimensions, and appended bytes in a PNG pass as an ignored image. It and the `verify_inputs()` addition to `_common.R` both parse under R 4.5.3. **Still open:** R5/E9 (the C05/C17 exclusion has no committed evidence, and `ALIGNMENT_SUMMARY.md` cites MultiQC directories containing only `.gitkeep`, attributes alignment statistics to uncommitted logs, and reports n=15 per tissue where both matrices have 30 columns), step 3 (no independent re-render has been performed), E10 (no SRA accession), E14 (upstream toolchain pinned only in prose), and all of phase 5 below. --- ## First clean-machine render (2026-09-06, PR #4) Neither run of the render workflow on 2026-08-30 had got past its first step: notebook 01 has only bash chunks, Quarto inferred the Jupyter engine, and the runner had no `nbformat`. With `engine: knitr` declared, and the parameterised 02 renders named through `output-file` rather than `--output` (which Quarto 1.10 mishandled; CI now pins 1.6.40), the analysis rendered end to end on Ubuntu for the first time. **What reproduced.** Every DESeq2 table: same rows, same `padj < 0.05` calls in both tissues. baseMean and normalised counts agree to 1e-14 relative; apeglm `log2FoldChange` to 0.0008 and `lfcSE` to 0.0013 log2 units. Step 3 of the post-phase-4 audit, "no independent re-render has been performed", is closed. **What did not.** The gonad enrichment. Gene-to-category counts were identical, but goseq's Wallenius p-values moved because the length-bias spline fitted by `nullp()` is platform-sensitive, and the three enriched terms sat at the threshold: FDR 0.004 / 0.039 / 0.045 on macOS, 0.060 / 0.057 / 0.081 on Linux. Investigating why the weighting mattered so much showed that it should not have been applied at all. TagSeq is 3′-end counting, so the length bias goseq corrects for is not present by design, and the length association in these data is strong and negative (gonad: 27% of the shortest-decile genes DE against about 5% in every other decile). The test is now hypergeometric, exact and deterministic. Gonad returns 35 GO terms and 6 KEGG pathways at FDR < 0.05, ribosomal and mitochondrial; liver still none. The Wallenius tables are in the history at `bf5bf06` and before. **Checker.** Ten of the fourteen failures the run reported were the checker's: a non-unique key on the single-gene-count tables, no absolute floor under the relative tolerance, and 0.05 call checks on raw p-values. All three fixed and verified against the Linux outputs and against perturbation probes. Also closed from the earlier list: E12's `LICENSE` and `CITATION.cff` remain open; the eleven count tables changed to comma-separated by `8ea8f80` but never regenerated are now regenerated and committed; GitHub Pages is enabled with the Actions source. --- ## One significance convention (2026-09-26, submission checklist B1) **The problem.** Notebook 02 computed the unshrunken results with `results(..., alpha = 0.05)` but called `lfcShrink()` without `res =`, so each shrinkage estimator computed its own results internally at DESeq2's default `alpha = 0.1`. `alpha` is the target independent filtering optimises for, so the two conventions filtered different low-count genes and produced different adjusted p-values. At the same `padj < 0.05` threshold the unshrunken table called 66 liver and 1,653 gonad genes and the shrunken tables 31 and 1,630. Every downstream product read the shrunken tables. The draft manuscript carried a Limitations paragraph explaining the discrepancy. **The fix.** The alpha-0.05 results object is passed to all three shrinkage calls through `res =`, and the notebook asserts that every shrunken table carries exactly the same `padj` as the unshrunken one. Shrinkage now changes fold changes and their standard errors only. **Re-render.** Run on Hyak in `bioconductor/bioconductor_docker:3.18-R-4.3.2` with Quarto 1.6.40, restoring `renv.lock` (R 4.3.2, DESeq2 1.42.1, apeglm 1.24.0, goseq 1.54.0); all eight input checksums matched. The unshrunken tables came back unchanged, which confirms the environment reproduces CI. The new tables were committed as the baseline for `check-reproduction.R`. | | Before | After | |---|---|---| | Liver genes tested / significant | 5,646 / 31 | 4,076 / 66 | | Gonad genes tested / significant | 18,935 / 1,630 | 18,156 / 1,653 | | Gonad KEGG / GO enriched (FDR < 0.05) | 6 / 35 | 5 / 36 | | Liver KEGG / GO enriched | 0 / 0 | 0 / 0 | | Genes significant in both tissues | 11 | 19 | Liver lost 3 genes (they are no longer tested under the stricter filtering; two were named in the draft, coagulation factor X and MARCKS-related protein) and gained 38, mostly respiratory-chain, ATP-synthase and ribosomal subunits. Gonad lost 14 and gained 37, none with |log2FC| above 0.93; the 20 largest effects are identical and in the same order. KEGG RNA polymerase (`one03020`) fell from FDR 0.019 to 0.075. All 19 shared genes change in the same direction in both tissues. Every hand-typed number in `index.qmd` and the manuscript was updated or replaced with an inline computation, and the Limitations paragraph was removed. `figures/figure_X` is hand-built and still shows the old counts; see `figures/README.md`. **Wallenius, again.** On the new gonad gene set goseq's `nullp()` spline fit fails outright: every fitted probability is the same constant, so a Wallenius test would silently equal the hypergeometric one. On the old gene set in the same environment the fit is sensible and Wallenius returns no KEGG pathway at FDR < 0.05. The decision to use the hypergeometric test stands, with more evidence behind it; notebook 04 and the manuscript now say this rather than quoting the old-gene-set figures as current. **New finding: the CRAN snapshot does not cover the lockfile.** `.Rprofile` pins CRAN to the 2024-04-24 Posit Package Manager snapshot, but `renv.lock` records many packages from a 2025-05-18 snapshot, and some (`nlme` 3.1-166) exist only in the later one. A restore limited to the pinned snapshot fails on `nlme` and `mgcv`. CI succeeds only because its runner also configures "latest" Package Manager as a fallback. The re-render above passed both snapshots to `renv::restore(repos = ...)`. Closed the next day; see below. ## The pinned snapshot now restores the lockfile (2026-09-27, checklist B4) Two ways to make the pinned snapshot and the lockfile agree were considered. Re-resolving the lockfile against 2024-04-24 would have downgraded 77 of the 146 locked CRAN packages, since that many were newer than anything the snapshot held, and forced another result rebaseline. Moving the pin to 2025-05-18 changes no package version: 137 locked versions are current in that snapshot and the other 9 (MASS, Matrix, XML, lattice, mgcv, nlme, nnet, renv, survival) are in its Archive, each confirmed by URL. **Changes.** - `renv.lock`: the CRAN repository URL moves to `cran/2025-05-18`. The 49 packages recorded against that snapshot's URL, and `renv` recorded as `RSPM`, are relabelled `CRAN`, which now names the same repository. No version changed; the diff is those 51 lines. - `.Rprofile`: the pin moves to the same date, and the comment explains why it must match the lockfile. - `.github/workflows/render.yml`: `use-public-rspm: false`, so CI no longer adds "latest" Package Manager as a second repository, which is what hid the mismatch. The renv cache key gains a `v2` prefix: a warm cache installs from the cache and never contacts a repository, so without the bump CI would not have exercised the change. **Verification.** A copy of the project containing only `renv.lock`, `.Rprofile` and `renv/`, with an empty renv cache and no library, ran a bare `renv::restore(prompt = FALSE)` in the same `bioconductor_docker:3.18-R-4.3.2` container whose restore had failed on `nlme` before. The only repositories the session saw were the six Bioconductor 3.18 repositories and `cran/2025-05-18`. It completed in 16 minutes, and all 184 locked packages installed at exactly their locked versions. Because no version changed, the committed results stand without a re-render. **Caveat, unchanged.** The date now governs anything newly installed, and 130 packages current in that snapshot require R >= 4.4. A future `renv::install()` may need an explicit older version; `.Rprofile` says so. --- ## Scripted manuscript figures (2026-09-27, checklist C2) `tag-seq/code/05-figures.qmd` draws Figures 1–3 and replaces the hand-built `figures/figure_X` (finding E11), which is removed. The model fit moved into `fit_deseq()` in `_common.R`, called by both 02 and 05, so the figures are drawn from the same fitted object that produces the tables. A full `render-all.sh` after that refactor passed `check-reproduction.R` with every table unchanged. The notebook's PCA reproduces notebook 02's variance explained exactly (liver 13.82% and 10.86%, gonad 21.37%), and the manuscript now reads those numbers from `figures/source-data/pca-variance.csv` instead of typing them. --- ## Supplementary tables (2026-09-27, checklist C3) `tag-seq/code/06-supplementary-tables.qmd` writes Tables S1–S6 to `manuscript/supplementary/`. It computes no statistics; it joins the committed tables and asserts that the four shrinkage tables carry identical p-values and adjusted p-values, that the significant sets match the `SIG-DEG` tables, and that the shared-gene set matches Figure 3 B's source data. **Bearing on R5/E9 (the C05/C17 exclusion).** Table S1 lists every library, including the excluded ones. C05 and C17 are the two smallest gonad libraries: 1.97 and 2.14 million gene-assigned counts, against 3.39–5.77 million for the 28 retained. This is the first evidence for the exclusion that can be computed from the repository. It is not the MultiQC correlation evidence the decision was originally based on, which is still absent, so R5/E9 stays open; but a library-size criterion is now available if the authors choose to state the exclusion that way. --- ## Tracked outputs trimmed (2026-09-27, checklist D2) `tag-seq/DESEQ_output/` tracked 53 files and 37 MB; it now tracks 20 and 11 MB. What left git and why: - normal and ashr tables (all and significant genes): regenerated into the untracked `alternatives/` folder. Their fold changes are published in Tables S2 and S3, which notebook 06 builds from that folder. - four per-estimator tables of normalised counts per significant gene: after B1 the four estimators flag the same genes, and the four tables were shown to hold identical rows in different orders. One is regenerated into `alternatives/`. - seven of eight volcano plots per tissue, the PCA pairs plot, and the MA plots as PDF (2-3 MB of vector points each); the MA plots are now PNG. - `gonad/MA_plots.png`, committed in 2022 and written by no code. `check-reproduction.R` now also checks `manuscript/supplementary/`, keying rows on tissue as well as gene or category, and skips Markdown files, so dropping the normal and ashr tables from git does not drop them from verification. --- ## Unread inputs out of the tree (2026-09-27, checklist D3; closes E11) The last open decisions of finding E11 are closed. The two transcript-level count matrices, the liver merged GTF (33 MB together) and the superseded `Onerka_LOCID_gene_table.txt` (3.8 MB) are removed from the working tree. None was read by any notebook; a full render without them reproduced every tracked result byte for byte. Each is recorded in `tag-seq/data/README.md` or `tag-seq/genome/README.md` with its sha256 and a commit-pinned GitHub address, downloaded and checked before removal; the LOC table is also on Gannet (CRLF line endings, otherwise identical). Its replacement's superset claim was re-verified: all 33,211 of its genes are in the feature table with matching descriptions (the earlier note said 33,210). Tracked content is now 67 MB, from 129 MB before D2 and D3. Also found: `index.qmd` gave the raw-data location as a Gannet `seashell/bu-github/` address that returns 404; it now gives the `panopea/berdahl-sockeye-salmon/` folder, which holds all 95 FASTQs (gonad 60, liver 30, brain 5) and which every other document already cited. --- ## Notes moved from the notebooks and READMEs (2026-09-27, checklist D4) The notebooks, `_common.R` and the directory READMEs now describe the analysis as it is. The audit narrative they carried (comparisons with the deleted `.Rmd` code, how each finding was found and fixed, and dated history) was removed from them and is kept here verbatim, so nothing that was written down is lost. Each passage is labelled with the finding or checklist item it belongs to and the place it came from; the findings themselves are described in the Findings section above. **E5, D1 (_common.R, tissue_config)** ```text ## Everything that differed between the two near-identical DESeq2 notebooks ## (finding E5) is collected here. Every output file is prefixed with the ## tissue name in lower case (liver-PCA.png, gonad-SIG-DEG-apeglm.csv, ...). ## Gonad files were "GONAD-*" until 2026-09-27, a historical accident of the two ## original notebooks; they were renamed so one rule covers both tissues ## (submission checklist D1). ``` **R5 (_common.R, gonad exclude). Note: the claim that phase 4 committed the MultiQC output was never true; corrected in the new comment** ```text ## C05 and C17 were dropped on the basis of the MultiQC sample-correlation ## heatmap. In the .Rmd this was an uncommented pair of lines in gonad and ## the same pair commented out in liver, with no recorded criterion ## (finding R5). Phase 4 commits the MultiQC output that justifies it. ``` **R1 (_common.R, sample_ids_from_columns)** ```text ## Derive the treatment-table sample ID from a StringTie count-matrix column. ## ## prepDE.py names columns after the GTF files: "01B_S43_R1.gtf". The treatment ## tables use "B01". The mapping is: take the token before the first underscore ## ("01B") and swap its number and tissue letter ("B01"). ## ## The .Rmd did this positionally -- colnames(cts) <- row.names(trt_list) -- and ## then "checked" the result with all(colnames(cts) %in% rownames(coldata)), ## which is unconditionally TRUE after the rename (finding R1). Any prepDE.py ## re-run that emitted columns in a different order would have silently ## mislabelled every sample. ``` **R8 (_common.R, gene_ids_from_rownames)** ```text ## Take the gene identifier from a StringTie row name. ## ## Row names are "gene-|": the identifier is a LOC number for genes ## with no assigned symbol ("gene-LOC115144855|LOC115144855") and a real symbol ## for the rest ("gene-arhgap8|arhgap8"). The part after the pipe is that ## identifier in both cases, and it is unique across all 37,942 genes. ## ## The .Rmd instead split each name on the literal "LOC" and kept the third ## piece, which works for LOC-numbered genes and produces the string "LOCNA" ## for every symbol-named one -- 4,731 genes, 12.5% of the annotation. Because ## data.frame() then de-duplicates row names, those genes reached the published ## tables as LOCNA.1, LOCNA.2292 and so on: identifiers that correspond to ## nothing and join to nothing, and that the gene tables consequently dropped. ## In the committed results that was 3 of 31 significant liver genes and 309 of ## 1,630 significant gonad genes (finding R8). ``` **R1 (_common.R, load_counts check.names)** ```text ## check.names = FALSE: the StringTie columns start with a digit ## ("01B_S43_R1.gtf"), and read.csv would otherwise silently rename them to ## "X01B_S43_R1.gtf". The .Rmd never noticed because it overwrote the column ## names positionally before ever reading them. ``` **R1 (_common.R, load_counts assertions)** ```text ## The three assertions the old tautological check should have been. ``` **R3, D1 (_common.R, reference level; trt renamed phenotype)** ```text ## Behavioural phenotype, observed at capture: territorial is compared with ## social, the reference. Made explicit so the sign of every log2FoldChange no ## longer depends on "social" happening to sort before "territorial" ## alphabetically (finding R3). The column was called "trt" until 2026-09-27; ## it is not a treatment, since phenotype was observed, not assigned. ``` **R8, E7 (_common.R, load_gene_annotation; LOC table and zero-byte GFF)** ```text ## Built from the assembly feature table rather than Onerka_LOCID_gene_table.txt, ## which is keyed entirely by LOC identifiers (33,210 of our 37,942 genes) and so ## could never annotate a symbol-named gene. The feature table covers 37,929 and ## agrees with the old table on every one of the 33,210 they share -- verified, ## so this is a strict superset, not a different annotation. ## ## Length is the median mRNA interval per gene. goseq's nullp() takes it as ## bias data; the test itself is hypergeometric (see 04), so the length enters ## the diagnostic plot but not the p-values. It comes from here rather than from ## sequences/GCF_006149115.2_Oner_1.1_mRNA.gff, which was tracked at zero bytes ## and has since been untracked (finding E7); the feature table is committed and ## complete, so no upstream step has to be re-run to get lengths. ``` **R3 (_common.R, fit_deseq)** ```text ## social is the reference level (set in load_counts), so coefficient 2 is the ## territorial-vs-social effect; that is asserted, not assumed (finding R3). ``` **R4 (_common.R, colours)** ```text ## Phenotype colours, derived from coldata rather than hardcoded run lengths. ## The .Rmd wrote rep("royalblue1", 15) / rep("red3", 15) for liver and 14 / 14 ## for gonad, assuming both the group sizes and the column order (finding R4). ## ## The hues ``` **E12, R1 (_common.R, verify_inputs)** ```text ## Check the input files against CHECKSUMS.sha256 (finding E12). ## ## Finding R1 was a count matrix whose columns had been reordered by a re-run of ## prepDE.py, silently mislabelling every sample; the loader assertions above ## catch that case now. This is the same failure mode one layer further out: an ## input file that changed without anyone noticing, so that the committed ## results no longer correspond to the committed inputs. ## ``` **E5 (02, introduction)** ```text This notebook runs the DESeq2 analysis for **one tissue**, chosen by the `tissue` parameter. It replaces `2_DESeq2_analysis-salmon-liver.Rmd` and `2_DESeq2_analysis-salmon-gonad.Rmd`, which were 571 and 582 lines differing in five places (finding E5). ``` **R1 (02, load data)** ```text Samples are matched between the count matrix and the treatment table **by name**, not by column position. `prepDE.py` names columns after the GTF files (`01B_S43_R1.gtf`); the treatment tables use `B01`. `load_counts()` derives one from the other and asserts the two sets are identical before aligning them. The previous code assigned `colnames(cts) <- row.names(trt_list)` positionally and then "checked" it with `all(colnames(cts) %in% rownames(coldata))`, which is unconditionally `TRUE` after such a rename (finding R1). The assignment happened to be correct for the committed matrices — verified — but nothing would have caught a reordered re-run. ``` **R8 (02, callout 'Gene identifiers were repaired')** ```text ::: {.callout-important} ## Gene identifiers were repaired StringTie row names are `gene-|`, and the identifier after the pipe is a `LOC` number for genes with no assigned symbol and a real symbol for the rest. The `.Rmd` split each name on the literal string `LOC` and kept the third piece, which works for `LOC115144855` and turns `gene-arhgap8|arhgap8` into the string `LOCNA` — for **4,731 genes, 12.5% of the annotation**, and disproportionately the well-characterised ones. Because `data.frame()` de-duplicates row names, those genes did not disappear: they reached the published tables as `LOCNA.1`, `LOCNA.2292` and so on — identifiers that correspond to nothing, join to nothing, and were therefore dropped from the gene tables as "no description". In the committed results that was 3 of 31 significant liver genes and **309 of 1,630 significant gonad genes**. Identifiers are now taken from after the pipe, which is unique across all 37,942 genes, and descriptions come from the assembly feature table, which covers 37,929 of them and agrees with the old LOCID table on all 33,210 they share. ::: ``` **R3 (02, fit the model)** ```text 2. Asserts that coefficient 2 is `phenotype_territorial_vs_social`. `social` is set as the reference level explicitly in `load_counts()`, so this holds by construction rather than by alphabetical accident (finding R3). 3. Keeps genes with 10 or more counts in at least a third of samples. This filtering happens *after* `DESeq()`, so dispersions were estimated on the full matrix — preserved from the original, which the committed results depend on. ``` **R7 (02, exports)** ```text Exports are written to full paths. The original wrapped these in `setwd("DESEQ_output/")` followed by two `setwd('..')` calls, which knitr resets at every chunk boundary (finding R7). ``` **B1-checklist (02, contrasts: alpha 0.05 vs 0.1 history)** ```text genes and so produces different adjusted p-values. Before this was fixed the unshrunken table (alpha 0.05) and the three shrunken tables (alpha 0.1) called different numbers of genes significant at the same 0.05 threshold. ``` **D2 (02, MA plots)** ```text ## PNG only. This was a PDF of 2-3 MB per tissue holding tens of thousands of ## vector points; gonad also carried a 2022 MA_plots.png that no code wrote. ``` **R2 (02, results tables)** ```text Significance is `padj < 0.05`, and only that. The previous code also defined `lfc.cutoff <- 1.5` and `svalue.cutoff <- 0.005`, never used them, and then labelled its outputs `DEG_apeglm-s0.005_lfc1.5` — describing a filter that never ran (finding R2). The unused variables are gone and the labels now say what was actually done. The counts themselves are unchanged. ``` **D2 (02, single-gene counts)** ```text The four estimators flag exactly the same genes, because they share one set of adjusted p-values. That is asserted here, and it is why one table of normalised counts serves all four: the counts depend on the gene, not on the estimator. Until 2026-09-27 there were four copies of that table, identical but for row order (submission checklist D2). ``` **D1 (02, summary rows)** ```text have `padj < 0.05` (`significant_*`). Until 2026-09-27 the first set was labelled `DEGs_all-genes*`, which read as a count of differentially expressed genes; it never was (submission checklist D1). ``` **D2 (02, volcano plots)** ```text ## Unshrunken and apeglm, shown here; only apeglm, the estimate the manuscript ## uses, is saved. Eight volcano PNGs per tissue were written until 2026-09-27, ## one per estimator for all and for significant genes (submission checklist ## D2); the manuscript's volcano plots are drawn by notebook 05. ``` **R4 (02, expression heatmap)** ```text Column annotation colours are derived from `coldata`. The previous code wrote `rep("royalblue1", 15)` / `rep("red3", 15)` for liver and `14` / `14` for gonad, hardcoding both the group sizes and the column order (finding R4). ``` **B2 (03, introduction; rule reconstructed from the committed outputs)** ```text Replaces `3_generate_gene_tables-salmon.Rmd`, which loaded `liver` and `gonad` and then operated on `control_ploidy`, `heat_ploidy`, `diploid_heat` and five other objects that never existed — leftovers from an oyster ploidy × multistressor experiment. It also read the DESeq2 CSVs with `sep=" "` when they were written with `sep=","`, and its final chunk did not parse (finding B2). Whatever produced the committed `gene_tables/` was therefore not in the repository. The rule has been reconstructed from those outputs and verified against both tissues: 1. Drop genes with **no description** — no match in the LOCID table. These are the rows the old oyster code called `LOCNA`. 2. Of the rest, **uncharacterized** are those whose description begins with `uncharacterized`; everything else is **characterized**. ``` **phase 3 (03, callout 'Output format changed deliberately': unquoted CSV in the 2023 gene tables)** ```text ::: {.callout-important} ## Output format changed deliberately The committed tables were written with `quote = FALSE` and `sep = ","`, but descriptions contain commas — `uncharacterized LOC115129809, transcript variant X3`. In `gonad_characterized.txt` only 754 of 1,280 rows have the expected 8 fields; 526 have 9, 10 or 11 and cannot be read back correctly. These are now written as properly quoted CSV. The rows and values are unchanged; only the quoting is. ::: ``` **B3, E13 (04, introduction; DAVID output and the other experiment's notebooks)** ```text Replaces the DAVID web-tool step that produced `DAVID_GOterms.txt` and `DAVID_KEGG_pathway.txt`, for which no script, input list, background or query date was recorded (both now in `archive/david-2023/`, see the README there). It also replaces `X_GOterm_Annotation.Rmd` and `X_GOSeq_annotation.Rmd`, which belonged to a different experiment entirely (finding B3). ``` **phase 3 (04, 'KEGG is available after all': the earlier misread organism list)** ```text **KEGG is available after all.** An earlier check of the KEGG REST organism list came back empty and was misread as "no such organism". KEGG does carry *Oncorhynchus nerka* as organism `one` (genome T07952, 37,942 proteins — the same gene count as our matrix). Its gene identifiers are NCBI GeneIDs, so they join to our `LOC` ids directly. The two REST responses are vendored alongside the GAF; see `tag-seq/genome/kegg_one_SOURCE.txt` for the URLs and retrieval date. No ortholog transfer is needed. ``` **PR #4 (04, gene lengths intro)** ```text This analysis uses the plain **hypergeometric** test instead, for three reasons that came out of the first clean-machine render (PR #4, 2026-09-06): ``` **PR #4, B1-checklist (04, reasons 2 and 3: the Wallenius numbers on the earlier gene set and platforms)** ```text against 5–13% in every other decile (Wilcoxon p ≈ 3e-119). Where the Wallenius weighting could be fitted it therefore down-weighted short genes heavily, and the pathways it suppressed were the ones made of short genes, oxidative phosphorylation first among them (on the gene set of PR #4, before the independent-filtering fix, on Linux: Wallenius p 3e-4, and no KEGG pathway at FDR < 0.05; hypergeometric, 59 of 186 genes DE, p 1e-17, on the current set). 3. **It was not reproducible.** The spline fitted by `nullp()` is numerically unstable against a decreasing relationship. On the PR #4 gene set it was platform-sensitive, and the Wallenius results sat at the threshold: the same inputs gave FDR 0.004 / 0.039 / 0.045 for the three enriched terms on macOS and 0.060 / 0.057 / 0.081 on Linux, so `check-reproduction.R` failed on the gonad enrichment tables and nothing else. On the current gene set, which differs from that one by a few hundred low-count genes, the fit fails outright: every fitted probability is the same constant, so a Wallenius test would silently return the hypergeometric result. The hypergeometric test is exact and deterministic. ``` **E7 (04, gene lengths source: the zero-byte GFF)** ```text The plan was to take lengths from `sequences/GCF_006149115.2_Oner_1.1_mRNA.gff`, but that file was tracked at **zero bytes** — the `bedtools` step that should have written it (documented in `01-upstream-alignment.qmd`, "Generate mRNA feature track") never produced anything that was committed. It has since been untracked rather than left in place as an empty file asserting that it exists (finding E7). Lengths come instead from the assembly feature table, which is committed and complete: the median mRNA interval length per gene, keyed by `GeneID`, which is our `LOC` identifier without the prefix. ``` **E11, C2 (05, introduction: figure_X)** ```text This notebook draws the three manuscript figures. It replaces `figures/figure_X`, which was assembled by hand in PowerPoint and so could not follow the data when the results changed. ``` **E11 (data README, introduction)** ```text What each tracked file is, and which notebook reads it. Written to close finding E11: 39 MB of this repository was tracked input that no current notebook reads, with nothing to say whether it was still needed. ``` **D3 (data README: 'tracked here until 2026-09-27')** ```text ## Moved out of the tree: the three provenance-only files Three files that no notebook reads were tracked here until 2026-09-27 (submission checklist D3). They are the only record of what the upstream pipeline produced besides the gene-level matrices above, so they were not deleted: they were removed from the working tree and remain in the repository's history, retrievable byte-for-byte from these commit-pinned addresses. ``` **E11, R8, D3 (genome README: section 'Removed: Onerka_LOCID_gene_table.txt')** ```text ## Removed: Onerka_LOCID_gene_table.txt The pre-2026 annotation table, removed from the tree on 2026-09-27 (submission checklist D3). No notebook read it. It is keyed entirely by LOC identifiers, so it can never annotate a symbol-named gene, which is what produced the fabricated `LOCNA.*` identifiers of finding R8. The feature table above replaces it and is a strict superset, re-verified before removal: all 33,211 genes in the LOC table are in the feature table, and every one of their descriptions matches one of that gene's mRNA names there. The only differences are CSV quoting and which transcript variant's name was taken (", transcript variant X1" against "X2"). Two copies remain: | Copy | sha256 | |---|---| | [commit-pinned, as tracked here](https://raw.githubusercontent.com/RobertsLab/project-sockeye-tagseq/97e1a1eafdef1e241aad325a3f4626cf30bf4f83/tag-seq/genome/Onerka_LOCID_gene_table.txt) | `931450449975391ce6135c0e06ed5d146ed783447769c2e35dae17fe21f5dd32` | | [Gannet](https://gannet.fish.washington.edu/panopea/berdahl-sockeye-salmon/genome/Onerka_LOCID_gene_table.txt) | `e60dbd4b9dbc3c0e06bfad3e71fe416ac2020ea9acb44236f39f8458953ffea5` | The Gannet copy has Windows (CRLF) line endings and is otherwise identical; that is the only reason the checksums differ. ``` **E11 (genome README, introduction)** ```text Written to close finding E11. ``` **E7 (sequences README: section 'Removed: GCF_006149115.2_Oner_1.1_mRNA.gff')** ```text ## Removed: GCF_006149115.2_Oner_1.1_mRNA.gff This file was tracked at **zero bytes** (finding E7). It is a derived intermediate: `01-upstream-alignment.qmd` ("Generate mRNA feature track") writes it by grepping mRNA features out of the genomic GFF and piping them through `bedtools`, and StringTie then consumes it as `-G`. An empty tracked file asserts that the artefact exists when it does not, so it has been untracked rather than left in place. To regenerate it, run the "Generate mRNA feature track" chunk of `01-upstream-alignment.qmd` on Raven. Nothing downstream needs it: goseq's gene lengths come from `genome/GCF_006149115.2_Oner_1.1_feature_table.txt`, which is committed and complete. ``` **E11, C2 (figures README: section 'History')** ```text ## History This directory used to hold `figure_X.png`, exported from a hand-assembled PowerPoint file (`figure_X.pptx`) that no script produced (finding E11). After the independent-filtering fix of 2026-09-26 its heatmap panels showed gene sets that no longer matched the results. Both files were removed when notebook 05 replaced them; they remain in the git history. ``` **R5, E9 (index.qmd exclusion callout, before update)** ```text ::: {.callout-warning} ## The two excluded gonad samples are not yet reproducibly justified `C05` and `C17` are dropped by a hardcoded `exclude` in `tag-seq/code/_common.R`, on the basis of a MultiQC sample-correlation inspection whose output is **not in this repository** — `tag-seq/QC/multiqc_gonad/` contains only a placeholder. This is the one analytical choice with a discretionary effect on the gonad results, and it currently rests on a statement rather than on evidence (audit findings R5 and E9). Treat the gonad counts as provisional until either the MultiQC output is committed or the exclusion criterion is computed inside notebook 02. ::: ``` **E10, D3 (index.qmd data availability, before update)** ```text `20220917-tagseq-liver/` (30 libraries), with FastQC and MultiQC reports alongside. This page previously linked a `seashell/bu-github/` address, which no longer resolves. They have **no archival accession yet** — no BioProject, no SRA (audit finding E10). ``` --- ## Licence and citation metadata (2026-09-27, checklist D7; closes E12) `LICENSE` and `CITATION.cff` exist, the last items open under E12. The licence is MIT for code and documentation and CC0 for data and results; the manuscript is not licensed; third-party files keep their own terms. The KEGG REST files are explicitly excluded from any licence, because KEGG provides its API for academic use only; whether they may stay in the repository is an open decision (submission checklist A4). --- ## KEGG downloaded at render time (2026-09-27, checklist A4) KEGG states it is not a public database, provides its REST API for academic use by academic users, and requires a licence for other use; Bioconductor's KEGG.db manual says the package "is no longer allowed to be current" for the same reason. The two vendored KEGG responses are therefore removed from the repository (and so from any Zenodo archive; they remain in git history, which was not rewritten). - `ensure_kegg()` in `_common.R` downloads them into `tag-seq/genome/` on first use (git ignores them), records the date and KEGG release in `kegg_one_RETRIEVED.txt`, and reuses them after that; `KEGG_REFRESH=1` forces a new download. Notebook 04 and the manuscript call it, and 04 reports whether the responses match those recorded for the committed results. - `kegg_one_SOURCE.txt` records the release (2026/09/25) and the sha256 of the two responses the committed results were computed from. They left `CHECKSUMS.sha256`, where a new KEGG release would have looked like a corrupted input. - `check-reproduction.R` accepts numeric drift in KEGG-derived tables, listed in their own section, only when the set of pathways enriched at FDR < 0.05 is unchanged. Rebaseline: the live KEGG links differed from the 2026-08-30 copy by one removed and five added gene-pathway links, none involving a significant gene. Rendering against them left all five enriched gonad pathways and their gene counts unchanged, moved p-values in the third significant figure, and changed one tested pathway's gene count by 2. The manuscript's KEGG counts in the methods are now computed from the downloaded data (12,284 genes in 196 pathways; the hand-typed figure was 12,285). --- ## How you will know it worked 1. **It renders.** `quarto render` completes from a clean checkout with no manual steps. This is the bar the repository failed at audit time. 2. **It renders the same.** Diff phase-2 outputs against `v0-preexisting`. Every difference should be one you can name and defend. 3. **It renders elsewhere.** The GitHub Action is the real test: a machine that has never seen this project, restoring from `renv.lock`, producing the same numbers.