Decision log
Decision log
Dated, numbered, append-only. Each entry: the problem, the choice, why.
D-001 (2026-10-06) Species codes and sample ids
Upstream uses Ptua/ptua in file names and Ptuh in counts.qmd. This repo uses Ptuh everywhere (it is the species, P. tuahiniensis; the genome is P. meandrina). Canonical sample id is {ACR|POR|POC}-{colony}-TP{n}; norm_sample_id() in R/harmonize.R accepts every upstream variant seen (plate prefixes, underscores, BAM paths, readr dedup suffixes, timepoint1).
D-002 (2026-10-06) Technical replicates in the Apul gene matrix
MOSAiC’s docs/corrected_sample_counts.md reports 80 columns (technical replicates) in the remote Apul gene matrix. On 2026-10-06 the remote file on gannet is byte-identical (MD5 6d9c7c720951658ab7c4a552c9bfc926) to the 40-column copy in MOSAiC/data/orthologs/. The replicate issue is therefore resolved upstream and no collapsing is done here. harmonize_matrix() stops on duplicate sample ids, so a reappearance will fail loudly rather than silently averaging.
D-003 (2026-10-06) Gene id normalization
Apul ortholog ids carry a -T\d+ transcript suffix; Peve and Ptuh count matrices carry a gene- prefix. Both are stripped (norm_gene_id()). After stripping, an Apul gene can appear in more than one ortholog group if two of its transcripts were assigned separately; genes_by_ortholog() keeps the first group per gene. Revisit if a hypothesis depends on isoform resolution.
D-004 (2026-10-06) Where compute runs
GitHub Actions runs only the fetch, harmonize, reproduction gate, and report render. Hypothesis analyses that need Bioconductor (DESeq2, variancePartition, MOFA2) or the 25 MB CpG matrices run locally or on raven via targets::tar_make(); their small outputs are committed. Revisit if Actions runtime proves sufficient.
D-005 (2026-10-06) What is committed under data/
data/raw/ is never committed. data/derived/ is committed except the CpG matrices (ignored by pattern). config/design.csv is committed and is the single sample-metadata table; colony-level site, nutrient, and haplotype are taken from the physiology master. Samples present in molecular layers but absent from the physiology master get NA for those fields rather than being dropped.
D-006 (2026-10-06) Timepoint dates
Exact collection dates per colony are not in any upstream table. The design uses mid-month nominal dates (15 Jan, 15 Mar, 15 Sep, 15 Nov 2020), as temperature.qmd in MOSAiC does. Hypotheses that compute thermal history (H11) must treat the window as approximate.
D-007 (2026-10-06) lncRNA sample counts differ from MOSAiC’s documentation
MOSAiC docs/corrected_sample_counts.md lists 45 / 43 / 44 lncRNA samples (Apul / Peve / Ptuh). The lncRNA files pinned in the lock on 2026-10-06 have 40 / 38 / 39 sample columns, identical to the gene matrices, plus the five featureCounts metadata columns (Chr, Start, End, Strand, Length). The documented numbers evidently came from an earlier upstream version. The reproduction gate asserts the current values; if upstream changes again the weekly drift check will catch it.
D-008 (2026-10-06) CpG matrix sizes are very unequal across species
Sites passing the upstream filter (>= 10x in every sample): Apul 96,772; Peve 242,100; Ptuh 1,991,606. The Ptuh matrix is 397 MB raw. Part of this is fewer Ptuh samples (32) making the all-samples filter easier to pass, but a 20-fold difference suggests the upstream filtering may not be identical across species. Methylation hypotheses (H04 to H06) must report per-species site counts and should consider a common coverage filter applied here.
D-009 (2026-10-07) lme4 pinned to 1.1-37
variancePartition 1.32.5 (Bioconductor 3.18, R 4.3) imports lme4::findbars, which lme4 2.0 moved to the reformulas package; with lme4 2.0.6 every fitExtractVarPartModel call fails. lme4 1.1-37 is installed from the CRAN archive, from source so it matches the local Matrix ABI. Do not run update.packages() on a machine that runs H01 or H04 without re-pinning: install.packages("https://cloud.r-project.org/src/contrib/Archive/lme4/lme4_1.1-37.tar.gz", repos = NULL, type = "source"). Revisit when moving to a Bioconductor release whose variancePartition supports lme4 2.
D-010 (2026-10-07) H04 runs on raven; varpart is chunked
A desktop run of H04 (16 GB, 7 forked workers) swapped once it reached the Ptuh CpGs: each worker grew to about 2 GB, 10.6 GB of swap was in use, and the run was stopped after 8.5 h without finishing. fit_varpart() now fits features in blocks of 100k with gc() between blocks (results verified identical to a single call), and H04 reads H04_WORKERS to cap parallelism. H04 runs on raven, as the README compute tier says. Ptuh accounts for about 85% of the CpG fits (see D-008).
D-011 (2026-10-07) miRNA libraries below 1,000 reads are dropped
The upstream miRNA count matrices contain near-empty libraries. Peve, eight below 1,000 total reads against a median of about 20,000: POR-236-TP2 (0), POR-74-TP4 (3), POR-262-TP3 (5), POR-262-TP4 (5), POR-69-TP2 (55), POR-83-TP3 (69), POR-262-TP2 (241), POR-216-TP4 (351); the next lowest is 1,368. Ptuh, one: POC-52-TP1 (401, median about 72,000). DESeq2 cannot size-factor an all-zero sample and the others are noise. Analyses using miRNA call mirna_qc_samples() (R/expression.R), which keeps samples with >= 1,000 total reads. The design table’s has_mirna flags are unchanged, so the reproduction gate still matches MOSAiC’s published sample counts. Worth an upstream issue on urol-e5/timeseries_molecular.
Peve miRNA depth is also confounded with timepoint among the libraries that pass the 1,000-read filter (29 of 37): median total reads TP1 182k, TP2 82k, TP3 52k, TP4 11k (n = 10 / 7 / 6 / 6; median poscounts size factors 3.5, 0.56, 0.57, 0.31). All eight failed libraries are TP2 to TP4. Peve miRNA timepoint effects (H07, H13) should be read as possibly technical until this is checked upstream.
D-012 (2026-10-07) R environment on raven
Raven’s system R is 4.3.3 with no project packages and no 4.3 user library. Packages were installed into ~/R/x86_64-pc-linux-gnu-library/4.3 (R adds it automatically once it exists) from Bioconductor 3.18, matching the desktop, with lme4 pinned to 1.1-37 (D-009). Two workarounds were needed, neither a system change:
fs: raven has no libuv headers. Install withUSE_BUNDLED_LIBUV=1.Deriv(needed by doBy -> pbkrtest -> lmerTest -> variancePartition): the current CRAN release requires R >= 4.5. Install the archived 4.2.0:install.packages("https://cloud.r-project.org/src/contrib/Archive/Deriv/Deriv_4.2.0.tar.gz", repos = NULL, type = "source").
A fresh fetch on raven reproduced every md5 in config/upstream.lock.yml; only fetched_on dates changed, and that diff was discarded. The reproduction gate passed before H04 ran.
D-013 (2026-10-07) H04 PERMANOVA bootstrap runs in parallel
In the first raven run the R2 bootstrap (200 resamples, one Euclidean distance over up to 1.8 M CpGs each) ran serially and took about 1 h 45 min of the 3 h 44 min total. permanova_one() now draws all bootstrap indices serially in the main process, in the same order replicate() did, and forks only the adonis2 fits with parallel::mclapply (H04_WORKERS cores). adonis2(..., permutations = 0) uses no random numbers, so the draws and the RNG state afterwards are unchanged: rerunning the purled qmd for Apul reproduced the committed permanova.csv rows exactly (tolerance 0). Each Ptuh worker holds one resampled copy of the matrix, about 1 GB.
Open: variancePartition fits ran at about 0.24 s per feature per worker with 47 workers, against 0.026 s in a single-process benchmark. Check whether threaded BLAS inside forked workers is the cause before the next large run.
D-014 (2026-10-07) lncRNA features have no strand
lncrna_<species>_features comes from the featureCounts table header (Chr, Start, End, Strand, Length) of the upstream *_lncRNA_counts.clean.filtered.txt files, and Strand is + for every lncRNA in all three species (15,559 / 8,319 / 11,236). Gene coordinates from the GFFs are about half -, so this is a property of the lncRNA SAF/GTF used upstream, not of harmonization. Treat lncRNA strand as unknown. H17 drops its orientation split and excludes all pairs within 2 kb in its read-through check (H17 Amendment 1). Asked upstream whether stranded lncRNA coordinates exist: https://github.com/urol-e5/timeseries_molecular/issues/139.
D-015 (2026-10-08) CpG matrices carry no coverage
cpg_<species> comes from merged-WGBS-CpG-counts_filtered.csv, which holds percent methylation for CpGs with >= 10x in every sample, and no read counts. Anything that needs per-sample coverage (binomial noise models, depth filters above 10x, beta-binomial tests) cannot be done from the manifest as it stands. H19 uses the 10x floor as a fixed noise bound (H19 Amendment 1). If coverage becomes necessary, the per-sample Bismark coverage files would have to be added to config/upstream.yml; they are large and should be read on raven.
D-016 (2026-10-08) proposed status for incomplete pre-registrations
H23–H27 are listed on the status board before their pre-registrations are complete (open items include matching bins, bootstrap counts, and, for H27, an effect threshold and power simulation). planned means the decision rule is frozen, so these use a new status, proposed, which comes before planned. A proposed hypothesis has no analysis.qmd and no target in _targets.R. Moving it to planned is a commit that resolves its “Before moving to planned” list; from then on CLAUDE.md section 2 applies unchanged.
D-017 (2026-10-09) Peve reads are dominated by one scaffold_167 locus
In Peve, a median 69% of raw gene + lncRNA reads per library fall on one locus, Porites_evermani_scaffold_167 about 401.5–408.4 kb: five nested lncRNAs (lncRNA_4889–lncRNA_4893, the three largest each 15.8% of all reads, i.e. the same reads counted more than once) and four short adjacent genes (Peve_00009616–Peve_00009619, 200–300 bp, no ortholog or SwissProt annotation). The share rises across the year (median TP1 0.57, TP2 0.72, TP3 0.70, TP4 0.78). The top-10-feature read fraction is 0.70 in Peve against 0.07–0.08 in Apul and Ptuh; with the locus removed it is 0.13. This looks like a highly expressed repetitive or rRNA-like element annotated as both genes and lncRNAs. Found while completing H26’s pre-registration (hypotheses/H26/feasibility_check.R; per-library totals only).
H26 excludes these nine features in Peve and computes library-composition covariates without them. Earlier Peve analyses used the matrices as delivered; DESeq2’s median-of-ratios size factors are robust to a few dominant features, but a time-varying share of this size could still leak into Peve timepoint effects through normalization. Not yet raised upstream (urol-e5/timeseries_molecular); decided 2026-10-09 to record it here first.
D-018 (2026-10-10) Threaded BLAS on raven: worker count matters more
Raven R links multithreaded OpenBLAS (48 threads by default) on 24 physical cores (Xeon Gold 5220R, 2 threads per core). Forked analyses use detectCores() - 1 = 47 workers, and the H08 rerun reached load average ~136. R/parallel.R adds pin_blas_threads(n = 1): it sets OPENBLAS_NUM_THREADS/OMP_NUM_THREADS/MKL_NUM_THREADS (for child R processes) and calls RhpcBLASctl::blas_set_num_threads() (for this process and its forks), and skips the latter if RhpcBLASctl is missing.
Benchmarks on raven, 2026-10-10, machine otherwise idle (load < 10), seconds of wall time, default = 48 BLAS threads, pinned = 1:
| piece | workers | default | pinned |
|---|---|---|---|
H08 count_edges, Apul 15,559 x 21,108, one pass, no fork |
1 | 16.7 | 9.4 |
H08 count_edges, 47 or 48 permutations |
47 | 141, 143, 145 | 160, 161, 161 |
| same | 24 | 77 | 74 |
| same | 16 | 50 | |
| same | 12 | 90 | 48 |
| same | 6 | 76 | |
fit_varpart(), Apul 20,000 genes x 40 samples |
47 | 31, 27, 28 | 30, 28, 28 |
| same | 24 | 30 | 30 |
| same | 12 | 46 | 46 |
| same, 500 genes serial | 1 | 11.1 | 11.1 |
H18 full render (quarto render) |
47 | 19 min (committed run) | 16 min |
Findings:
- The variancePartition / lme4 fits never start BLAS threads (2 threads per worker in both modes; 40 x few-column matrices are below OpenBLAS’s threading cutoff). Pinning does not change their speed, so it is not the cause of H04’s 0.24 s/feature/worker (the open item in D-013). Here the 47-worker rate was 0.065 s/feature/worker against 0.022 serial, and throughput stops improving at about 24 workers, the physical core count. H04’s slower rate on the 1.8 M Ptuh CpGs is more likely memory per fork.
- H08’s correlation pass is memory-bandwidth bound. 47 workers are about 3x slower than 12–16 pinned workers for the same 48 permutations. Pinning helps only once workers are reduced (12 workers: 90 s default, 48 s pinned); at 47 workers it is about 12% slower.
- Pinning changes lme4 results beyond floating-point level. Default vs pinned
fit_varpart()on 20,000 genes: 5,738 features differ by > 1e-6, 7 by > 1e-3, max 0.024 (fits on the zero-variance boundary); medians of every component agree to 1e-5. Each mode is bit-reproducible on rerun. H18 rendered pinned matched every committed output CSV to <= 6e-10 absolute, with identical text columns and verdicts.
Choice: pin_blas_threads() is opt-in, and neither _targets.R nor a project .Rprofile/.Renviron calls it. A global pin on its own makes no analysis faster at the current worker counts (H08 gets slower). The gain needs fewer workers, and most analyses hard-code detectCores() - 1, so that change means editing their analysis.qmd and marking the targets outdated. A global pin would also make H01/H04 reruns drift from their committed varpart outputs with no speed benefit. For a raven rerun, use the existing overrides:
OPENBLAS_NUM_THREADS=1 H08_WORKERS=12 Rscript -e 'targets::tar_make(H08)'
H04_WORKERS=24 Rscript -e 'targets::tar_make(H04)'(Environment variables set in the calling shell reach the knitr R process; checked in /proc/<pid>/environ during the H18 render.) When an analysis that forks is next edited for other reasons, cap its workers (about 16 for memory-bound passes, 24 for CPU-bound fits) and call pin_blas_threads() before forking. The change does not touch decision rules. H09’s WGCNA threads were not benchmarked.