Gc bias genomes

.github/workflows/test.yml

@@ -25,11 +25,15 @@ jobs:
           micromamba-version: '1.5.10-0' # any version from https://github.com/mamba-org/micromamba-releases
           post-cleanup: 'all'
 
-      - name: Setup Nextflow latest-edge
+      # Pinned, not 'latest': main.nf needs `nextflow.preview.topic`, which 25.x removed once
+      # topic channels became stable. 24.10.5 is what the cluster runs, so CI matches production.
+      - name: Setup Nextflow
         uses: nf-core/setup-nextflow@v2
         with:
-          version: latest
+          version: "24.10.5"
 
+      # Pinned: nf-test serialises the workflow.trace map in a version-dependent key order, so
+      # every tests/*.snap must be generated with the same version this installs.
       - name: Install nf-test
         uses: nf-core/setup-nf-test@v1
         with:

README.md

@@ -88,6 +88,7 @@ nextflow run main.nf \
 | `--outputDir` | Output directory | `em-seq_output` |
 | `--enable_neb_agg` | Enable NEB aggregation reporting | `False` |
 | `--adapter_set` | adapter set found in conf/adapter_sequences.config (`Illumina`, `Element`, `MGI`,  `Nextera-mosaic`) | Required |
+| `--multiqc_gc_group` | contig group whose GC bias curve MultiQC plots (see [GC bias curves](#gc-bias-curves)) | largest group |
 
 ### References Config
 
@@ -96,6 +97,28 @@ Modify the conf/references.config file to specify your genome files
 - `genome_fai` path to your genome fasta fai file 
 - `bwameth_index` path to your genome fasta file where bwameth indices exist 
 - `target_bed` BED file for targeted analysis, Optional 
+- `gc_groups_dir` directory of per-organism GC bias references, Optional (see below)
+
+### GC bias curves
+
+Picard measures GC bias by comparing where reads start against the GC of every window in the
+reference it is given, so a composite reference reports one curve blended across every organism in
+it. Setting `gc_groups_dir` for a genome adds a separate curve per organism alongside the
+whole-reference one. The directory holds a `contig_groups.tsv` assigning each contig to a group,
+plus a `<group>.fa`/`.fa.fai`/`.dict`/`.bed`/`.contigs.txt` per group; it is validated against
+`genome_fai` at startup.
+
+MultiQC plots a single curve per library. By default that is the curve of the group with the most
+windows -- the host organism for the `+meth_controls` references -- because the whole-reference curve
+counts the spike-in controls, and pUC19 (51% GC), lambda (50%) and Xp12 (68%) pull it noticeably
+away from a host-only curve. Use `--multiqc_gc_group <name>` to plot a different group.
+
+Without `gc_groups_dir` you get the whole-reference curve alone, controls included; the run warns
+that it is not directly comparable to a host-only curve.
+
+Group definitions for the published references ship in
+[`assets/contig_groups/`](assets/contig_groups/); the per-group FASTA subsets they refer to are
+built separately and are not part of this repository.
 
 ### Advanced Options
 - `--tmp_dir` - Temporary directory (default: `/tmp`)
@@ -104,10 +127,35 @@ Modify the conf/references.config file to specify your genome files
 
 
 ## Reference Genomes
-Pre-built reference genomes with methylation spike-in controls + CpG Island Annotation files:
-- **T2T CHM13**: [T2T_chm13v2.0+  meth_controls](https://neb-em-seq-sra.s3.us-east-1.amazonaws.com/T2T_chm13v2.0%2Bmeth_controls.fa.gz) [CpG Islands](https://neb-em-seq-sra.s3.us-east-1.amazonaws.com/human_t2t_cpg_islands.gtf)
-- **GRCh38**: [GRCh39 + meth controls](https://neb-em-seq-sra.s3.us-east-1.amazonaws.com/grch38_core%2Bmeth_controls.fa.gz) [CpG Islands](https://neb-em-seq-sra.s3.us-east-1.amazonaws.com/human_grch38_cpg_islands.gtf)
-- **GRCm39**: [GRCm39 + meth_controls](https://neb-em-seq-sra.s3.us-east-1.amazonaws.com/grcm39%2Bmeth_controls.fa.gz) [CpG Islands](https://neb-em-seq-sra.s3.us-east-1.amazonaws.com/grcm39_cpg_islands.gtf)
+
+Pre-built reference genomes with methylation spike-in controls, plus matching UCSC CpG island
+annotations. These are the same references NEB uses internally:
+
+| Species | Assembly | Genome FASTA | CpG islands gtf |
+|---|---|---|---|
+| Human | T2T CHM13v2.0 | [T2T_chm13v2.0+meth_controls](https://neb-em-seq-sra.s3.us-east-1.amazonaws.com/T2T_chm13v2.0%2Bmeth_controls.fa.gz) | [t2t_cpg_islands](https://neb-em-seq-sra.s3.us-east-1.amazonaws.com/human_t2t_cpg_islands.gtf) |
+| Human | GRCh38 | [grch38_core+meth_controls](https://neb-em-seq-sra.s3.us-east-1.amazonaws.com/grch38_core%2Bmeth_controls.fa.gz) | [grch38_cpg_islands](https://neb-em-seq-sra.s3.us-east-1.amazonaws.com/human_grch38_cpg_islands.gtf) |
+| Mouse | GRCm39 | [grcm39+meth_controls](https://neb-em-seq-sra.s3.us-east-1.amazonaws.com/grcm39%2Bmeth_controls.fa.gz) | [grcm39_cpg_islands](https://neb-em-seq-sra.s3.us-east-1.amazonaws.com/grcm39_cpg_islands.gtf) |
+
+Download and unzip the FASTA, then add its paths to `conf/references.config` (see
+[References Config](#references-config)). The pipeline reads the reference from disk and never
+fetches it, so it also needs a `.fai` and a bwameth index alongside; `main.nf` checks for both at
+startup and prints the command to build whichever is missing.
+
+The CpG island files are only used by `feature_cov_meth.nf`. Pass one with
+`--ucsc_cpg_islands_gtf`, or save the downloads under `--local_ref_files_path` (default
+`~/nebnext_projects/em-seq/em-seq_ref_files`) keeping their names, and the `--mouse` and
+`--human_t2t2` shortcuts will find them. GRCh38 has no shortcut, so it always needs these
+parameters passed explicitly.
+
+Per-organism GC bias curves for these references are configured by the group definitions in
+[`assets/contig_groups/`](assets/contig_groups/) — see [GC bias curves](#gc-bias-curves).
+
+> If you already have one of these references from before September 2026, re-download it. The
+> `phage_lambda`, `phage_T4` and `phage_Xp12` control sequences were corrected to match RefSeq, so
+> an older copy gives control metrics that are not comparable. You must rebuild the bwameth index
+> and the `.fai` after replacing the FASTA.
+
 - Create your own reference by appending the [control sequences](assets/methylation_controls.fa) to your preferred genome fasta (e.g. `cat genome.fa methylation_controls.fa > genome+methylation_controls.fa`)
    | Sequence | Methylation State                                | Purpose |
    |----------|--------------------------------------------------|---------|

assets/contig_groups/README.md

@@ -0,0 +1,68 @@
+# Contig group definitions
+
+One TSV per reference, assigning every contig to a group so the pipeline can report a GC bias
+curve per organism instead of a single curve blended across the whole composite. Setting a
+genome's `gc_groups_dir` in `conf/references.config` turns the feature on; see **GC bias curves**
+in the top-level README.
+
+These cover the references linked under [Reference
+Genomes](../../README.md#reference-genomes), and are checked against the genome's `.fai` at
+startup, so they only apply to those exact FASTAs.
+
+| Column | Meaning |
+|---|---|
+| `contig_name` | must match the reference `.fai` exactly |
+| `sequence_length` | must match the `.fai`; a mismatch means the TSV was built against a different reference |
+| `frac_gc` | GC over ACGT positions only, so an N-heavy contig is not reported as low-GC. Informational — the pipeline does not read it |
+| `contig_group` | the curve this contig contributes to. `[A-Za-z0-9_.-]+`. Blank means no curve of its own |
+| `contig_type` | one of `autosomal`, `sex_chromosome`, `organelle`, `unplaced`, `control`, `unknown` |
+
+Groups are named `<organism>_autosome`, `<organism>_sex`, `<organism>_mito` and
+`<organism>_unplaced` for the host, and `lambda`, `pUC19`, `T4`, `Xp12` (plus `EBV` for the human
+references) for the non-host contigs:
+
+| Reference | Groups |
+|---|---|
+| `T2T_chm13v2.0+meth_controls` | human_autosome (22), human_sex (2), human_mito, EBV, lambda, pUC19, T4, Xp12 |
+| `grch38_core+meth_controls` | human_autosome (63), human_sex (3), human_mito, human_unplaced (127), EBV, lambda, pUC19, T4, Xp12 |
+| `grcm39+meth_controls` | mouse_autosome (58), mouse_sex (2), mouse_mito, lambda, pUC19, T4, Xp12 |
+
+## Provenance
+
+These are the group assignments NEB uses internally, and match published references.
+
+One caveat carried over from how the references were built: **`grcm39+meth_controls` names its
+chromosomes by GenBank accession** (`CM000994.3` for chr1, and so on) rather than `chr1`-style,
+unlike the two human references. This is deliberate — it is why `feature_cov_meth.nf`'s mouse
+preset sets `cpg_chr_lookup = '$10,$5'` to translate assembly-report chromosome names to
+accessions, where the T2T preset uses `'$10,$10'`.
+
+## Adding another reference
+
+Write a TSV with the five columns above: one row per contig in the reference's `.fai`, in `.fai`
+order, with `sequence_length` matching the `.fai` exactly. Give every contig that should get its own
+curve a `contig_group`, and leave the column blank for contigs that should not.
+
+Then build one subset reference per group and point that genome's `gc_groups_dir` at the directory
+holding them:
+
+```bash
+TSV=contig_groups.tsv
+FA=genome.fa            # must already have a .fai alongside it
+
+for grp in $(awk -F'\t' 'NR>1 && $4!=""{print $4}' "$TSV" | sort -u); do
+    # contig list and BED are taken in .fai order, which is the order picard requires
+    awk -F'\t' -v g="$grp" 'NR==FNR{if(FNR>1&&$4==g)k[$1];next} $1 in k{print $1}' \
+        "$TSV" "$FA.fai" > "$grp.contigs.txt"
+    awk -F'\t' -v g="$grp" 'NR==FNR{if(FNR>1&&$4==g)k[$1];next} $1 in k{print $1"\t0\t"$2}' \
+        "$TSV" "$FA.fai" > "$grp.bed"
+    samtools faidx "$FA" -r "$grp.contigs.txt" > "$grp.fa"
+    samtools faidx "$grp.fa"
+    samtools dict  "$grp.fa" > "$grp.dict"
+done
+```
+
+Each group's subset FASTA is a copy of those contigs' sequence, so the directory costs roughly as
+much disk as the reference itself. The pipeline validates the TSV against `genome_fai` and checks
+each group's contig order against its subset reference before running picard, so a mistake here
+fails at startup rather than part-way through a run.

assets/contig_groups/T2T_chm13v2.0+meth_controls.contig_groups.tsv

@@ -0,0 +1,31 @@
+contig_name	sequence_length	frac_gc	contig_group	contig_type
+chr1	248387328	0.415003	human_autosome	autosomal
+chr2	242696752	0.402384	human_autosome	autosomal
+chr3	201105948	0.394389	human_autosome	autosomal
+chr4	193574945	0.381309	human_autosome	autosomal
+chr5	182045439	0.395244	human_autosome	autosomal
+chr6	172126628	0.396116	human_autosome	autosomal
+chr7	160567428	0.407184	human_autosome	autosomal
+chr8	146259331	0.402332	human_autosome	autosomal
+chr9	150617247	0.413461	human_autosome	autosomal
+chr10	134758134	0.415390	human_autosome	autosomal
+chr11	135127769	0.415531	human_autosome	autosomal
+chr12	133324548	0.407746	human_autosome	autosomal
+chr13	113566686	0.383205	human_autosome	autosomal
+chr14	101161492	0.410670	human_autosome	autosomal
+chr15	99753195	0.421158	human_autosome	autosomal
+chr16	96330374	0.435331	human_autosome	autosomal
+chr17	84276897	0.452839	human_autosome	autosomal
+chr18	80542538	0.398008	human_autosome	autosomal
+chr19	61707364	0.476621	human_autosome	autosomal
+chr20	66210255	0.437182	human_autosome	autosomal
+chr21	45090682	0.416452	human_autosome	autosomal
+chr22	51324926	0.458491	human_autosome	autosomal
+chrX	154259566	0.395106	human_sex	sex_chromosome
+chrY	62460029	0.361454	human_sex	sex_chromosome
+chrM	16569	0.443720	human_mito	organelle
+chrEBV	171823	0.594967	EBV	control
+phage_lambda	48502	0.498577	lambda	control
+plasmid_puc19c	2686	0.506329	pUC19	control
+phage_T4	168903	0.352984	T4	control
+phage_Xp12	64272	0.681868	Xp12	control