Sockeye salmon TagSeq: territorial vs social
Gene expression in gonad and liver, Oncorhynchus nerka
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 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): 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.
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 | FASTQs on Gannet | count matrices |
| 02 — Differential expression (liver) | tag-seq/data/ |
tag-seq/DESEQ_output/liver/ |
| 02 — Differential expression (gonad) | tag-seq/data/ |
tag-seq/DESEQ_output/gonad/ |
| 03 — Gene tables | significant DEGs | tag-seq/gene_tables/ |
| 04 — GO and KEGG enrichment | all DEGs, GAF, KEGG | tag-seq/GO_output/ |
| 05 — Manuscript figures | counts, result tables | figures/ |
| 06 — Supplementary tables | 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/— per-tissue DESeq2 tables (unshrunken andapeglm), the summary counts, PCA, sample-correlation and expression heatmaps, a volcano plot and MA plots. Thenormalandashrestimates are regenerated on every render but not tracked; their fold changes are in Supplementary Tables S2 and S3.tag-seq/gene_tables/— significant genes split by annotation status, and the LOC identifier liststag-seq/GO_output/— full and FDR-filtered GO and KEGG results
Data availability
Raw FASTQ files are on Gannet, at https://gannet.fish.washington.edu/panopea/berdahl-sockeye-salmon/: 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.