[study] CIBNOR2021_HEAT_RRSS (PRJNA690951): first reanalysis

config/CIBNOR2021_HEAT_RRSS.config

@@ -22,11 +22,11 @@
  *     --attribute family --attribute replicate --attribute tissue --attribute temp \
  *     --methylkit-levels SS RR
  *
- * NOT YET RUNNABLE AS A VALID ANALYSIS. DMR_METHYLKIT fits RR vs SS with no
- * family term, so it would treat 6 vs 6 libraries as independent when the
- * phenotype is assigned per family (2 vs 2). See the study YAML limitations.
- * Model family, or aggregate to family level, before computing results that
- * enter a pool.
+ * The registered result (2026-10-16) was produced with these settings:
+ * family as the replicate unit, relaxed Bismark scoring, read-2 M-bias trim,
+ * and limma on family-pooled, destranded 1 kb tiles (PR #67). The methylation
+ * work dir must be on a volume with no space in its name: Trim Galore and
+ * Bismark resolve and pass it unquoted.
  */
 
 includeConfig 'base.config'
@@ -52,6 +52,24 @@ params {
 
         // SS = 0 (reference), RR = 1, keyed by run accession.
         sample_sheet = 'data/studies/CIBNOR2021_HEAT_RRSS/methylkit_samplesheet.csv'
+
+        // Thermal tolerance is a family property: 2 families per class.
+        replicate_unit = 'family'
+        design_sheet   = 'data/studies/CIBNOR2021_HEAT_RRSS/samplesheet.tsv'
+
+        // 53.9% vs 15.5% unique mapping on a 50k-pair test
+        // (data/studies/CIBNOR2021_HEAT_RRSS/bismark_scoring_diagnostic.txt).
+        bismark_align_args = '--score_min L,0,-0.6'
+        // M-bias: read 2 positions 1-20 are end-repair biased; position 150
+        // of both reads is an artifact.
+        extractor_args     = '--ignore_r2 20 --ignore_3prime 1 --ignore_3prime_r2 1'
+
+        // limma on family-pooled tiles; coverage applied after pooling.
+        dmr_engine        = 'limma'
+        dmr_features      = 'genes'
+        min_coverage      = 1
+        unit_min_coverage = 10
+        tile_min_cpgs     = 10
     }
 
     outdir = 'results/methylation/CIBNOR2021_HEAT_RRSS'
@@ -64,11 +82,17 @@ params {
 
 process {
     // 12 HiSeq 2000 PE libraries, 54-73 M read pairs each (~30x of a ~0.6 Gb
-    // genome). Bismark runs several bowtie2 instances per --multicore.
-    withName: 'BISMARK_GENOME_PREPARATION'    { cpus = 8;  memory = '32 GB'; time = '6h'  }
-    withName: 'BISMARK_ALIGN'                 { cpus = 12; memory = '48 GB'; time = '24h' }
-    withName: 'BISMARK_DEDUPLICATE'           { cpus = 2;  memory = '16 GB'; time = '6h'  }
-    withName: 'BISMARK_METHYLATION_EXTRACTOR' { cpus = 8;  memory = '16 GB'; time = '12h' }
-    withName: 'TRIM_GALORE'                   { cpus = 4;  memory = '8 GB';  time = '4h'  }
-    withName: 'DMR_METHYLKIT'                 { cpus = 4;  memory = '32 GB'; time = '6h'  }
+    // genome). Bismark runs several bowtie2 instances per --multicore. These
+    // fit a 10-CPU, 26 GB workstation (Chlo), one alignment at a time, where
+    // each library took ~4-7 h at the relaxed score.
+    withName: 'FASTQC'                        { cpus = 2;  memory = '4 GB';  maxForks = 2 }
+    withName: 'TRIM_GALORE'                   { cpus = 4;  memory = '8 GB';  maxForks = 1 }
+    withName: 'BISMARK_GENOME_PREPARATION'    { cpus = 8;  memory = '16 GB'; time = '12h' }
+    // --multicore = cpus / 3 = 3 Bismark instances; one library at a time.
+    withName: 'BISMARK_ALIGN'                 { cpus = 9;  memory = '22 GB'; time = '48h'; maxForks = 1 }
+    withName: 'BISMARK_DEDUPLICATE'           { cpus = 2;  memory = '12 GB'; time = '12h'; maxForks = 2 }
+    withName: 'BISMARK_METHYLATION_EXTRACTOR' { cpus = 6;  memory = '12 GB'; time = '24h'; maxForks = 1 }
+    withName: 'COVERAGE_FILTER'               { cpus = 1;  memory = '2 GB' }
+    withName: 'DMR_METHYLKIT'                 { cpus = 4;  memory = '24 GB'; time = '12h' }
+    withName: 'DMR_POOLED_TILES'              { cpus = 2;  memory = '22 GB'; time = '6h' }
 }

data/studies/CIBNOR2021_HEAT_RRSS/CIBNOR2021_HEAT_RRSS_rr_vs_ss_day30_oscillating_heat_annotation_qc.tsv

@@ -0,0 +1,6 @@
+study_id	comparison_id	annotation_context	n_regions
+CIBNOR2021_HEAT_RRSS	rr_vs_ss_day30_oscillating_heat	promoter	0
+CIBNOR2021_HEAT_RRSS	rr_vs_ss_day30_oscillating_heat	exon	0
+CIBNOR2021_HEAT_RRSS	rr_vs_ss_day30_oscillating_heat	intron	0
+CIBNOR2021_HEAT_RRSS	rr_vs_ss_day30_oscillating_heat	gene_body	31166
+CIBNOR2021_HEAT_RRSS	rr_vs_ss_day30_oscillating_heat	intergenic	101421

data/studies/CIBNOR2021_HEAT_RRSS/CIBNOR2021_HEAT_RRSS_rr_vs_ss_day30_oscillating_heat_dmr_qc.tsv

@@ -0,0 +1,2 @@
+study_id	comparison_id	n_samples	n_candidate_regions	n_significant_regions	dmr_mode	tile_size	min_coverage	qvalue_cutoff	meth_diff_cutoff	replicate_unit	n_replicate_units	overdispersion	test	engine	df_prior	destranded	tile_min_cpgs	n_cpgs_all_units
+CIBNOR2021_HEAT_RRSS	rr_vs_ss_day30_oscillating_heat	12	132587	0	genes	1000	10	0.05	10	family	4	empirical-Bayes moderated variance	limma moderated t on unit M-values, BH	limma	2.33914431336831	TRUE	10	7269210

data/studies/CIBNOR2021_HEAT_RRSS/CIBNOR2021_HEAT_RRSS_rr_vs_ss_day30_oscillating_heat_workflow_manifest.json

@@ -0,0 +1,104 @@
+{
+  "workflow_name": "aree-methylation",
+  "workflow_version": "0.1.0-scaffold",
+  "study_id": "CIBNOR2021_HEAT_RRSS",
+  "comparison_id": "rr_vs_ss_day30_oscillating_heat",
+  "mode": "raw_reanalysis",
+  "date_generated": "2026-10-16T15:53:23.785171-07:00",
+  "generated_by": "automated:aree-methylation@0.1.0-scaffold",
+  "parameters": {
+    "mode": "raw_reanalysis",
+    "study_id": "CIBNOR2021_HEAT_RRSS",
+    "comparison_id": "rr_vs_ss_day30_oscillating_heat",
+    "min_coverage": 1,
+    "qvalue_cutoff": 0.05,
+    "meth_diff_cutoff": 10.0,
+    "dmr_mode": "tile",
+    "tile_size": 1000,
+    "tile_step": 1000,
+    "dmr_engine": "limma",
+    "dmr_features": "genes",
+    "unit_min_coverage": 10,
+    "tile_min_cpgs": 10,
+    "replicate_unit": "family",
+    "design_sheet": "data/studies/CIBNOR2021_HEAT_RRSS/samplesheet.tsv",
+    "bismark_align_args": "--score_min L,0,-0.6",
+    "extractor_args": "--ignore_r2 20 --ignore_3prime 1 --ignore_3prime_r2 1",
+    "promoter_upstream_bp": 1500,
+    "genome_assembly": "GCF_963853765.1"
+  },
+  "input_checksums": {
+    "CIBNOR2021_HEAT_RRSS_rr_vs_ss_day30_oscillating_heat_dmr_standardized.tsv": "sha256:3dd8f2e9c7911f71f2b709320e94cc487c23113c1317297e362e3a4969069bf3",
+    "CIBNOR2021_HEAT_RRSS_rr_vs_ss_day30_oscillating_heat_methylkit_qc.tsv": "sha256:021e34c1d28c75cb63c64c083e31d106e5057d91513accdb36575553ef6fc4c3",
+    "CIBNOR2021_HEAT_RRSS_rr_vs_ss_day30_oscillating_heat_annotation_qc.tsv": "sha256:9c9603f4e48acb7edd8e834abf21e3d8515b9a6891b245f81dd5f5f93f683829",
+    "CIBNOR2021_HEAT_RRSS_rr_vs_ss_day30_oscillating_heat_standardize_qc.tsv": "sha256:8b849b0d33fa8fe6e35cec6483378d5741bd3332ef67078d745c16459fa37ab8"
+  },
+  "software_versions": {
+    "ANNOTATE_REGIONS": {
+      "bioconductor-genomicranges": "1.62.1",
+      "bioconductor-rtracklayer": "1.70.1",
+      "r-base": "4.5.1"
+    },
+    "BISMARK_ALIGN": {
+      "bismark": "0.24.2"
+    },
+    "BISMARK_DEDUPLICATE": {
+      "bismark": ""
+    },
+    "BISMARK_GENOME_PREPARATION": {
+      "bismark": "0.24.2"
+    },
+    "BISMARK_METHYLATION_EXTRACTOR": {
+      "bismark": ""
+    },
+    "COVERAGE_FILTER": {
+      "awk": "awk version 20200816"
+    },
+    "DMR_POOLED_TILES": {
+      "bioconductor-bsseq": "1.46.0",
+      "bioconductor-dss": "2.58.0",
+      "bioconductor-limma": "3.66.0",
+      "r-base": "4.5.1",
+      "r-data.table": "1.17.8"
+    },
+    "FASTQC": {
+      "fastqc": "0.12.1"
+    },
+    "STANDARDIZE_OUTPUT": {
+      "python": "3.14.8"
+    },
+    "TRIM_GALORE": {
+      "trim_galore": "0.6.10"
+    },
+    "EMIT_MANIFEST": {
+      "python": "3.14.8"
+    }
+  },
+  "software_versions_not_captured": {
+    "RENDER_REPORT": "renders this manifest, so it runs afterwards; its versions.yml is published to methylation/report/"
+  },
+  "container_engine": null,
+  "qc_metrics": {
+    "n_standardized_regions": 132587,
+    "n_hyper": 62617,
+    "n_hypo": 69970,
+    "mean_coverage": null,
+    "pct_regions_passing_coverage_filter": null
+  },
+  "warnings": [
+    "software_version_unresolved: BISMARK_DEDUPLICATE wrote an empty version for bismark (version command failed or tool missing)",
+    "software_version_unresolved: BISMARK_METHYLATION_EXTRACTOR wrote an empty version for bismark (version command failed or tool missing)"
+  ],
+  "output_schema": [
+    "region_id",
+    "chrom",
+    "start",
+    "end",
+    "gene_id",
+    "annotation_context",
+    "meth_diff_percent",
+    "qvalue",
+    "direction"
+  ],
+  "harmonize_command": "aree harmonize --study CIBNOR2021_HEAT_RRSS --input CIBNOR2021_HEAT_RRSS_rr_vs_ss_day30_oscillating_heat_dmr_standardized.tsv"
+}

data/studies/CIBNOR2021_HEAT_RRSS/bismark_scoring_diagnostic.txt

@@ -0,0 +1,9 @@
+A_default            map=15.5%  CpG=16.5%  CHG=2.6%   CHH=3.0%  
+B_score_L0_-0.6      map=53.9%  CpG=14.4%  CHG=2.8%   CHH=2.9%  
+C_non_directional    map=15.5%  CpG=16.5%  CHG=2.6%   CHH=3.0%  
+D_maxins_1000        map=15.8%  CpG=16.5%  CHG=2.6%   CHH=3.0%  
+E_R1_single          map=37.2%  CpG=14.3%  CHG=0.8%   CHH=0.8%  
+F_R2_single_pbat     map=17.2%  CpG=17.7%  CHG=4.1%   CHH=5.0%  
+
+[exited with code 0]
+# 50k trimmed pairs of SRR13482032, Bismark 0.24.2 / bowtie2 2.5.5; script /Volumes/blue-block/aree/tmp/bismark_diag/diag.sh

data/studies/CIBNOR2021_HEAT_RRSS/library_qc.tsv

@@ -0,0 +1,13 @@
+run_accession	sample_id	family	phenotype_class	read_pairs	unique_mapping_pct	duplicates_removed_pct	cpg_methylation_pct	chg_methylation_pct	chh_methylation_pct
+SRR13482027	SS35-R2	35	SS	62692957	53.6	7.43	13.1	0.8	0.7
+SRR13482028	SS35-R1	35	SS	59408893	53.9	7.22	12.9	0.8	0.8
+SRR13482029	SS05-R3	05	SS	60561212	51.5	8.17	14.4	0.8	0.8
+SRR13482030	RR59-R3	59	RR	54048933	53.5	6.26	12.6	1.1	1.1
+SRR13482031	RR59-R2	59	RR	64457182	53.0	7.65	12.6	0.7	0.7
+SRR13482032	SS05-R2	05	SS	62518979	52.0	7.07	12.9	0.8	0.8
+SRR13482033	RR59-R1	59	RR	72642974	52.2	7.87	13.7	0.8	0.7
+SRR13482034	RR52-R3	52	RR	62839334	50.3	7.90	13.7	0.8	0.8
+SRR13482035	RR52-R2	52	RR	54810357	53.4	6.73	11.9	0.7	0.7
+SRR13482036	RR52-R1	52	RR	59967481	53.9	7.28	12.5	0.6	0.6
+SRR13482037	SS35-R3	35	SS	59620738	50.0	6.43	13.4	0.8	0.8
+SRR13482038	SS05-R1	05	SS	69404859	50.9	9.39	13.9	0.8	0.8

data/studies/CIBNOR2021_HEAT_RRSS/replicate_dispersion_qc.tsv

@@ -0,0 +1,5 @@
+family	class	n_tiles	median_between_library_variance_over_binomial
+05	SS	13668	8.10
+35	SS	13063	6.72
+52	RR	12453	6.83
+59	RR	13731	7.22

registry/studies/CIBNOR2021_HEAT_RRSS.yaml

@@ -17,10 +17,18 @@
 # (2021-01-09), following IOCAS2022_OA_ENERGY. Rename it only if a methylation
 # publication is found.
 #
-# Raw reads are public; no results have been computed yet. results_file is
-# null until the methylation workflow is run (see
-# config/CIBNOR2021_HEAT_RRSS.config), and that run should not happen until the
-# DMR step models family (see limitations).
+# REANALYZED 2026-10-10..16 on Chlo: AREE's first complete real methylation
+# run. All 12 libraries were aligned with Bismark (--score_min L,0,-0.6;
+# 50.0-53.9% unique mapping), deduplicated, and extracted with the read-2
+# end-repair bias trimmed (--ignore_r2 20). Then CpGs were destranded, pooled by
+# FAMILY, kept at >= 10x in all 4 families, summed into 1 kb tiles (>= 10
+# CpGs), and tested with limma's moderated t, 2 RR vs 2 SS families
+# (methylation.dmr_engine = 'limma', dmr_features = 'genes', PR #67): one unit
+# per gene (body + 1.5 kb promoter, 31,166 genes) plus 101,421 intergenic
+# 1 kb tiles. RESULT: none significant (smallest q = 0.57). At 2 families per class this design
+# cannot distinguish phenotype-class differences from family-to-family
+# variation, so this is a null at very low power, not evidence of no effect.
+# See limitations (1) and (9).
 
 study_id: CIBNOR2021_HEAT_RRSS
 simulated: false
@@ -55,8 +63,8 @@ data_availability:
   status: "available"
   access_notes: "Raw paired-end FASTQ for all 12 runs are public in ENA/SRA under PRJNA690951 (first public 2022-02-09). No methylation calls, DMR table, or methylation publication is available. Sample sheet, methylKit sample sheet, and FASTQ manifest with ENA MD5s were generated from ENA metadata; sequence data are not committed. ENA library_selection is RANDOM and the sample alias is 'RAC-Cg-Methyl-Seq'; the bisulfite library kit is not stated."
 
-qc_status: "not_started"
-analysis_status: "not_started"
+qc_status: "passed_with_warnings"
+analysis_status: "harmonized"
 
 # ---------------------------------------------------------------------------
 # One comparison: RR vs SS families, both after 30 d of the same challenge.
@@ -109,7 +117,7 @@ comparisons:
     resilience_classification: "resilience"
     sample_size: 12
     biological_replicates: 2
-    results_file: null
+    results_file: "data/studies/CIBNOR2021_HEAT_RRSS/CIBNOR2021_HEAT_RRSS_rr_vs_ss_day30_oscillating_heat_dmr_standardized.tsv.gz"
     meta_analysis_primary: true
 
 quality_flags:
@@ -123,14 +131,21 @@ limitations: >
   (1) Pseudo-replication risk: thermal tolerance is a family property and each
   class has only two families, so the effective sample size is 2 vs 2. The
   current DMR_METHYLKIT step fits RR vs SS with no family term and would treat
-  6 vs 6 libraries as independent, overstating precision. Family must be
-  modeled, or libraries aggregated to family, before any result is pooled.
+  6 vs 6 libraries as independent, overstating precision. The registered
+  result therefore pools each family's three libraries and tests 2 vs 2
+  families (limma moderated t on family M-values). That is the honest error,
+  and it leaves almost no power: limma's prior df is 3.3, so most tiles rest
+  on their own 2 residual df.
   Family genetic background is fully confounded with phenotype class, so
   family-specific methylation (including cytosine variation from SNPs, which
   bisulfite data cannot separate from C-to-T methylation calls without
   genotype-aware filtering) can look like a phenotype effect. (2) All 12
-  libraries share one BioSample; the deposit does not establish that
-  replicates are distinct animals. (3) There is no ambient-temperature control,
+  libraries share one BioSample, so the deposit does not establish that
+  replicates are distinct animals. The data are consistent with distinct
+  animals: within each family, tile methylation varies between libraries
+  6.7-8.1x more than binomial sampling predicts (replicate_dispersion_qc.tsv;
+  1 kb tiles on NC_088853.1 at >= 30x in all three libraries). Technical
+  replicates would give about 1. (3) There is no ambient-temperature control,
   so the comparison is tolerant vs susceptible genotype under heat, not heat
   response. (4) The phenotype is assigned by pedigree and breeding-program
   performance; no organismal outcome for these animals is deposited, and the
@@ -143,6 +158,32 @@ limitations: >
   to GCF_963853765.1 and must be lifted over, not relabeled, if the reference
   changes. (8) Until a second real thermal methylation study is registered
   (candidate: PRJNA562805), this study cannot contribute to any k >= 2 pool.
+  (9) No unit is significant at FDR 0.05 (smallest q 0.572). 809 of 31,166
+  genes and 3,363 of 101,421 intergenic tiles differ by >= 10 points between
+  classes. Units are genes (all CpGs in body + promoter) so that each gene
+  gives one record and links to RNA-seq. An earlier tile-only pass (319,088
+  1 kb tiles, smallest q 0.336) gave several records per gene, which collide
+  in meta-analysis. The 17 "DMRs" from
+  a first family-level methylKit pass (2026-10-16) were artifacts: with 2
+  families per arm, methylKit's per-region overdispersion collapsed to the
+  binomial for 296 of its 397 q < 0.05 tiles. Unadjusted p < 0.001: 25 genes
+  (~31 expected by chance, so no excess at gene level) and 197 intergenic
+  tiles (~101 expected), so what signal there is looks local and mostly
+  outside genes. Every unit carries a standard error
+  (se_meth_diff_percent = |diff / moderated t|) and an unadjusted p-value, so
+  it can pool: use the effect sizes and standard errors, not the per-study
+  calls. (10) Methylation calls
+  ignore read-2 positions 1-20 and the last base of each read (M-bias:
+  read-2 CpG methylation 84% at position 1, settling to ~14% by 15-20).
+  Without this, CHH methylation read 2.9%; after it, 0.6-1.1%, i.e. >= 99%
+  bisulfite conversion. (11) Unique mapping is 50-54% even with relaxed
+  scoring (13-15% at Bismark's default), so heterozygous or diverged regions
+  are under-represented. (12) The manifest records an empty Bismark version
+  for deduplication and extraction (a version-string mismatch, fixed in
+  PR #67 after this run). All three steps used the same Bismark 0.24.2. To
+  keep the deduplication and extraction cache, the run held those two module
+  files at their pre-fix version (faecd31~1). They differ from PR #67 only in
+  that version line.
 
 provenance:
   registered_by: sr320
@@ -168,6 +209,25 @@ provenance:
     Because sr320 is a co-author on that paper, the replicate identity and
     phenotyping questions in the limitations can be put directly to the
     corresponding author (C. Escobedo-Fregoso, CIBNOR); see issue #6.
+    REANALYSIS (Chlo, 2026-10-09..16; draft PR #14): 24 FASTQ (129 GB),
+    all matching ENA MD5. Twelve resumed transfers came out corrupt at full
+    size and were re-fetched. Nextflow session
+    cdc10126-dbdf-465e-b876-71efcae13f47. First real run of
+    workflows/methylation; first-run fixes went to main in PR #15 (zcat on
+    macOS, CpG-only extraction, .fa genome name, --basename vs --multicore,
+    per-study Bismark options) and PR #67 (limma/DSS pooled-tile engines,
+    report table cap, version capture). Tools: Bismark 0.24.2, bowtie2
+    2.5.5, Trim Galore 0.6.10 (cutadapt 5.2), limma 3.66.0, R 4.5.1. ALIGNMENT
+    SCORE: a 50k-pair test on SRR13482032 mapped 15.5% of pairs at the
+    default --score_min L,0,-0.2, 53.9% at L,0,-0.6, 15.5% non-directional
+    and 15.8% with -X 1000 (bismark_scoring_diagnostic.txt). The library is
+    directional, and the threshold was the cause. The full default run
+    (13.3-14.0%) was discarded. DMR TEST: chosen on synthetic 2-vs-2-family
+    data with real family-level variance (workflows/methylation/README.md).
+    limma was calibrated (1/380 false positives, 5.0% of null p < 0.05); DSS
+    gave 57/380 false positives; methylKit MN collapses. Per-library QC is in
+    library_qc.tsv. Decisions by sr320: native toolchain, relaxed scoring and
+    realignment (2026-10-12), DSS-or-better pooled test (2026-10-16).
   source_links:
     - "https://www.ncbi.nlm.nih.gov/bioproject/PRJNA690951"
     - "https://www.ebi.ac.uk/ena/browser/view/PRJNA690951"