--- title: "Sockeye salmon TagSeq: territorial vs social" subtitle: "Gene expression in gonad and liver, *Oncorhynchus nerka*" date: today --- This is the rendered analysis for the Berdahl sockeye salmon TagSeq experiment. Every number and figure below comes from re-running the notebooks in this repository; nothing is transcribed by hand. Reproduce it with `renv::restore()` followed by `./render-all.sh` from the project root. See the [repository README](https://github.com/RobertsLab/project-sockeye-tagseq#reproducibility) for the entry point and the expected-difference statement, and `REPRODUCIBILITY-PLAN.md` for the audit trail behind the current state of the code. ## Design TagSeq (3'-end mRNA sequencing) of two tissues from fish assigned to two behavioural treatments, territorial and social. | Tissue | Samples sequenced | Territorial | Social | Analysed | Excluded | |---|---|---|---|---|---| | Liver | 30 (`B01`–`B30`) | 15 | 15 | 30 | — | | Gonad | 30 (`C01`–`C30`) | 15 | 15 | 28 | `C05`, `C17` | Counts are gene-level, quantified against NCBI assembly [GCF_006149115.2 (Oner_1.1)](https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_006149115.2/): 37,942 genes before filtering. Treatment contrast is territorial vs social with `social` as the reference level, tested with DESeq2 and reported at `padj < 0.05`. ::: {.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`. The recorded reason is a sample-correlation inspection whose output is **not in this repository**; the gonad MultiQC reports, now in `tag-seq/QC/`, cover FastQC only and have no correlation heatmap. What the repository can show is that they are the two smallest gonad libraries by gene-assigned counts (1.97 and 2.14 million, against 3.39–5.77 million for the other 28; Supplementary Table S1). This is the one analytical choice with a discretionary effect on the gonad results. Treat the gonad counts as provisional until the exclusion is either supported by the original evidence or stated as a rule computed in notebook 02 (submission checklist A5). ::: ## Notebooks | | Reads | Writes | |---|---|---| | [01 — Upstream alignment](tag-seq/code/01-upstream-alignment.html) | FASTQs on Gannet | count matrices | | [02 — Differential expression (liver)](tag-seq/code/02-differential-expression-liver.html) | `tag-seq/data/` | `tag-seq/DESEQ_output/liver/` | | [02 — Differential expression (gonad)](tag-seq/code/02-differential-expression-gonad.html) | `tag-seq/data/` | `tag-seq/DESEQ_output/gonad/` | | [03 — Gene tables](tag-seq/code/03-gene-tables.html) | significant DEGs | `tag-seq/gene_tables/` | | [04 — GO and KEGG enrichment](tag-seq/code/04-enrichment.html) | all DEGs, GAF, KEGG | `tag-seq/GO_output/` | | [05 — Manuscript figures](tag-seq/code/05-figures.html) | counts, result tables | `figures/` | | [06 — Supplementary tables](tag-seq/code/06-supplementary-tables.html) | counts, result tables | `manuscript/supplementary/` | Notebook 01 documents the FASTQ-to-count-matrix pipeline as it was run on Raven (HISAT2, StringTie, `prepDE.py`). Every chunk in it is `eval: false` — it is a record, not a step this render performs. ## Results as committed Differential expression, territorial vs social, `padj < 0.05`, `apeglm` shrinkage. Adjusted p-values come from the Wald test with independent filtering at `alpha = 0.05`, and the shrunken tables carry those same values, so every table makes the same significance calls: | Tissue | Genes tested | Significant | Characterized | Uncharacterized | |---|---|---|---|---| | Liver | 4,076 | 66 | 66 | 0 | | Gonad | 18,156 | 1,653 | 1,609 | 42 | Rows with no annotation available: 11 of 4,076 (liver), 12 of 18,156 (gonad), including 2 of the significant gonad genes. "Uncharacterized" means NCBI describes the gene as `uncharacterized LOC…`; 2.5% of the significant gonad genes. Enrichment, FDR < 0.05: | Tissue | GO terms tested | GO enriched | KEGG pathways tested | KEGG enriched | |---|---|---|---|---| | Liver | 3,109 | 0 | 188 | 0 | | Gonad | 6,011 | 36 | 195 | 5 | Both use `goseq`'s hypergeometric over-representation test, GO against the NCBI GO annotation for this assembly (155,362 annotations, 32,629 genes, 6,937 terms) and KEGG against organism `one`, Benjamini–Hochberg corrected. The gonad signal is ribosomal and mitochondrial: the top KEGG pathways are oxidative phosphorylation (`one00190`, 59 of 186 genes, FDR 2e-15) and ribosome (`one03010`, 66 of 261, FDR 3e-12), and the top GO terms are structural constituent of ribosome, translation and the mitochondrial respiratory chain. `goseq`'s default Wallenius test, which re-weights genes by transcript length, is not used. TagSeq is 3′-end counting, so read count does not scale with length; the length association in these data runs the other way, with short genes far more likely to be differentially expressed; and the Wallenius result turned out to be numerically unstable: on an earlier gene set the three terms it called at FDR 0.004–0.045 on macOS fell to 0.057–0.081 on Linux, and on the current gene set its length-bias spline cannot be fitted and collapses to a constant. The reasoning is set out in notebook 04. That the liver comparison yields no enriched terms is a result, not a failure: 66 significant genes of small effect is little power for an over-representation test. The closest liver pathway is oxidative phosphorylation, at FDR 0.11. ## Result tables Tracked in git, so a re-render either reproduces them or shows a diff: - [`tag-seq/DESEQ_output/`](https://github.com/RobertsLab/project-sockeye-tagseq/tree/main/tag-seq/DESEQ_output) — per-tissue DESeq2 tables (unshrunken and `apeglm`), the summary counts, PCA, sample-correlation and expression heatmaps, a volcano plot and MA plots. The `normal` and `ashr` estimates are regenerated on every render but not tracked; their fold changes are in Supplementary Tables S2 and S3. - [`tag-seq/gene_tables/`](https://github.com/RobertsLab/project-sockeye-tagseq/tree/main/tag-seq/gene_tables) — significant genes split by annotation status, and the LOC identifier lists - [`tag-seq/GO_output/`](https://github.com/RobertsLab/project-sockeye-tagseq/tree/main/tag-seq/GO_output) — full and FDR-filtered GO and KEGG results ## Data availability Raw FASTQ files are on Gannet, at : gonad in `20220714-tagseq-gonad/` (30 libraries, two lanes each), liver in `20220917-tagseq-liver/` (30 libraries), with FastQC and MultiQC reports alongside. They have **no archival accession yet**: no BioProject and no SRA deposit (submission checklist A3). Until they are deposited, the raw data has neither a persistence guarantee nor a citable identifier, and the analysis cannot be re-derived from source by anyone outside the lab. The committed count matrices are checksummed in `CHECKSUMS.sha256` and verified at the start of every render, which pins what this analysis was run on but does not substitute for the deposit.