Batrachochytrium read depth, heterozygosity and parent-of-origin analysis

datasets/reads/sketch_store/reports/batrachochytrium_dendrobatidis/coverage_heterozygosity/ancestry_pca.md

@@ -0,0 +1,19 @@
+# Bd isolate ordination (PCA) from SNP allele fractions (2026-11-14)
+
+Status: exploratory. Script `scripts/bd_snp_pca.py` (test in `tests/test_bd_coverage.py`); tables `tables/pca_scores.tsv`, `tables/pca_variance.tsv`.
+Input: the 69 per-isolate SNP tables against CLFT044. Method: biallelic SNPs on the 14 scaffolds of 500 kb or more, DP >= 10 and within 0.5x to 2x of the
+isolate's median depth, every 10th position, sites carried by 3 or more isolates and by at most 90% of isolates (15,804 sites). Entry = alternative-allele
+fraction; an isolate with no SNP row at a site gets 0. Mean-centred, SVD. This is allele-fraction space, not a tree.
+
+## Result
+- PC1 holds 47.3% of the variance (PC2 to PC4: 4.3, 2.8, 2.4%). It is the CLFT-type against JEL-type axis: CLFT-type isolates at +26 to +52, JEL-type mostly at -23 to -9.
+- PC4 separates two isolates from all others: KRBOOR_331 (South Korea) and the Swiss isolate, both at about -31 (all others between -4 and +4). Both have high depth (171x and 81x), so this is not a depth effect. This is the third-lineage-like signal behind their low JEL423 recovery in the parent-of-origin analysis.
+- CLFT160 (HET-rich; 26x) lies at PC1 +16.5, between the class centres (JEL-type median -17.1, CLFT-type median +51.4), about halfway, and at PC4 +2.2, so not on the Swiss/Korean axis. JEL-type isolates at similar depth (22 to 30x) lie at -21.8 to -7.1, so the position is not a depth effect. This agrees with a mixture of the CLFT-type and JEL-type lineages and does not fit the third lineage.
+- CLFT039 (HET-rich; 15x) lies at PC1 +21.1 and PC4 -1.05. This position cannot be read: see the depth effect below.
+
+## The depth effect (a limit of this ordination)
+A missing SNP row is scored as the reference allele, so low depth pulls an isolate toward the middle of PC1. Within the JEL-type class, depth and PC1 correlate with Spearman rho = -0.78 (n = 52). JEL-type isolates of 17x or less have a median PC1 of -1.2 (range -8.6 to +23.6) against -17.6 at 25 to 35x and -20.8 at 60x or more. CLFT039 (15x) falls inside the range of those low-depth isolates, so the ordination does not separate it from them; six other isolates of 13 to 16x sit between the class centres for the same reason and are not evidence of mixture.
+
+## What this does not show
+- A tree, a formal ancestry model (no admixture proportions), or a third-lineage reference. The statement "the lineage mixture fits CLFT160 and not the Swiss/Korean third-lineage signal" rests on one isolate at adequate depth.
+- Anything about CLFT039 beyond the parent-of-origin windows (see `report.md`), which are not depth-corrected either.

datasets/reads/sketch_store/reports/batrachochytrium_dendrobatidis/coverage_heterozygosity/aneuploidy.md

@@ -0,0 +1,36 @@
+# Bd chromosome depth changes and allele balance (2026-11-14)
+
+Status: exploratory; nothing is curated. Source: `bd_depth_by_seq_JEL423.tsv` (69 isolates against the 20-sequence JEL423 assembly, job set of
+2026-11-13), scripts `scripts/bd_aneuploidy.py` and `scripts/bd_aneuploidy_af.py` (tests in `tests/test_bd_coverage.py`). Tables in `tables/`:
+`aneuploidy_flags.tsv`, `aneuploidy_af.tsv`, `genomewide_allele_balance.tsv`, `aneuploidy_by_sequence.tsv`.
+
+## Method
+Sequences of 500 kb or more (16 chromosomes). The depth ratio of a sequence (median depth over the genome median) is divided by the median ratio of the
+same sequence over the isolate's lineage class (CLFT-type 15 isolates, JEL-type 52; the two HET-rich isolates are too few to form a class). This removes
+lineage-specific offsets. An isolate-sequence is flagged when the normalised ratio is outside 1 +/- 0.25 and the isolate has a median depth of 20x or more.
+Then, for each flag, the heterozygous SNPs (minor allele fraction 0.15 to 0.5) of that sequence are compared with those of the isolate's other chromosomes:
+the share near 1/3 (0.27 to 0.40) and near 1/2 (0.43 to 0.50). An extra copy of a chromosome in a diploid should move sites from 1/2 to 1/3.
+
+## Result
+- Population-wide offsets (not isolate-specific): median ratio 1.28 on CP161924.1 and 1.23 on CP161925.1 across all isolates, 0.77 on CP161938.1. They are properties of the reference or of the species relative to JEL423, not aneuploidy. The small sequences CP161939 to CP161942 (57 to 141 kb) have ratios of 0, 0, 3.9 and 0.73: repeats or organelles, excluded from the analysis.
+- 48 isolate-sequence flags in 28 of 69 isolates: 35 gains (norm > 1) and 13 losses; 10 gains at 1.4 or more, 2 losses at about 0.52.
+- Gains (24 with norm 1.3 or more): 15 show the allele-balance shift (1/3 share up and 1/2 share down by more than 0.1 against the isolate's other chromosomes); 9 do not. The clearest cases are in high-depth isolates (50x or more): L2203 CP161926.1 (ratio 1.46; 66% of het sites near 1/3 against 22.5% on other chromosomes, 7% near 1/2 against 50%), KRBOOR_331 CP161926.1 (1.45; 76% against 22%) and RTP5 CP161931.1 (1.38; 68% against 24%). That is the signature of one extra copy of one chromosome in a diploid isolate. The same isolates have 50 to 57% of het sites near 1/2 on their other chromosomes, which is the diploid pattern.
+- CLFT065 CP161936.1 (ratio 1.59, 55x) shows the opposite shift (53% of het sites near 1/2 against 30% on the rest), so it is not the trisomy pattern; I cannot say what it is (a duplication, or a different copy state).
+- Losses (7 with norm 0.7 or less): in 4 the het SNP density on the flagged chromosome is below 0.75 of the rest, as expected when one copy is lost (CLFT078 and CLFT166 CP161938.1 at ratio about 0.52: 831 and 698 het per Mb against 1,866 and 1,078).
+- CP161926.1 is gained in 6 isolates from different lineage classes and countries (Montserrat L2203, South Korea KRBOOR_331, Sweden, and three Brazilian isolates); CP161936.1, CP161937.1 and CP161931.1 are gained in 4 isolates each. I did not test whether these gains are shared by descent.
+- Isolates with several flagged chromosomes: CLFT065 (5), CLFT080, CLFT156, CLFT165, CLFT167 (3 each) and the "Generic sample" from Switzerland (4). I did not examine them further.
+
+## What this does not show
+- A baseline ploidy for each isolate. The genome-wide allele balance puts 21 isolates in a "third-dominated" group, 15 "half-dominated" and 33 "mixed", but this label is confounded: the third-dominated isolates have a median depth of 21x against 105x for half-dominated ones. At low depth the allele fractions are noisy and reference divergence skews them, so I do not read the label as triploidy. Only the within-isolate comparison (flagged chromosome against the same isolate's other chromosomes) is controlled for depth.
+- A whole-genome ploidy change is invisible in a ratio to the genome median.
+- The 20x cutoff leaves noise at the lowest depths; 10 isolates are below 20x and were not flagged.
+- The two HET-rich isolates are covered in the next section; CLFT039 has a median depth of 16.2x, below the 20x cutoff, so it cannot be flagged.
+
+## The two HET-rich isolates (CLFT160 = SRR10698391, CLFT039 = SRR10698425): hybrid-like or aneuploid?
+`scripts/bd_aneuploidy.py --fallback-class JEL-type` compares them with the JEL-type median (two isolates cannot form a class; the JEL-type baseline is a rough
+reference). `scripts/bd_hetrich_check.py` gives per-sequence depth, allele balance (JEL423 reference) and parent-of-origin window classes (CLFT044 reference).
+Tables: `aneuploidy_flags_with_hetrich.tsv`, `hetrich_by_sequence.tsv`.
+- Depth: CLFT160 has two gain flags (CP161931.1 at 1.32 and CP161936.1 at 1.38 against the JEL-type median); CLFT039 (16.2x) is below the flag cutoff, and its ratios lie between 0.77 and 1.24. Neither isolate has a loss flag.
+- Parent-of-origin windows on the 14 scaffolds of 500 kb or more: heterozygous windows (JEL423 and CLFT044 alleles both present) are the majority on almost all of them: median share 0.70 in CLFT160 and 0.72 in CLFT039. In CLFT039 the share is 0.47 or more on every scaffold; in CLFT160 it is 0.40 or more on 12 of 14. The two exceptions are JAQIHF010000010.1 (share 0.105; 32% CLFT windows) and JAQIHF010000014.1 (share 0.000; 50% CLFT windows); their depth is normal (0.99, 0.93), so they are not copy-number changes. Windows homozygous for the JEL423 allele are almost absent (median share 0.01 and 0.05); CLFT-class windows are 29% and 20%.
+- So the heterozygosity of both isolates is genome-wide and is not tied to a gained chromosome: CLFT160's two gains are on scaffolds with ordinary HET shares. The pattern fits an isolate that carries both lineages in one genome with part of the genome homozygous for the CLFT044 allele (possible loss of heterozygosity on two scaffolds of CLFT160), and does not fit aneuploidy as the cause of the heterozygosity.
+- Limits: the allele-balance columns of the same table (median 0.37 and 0.33 of het sites near 1/2) are at depths of 27.5x and 16.2x and are too noisy to read; the "hybrid" reading rests on the window classes, which depend on the 0.5 recovery cutoff set from two control isolates; a third lineage that shares alleles with both references would give the same pattern; the two isolates are the only ones of this type, so no comparison group exists.

datasets/reads/sketch_store/reports/batrachochytrium_dendrobatidis/coverage_heterozygosity/idiomorph_scan.md

@@ -0,0 +1,28 @@
+# Bd cross-isolate depth scan for a mating-type idiomorph (2026-11-14)
+
+Status: negative result. No candidate idiomorph was found. Nothing here is curated.
+
+Script: `scripts/bd_idiomorph_scan.py` (tests in `tests/test_bd_coverage.py`). Input: the 10 kb window tables of 58 (CLFT044
+reference) and 59 (JEL423 reference) isolates with median depth of at least 20x. Genes come from the BFD funannotate GFF3
+(scaffold IDs identical to the references used). Tables: `tables/idiomorph_regions_<ref>.tsv`, `tables/idiomorph_null_<ref>.tsv`.
+
+## Method
+Each window ratio is divided by its sequence median, which removes whole-sequence copy-number change. A window is FULL
+(0.8 to 1.25), HALF (0.35 to 0.65) or LOW (below 0.25) per isolate. A candidate window has at least 5 FULL isolates and at least
+3 LOW or HALF isolates. Neighbouring candidate windows merge; regions of 20 kb or more are kept. Expected signature of a
+diploid idiomorph: FULL, HALF and LOW isolates in the same window, with few het SNPs in HALF isolates (only one copy aligns).
+
+## Result
+- CLFT044 coordinates: 71 regions. JEL423 coordinates: 70 regions.
+- Most regions are HALF in many isolates, with a median of 20 to 60 het SNPs per 10 kb window in those isolates. That fits
+  divergent sequence that maps poorly onto a reference of another lineage, not a hemizygous region.
+- Regions with n_half >= 3 and median het in HALF isolates <= 3: 3 against CLFT044 and 5 against JEL423. None has a gene whose
+  product or note matches HMG, homeodomain, mating, sex or pheromone. HALF isolates in these regions come from both CLFT-type and
+  JEL-type groups, so lineage does not explain them. I did not examine their gene content further.
+- The rotation null (each isolate rotated by a random offset, 100 permutations) gave a mean of 120 (CLFT044) and 110 (JEL423)
+  regions, more than observed. The region count is therefore not a usable test statistic, and I make no significance claim.
+
+## Limits
+- No evidence that Bd has a two-idiomorph system. A scan like this assumes a presence/absence or hemizygous idiomorph.
+- 10 kb windows. A small locus can be missed. Smaller windows need a new alignment, because no BAMs were kept.
+- Isolates below 20x were dropped. Two of the references are from different lineages, which drives the HALF windows.

datasets/reads/sketch_store/reports/batrachochytrium_dendrobatidis/coverage_heterozygosity/report.md

@@ -0,0 +1,101 @@
+# Batrachochytrium dendrobatidis: read depth, heterozygosity and parent of origin
+
+Status: exploratory analysis of public reads. Nothing here is curated or accepted. Bd isolates are diploid, and the
+allele balance is read that way.
+
+Pipeline commit: `6bfc6a3cad`. Jobs: JEL423 reference 32070127 and 32070630 (array 2-9); CLFT044 reference 32070128 and
+32070631; reference alignment 32089404 (minimap2 2.30, `-x asm10 -c --cs=short`). Result folder on the HPCC:
+`results/2026-11-13_bd_coverage/` (not committed; tables below are copies).
+
+## Data
+
+- Runs: 69 runs analysed of 72 listed (Brazil 52, other countries 17), chosen from PopSeq
+  `catalog/runs/fungi_all_runs.csv.gz` (ENA snapshot 2026-11-06) and ENA. All were Illumina WGS. See `tables/bd_runs_popseq.tsv`.
+- References: JEL423 (GCA_048537975.1) and UCR_BdCLFT044 (Brazil lineage per the curator). I did not verify the lineage of JEL423 here.
+- Each run: `prefetch`, `fasterq-dump`, read-pair count checked against ENA `read_count`, `bwa-mem2`, duplicate removal,
+  `mosdepth` 10 kb windows (MAPQ 20), `bcftools` SNP calls (QUAL 30).
+- Excluded runs (`tables/excluded_runs.tsv`):
+  - SRR13660529: mate2 7132682 vs ena 29477484
+  - ERR14883802: mates 4218822 each; ena 8437644 (= 2x mate count)
+
+## 1. Depth and SNP density
+
+Median depth of the 69 isolates: 32.9x (range
+10.8 to 178.6). 10 isolates have median depth below 20x; treat their numbers with more care.
+
+Isolates against their own reference (controls):
+
+| control | SNPs/Mb vs JEL423 | het fraction | SNPs/Mb vs CLFT044 | het fraction |
+|---|---|---|---|---|
+| CLFT044 reads | 5399.1 | 0.4391 | 1356.8 | 0.9345 |
+| JEL423 reads | 1851.2 | 0.9264 | 4611.6 | 0.4079 |
+
+An isolate does not give near-zero SNPs against its own reference. Each reference is a haploid representation of a diploid
+genome, so heterozygous sites are expected, and the het fraction is high (0.9). Repeats and paralogs can add het-like SNPs;
+I did not separate the two causes. A het fraction alone therefore cannot call a hybrid.
+
+Depth flags: sequences of 200 kb or more with median depth outside 0.75 to 1.25 of the genome median occur in most isolates (0 to 11 per isolate)
+(see `tables/bd_isolate_summary_*.tsv`, column `depth_flag_seqs`). I did not separate copy-number change from repeat or
+assembly artefacts. Ploidy variation, including aneuploidy, is also known in Bd; I did not test it here.
+
+![SNP density and parent-of-origin summary](figures/bd_summary.png) [PDF](figures/bd_summary.pdf)
+
+Figure source: `tables/isolate_classes.tsv`, `tables/parent_origin_summary.tsv`; script `scripts/bd_report.py`.
+
+## 2. Parent of origin
+
+Method (`scripts/bd_parent_of_origin.py`): align JEL423 to CLFT044 and take single-base differences in primary alignments
+(172,706 sites). In each isolate, a site is "called" when its CLFT044-reference SNP table has the JEL423 base. Per 50 kb
+window: recovery = called sites / usable sites; het share = called sites with alt fraction 0.2 to 0.8. A window is `CLFT`
+when recovery < 0.5, else `HET` when het share >= 0.5, else `JELHOM`.
+
+Why sites without a SNP row are not counted as CLFT044 allele: in the JEL423 isolate about a third of divergent sites have no row,
+so absence is not evidence. I first counted them as CLFT044 and the JEL423 isolate then scored only 0.57 JEL fraction; that
+version was dropped.
+
+Calibration: only two controls. In the CLFT044 isolate the window recovery has median 0.05 and 95th percentile 0.41. In the
+JEL423 isolate the median is 0.87 and the 5th percentile 0.16. I set the cutoff at 0.5 from these. It is not validated on isolates
+of known hybrid status. Isolate class: CLFT-type if over 60% CLFT windows, JEL-type if over 60% JELHOM windows, HET-rich if
+over 50% HET windows, else mixed.
+
+| group | CLFT-type | JEL-type | HET-rich | mixed |
+|---|---|---|---|---|
+| Brazil (n=52) | 13 | 37 | 2 | 0 |
+| Other countries (n=17) | 2 | 15 | 0 | 0 |
+
+Control isolates: CLFT044 is CLFT-type; JEL423 is JEL-type.
+
+### Isolates that are HET-rich or mixed
+
+| isolate | country | median depth | SNPs/Mb vs JEL423 | SNPs/Mb vs CLFT044 | recovery | het share of called | windows CLFT/HET/JELHOM | class |
+|---|---|---|---|---|---|---|---|---|
+| CLFT160 | Brazil | 27.5 | 4593.9 | 3939.5 | 0.5369 | 0.8059 | 152/277/28 | HET-rich |
+| CLFT039 | Brazil | 16.2 | 3874.4 | 3429.5 | 0.5664 | 0.7422 | 102/292/38 | HET-rich |
+
+These isolates carry JEL423 alleles in most windows, mostly in one copy. That is the pattern expected for a first-generation
+hybrid of the two lineages. It does not show that they are hybrids. Other explanations I did not exclude: a third lineage that shares
+alleles with both references; aneuploidy; mapping effects at divergent sites. A hybrid call needs a tree or ancestry analysis and
+a lineage reference beyond these two. Published context: a Bd-GPL x Bd-Brazil hybrid zone in the Serra da Graciosa, Parana (Jenkinson et al. 2016,
+doi 10.1111/mec.13599). I did not compare my isolates with that study's samples.
+
+### Non-Brazil isolates that are CLFT-type
+
+| isolate | country | median depth | SNPs/Mb vs JEL423 | SNPs/Mb vs CLFT044 | recovery | het share of called | windows CLFT/HET/JELHOM | class |
+|---|---|---|---|---|---|---|---|---|
+| KRBOOR_331 | South Korea | 175.0 | 6472.1 | 5103.9 | 0.3672 | 0.4199 | 350/8/122 | CLFT-type |
+| Generic sample from B. dendrobatis | Switzerland: Gamlikon | 86.2 | 5091.5 | 4115.0 | 0.3267 | 0.3347 | 377/1/95 | CLFT-type |
+
+Their recovery is 0.33 and 0.37, well below the JEL-type isolates (about 0.6 to 0.75). They carry few alleles of either reference.
+A third lineage is a possible reason; I did not check their lineage.
+
+## Limits
+
+- Two references only. Lineage names are the curator's for CLFT044 and my understanding for JEL423.
+- SNP tables keep QUAL >= 30 and no per-site depth was stored, so absence of a call cannot be told from no coverage.
+- Window class is noisy: the CLFT044 control has windows with recovery up to 0.6.
+- The Brazil set is 52 runs; PopSeq lists 58 Brazilian isolates, so the set is not all of them.
+
+## Next
+
+Add an ancestry or tree analysis with at least one more lineage reference (for example a BdASIA or BdCAPE assembly), and test
+the HET-rich isolates for aneuploidy with allele balance by sequence.