B2: 5' homogeneity as the Simpson index of mature-arm starts

AutoDeep/scripts/folding_properties.py

@@ -352,14 +352,16 @@ def calculate_homogeneity_score(five_prime_df, miRNA_id=None):
 drosha_overhangs = []
 dicer_overhangs = []
 for item in tqdm(loci):
-	homogeneity, loci_count = five_prime_features(mrd_loci[item].reads)
+	homogeneity, loci_count = five_prime_features(mrd_loci[item])
 	five_prime_process_signal_selected_row_count_sum_over_mature_count_sum.append(homogeneity)
 	num_locis.append(loci_count)
 	# Signed 3' overhangs; blank for loci without star reads
 	drosha, dicer = overhangs(mrd_loci[item])
 	drosha_overhangs.append(drosha)
 	dicer_overhangs.append(dicer)
 
+# 5' homogeneity (B2: Simpson index of the mature-arm starts) keeps its old column name
+# until B5 renames the features, since the model selects them by name
 signaling_df = pd.DataFrame({'provisional_id': loci, "five_prime_process_signal_selected_row_count_sum_over_mature_count_sum" : five_prime_process_signal_selected_row_count_sum_over_mature_count_sum,  'num_locis' : num_locis})
 intersected_data = pd.merge(signaling_df, intersected_data, on='provisional_id')
 # "Bot" is the Drosha end (base of the stem), "Top" the Dicer end (loop side). They

AutoDeep/scripts/mrd.py

@@ -17,7 +17,6 @@
 from typing import List, NamedTuple, Optional, Tuple
 
 import numpy as np
-import pandas as pd
 
 READ_LINE = re.compile(r'^(\S+)_x(\d+)\s+([.A-Za-z]+)(?:\s|$)')
 FIELD_LINE = re.compile(r'^(exp|obs|pri_seq|pri_struct)\s+(\S+)')
@@ -79,47 +78,56 @@ def parse_mrd(path):
             for pid, reads in loci.items()}
 
 
-def arm_counts(locus):
-    """Return {'mature', 'star', 'loop'} read counts, as miRDeep2 reports them in result.csv.
-
-    Reads that fit none of the three regions, such as reads spanning two arms,
-    count towards none.
-    """
+def reads_by_arm(locus):
+    """Yield (arm, start, count) for each read, start 1-based, arm None if it fits no region."""
     if locus.exp is None:
         raise ValueError("locus has no exp line")
     regions = [(arm, locus.exp.index(mark) + 1, locus.exp.rindex(mark) + 1)
                for arm, mark in ARMS if mark in locus.exp]
-    counts = {arm: 0 for arm, _ in ARMS}
     for aln, count in locus.reads:
         beg = len(aln) - len(aln.lstrip('.')) + 1
         end = len(aln.rstrip('.'))
-        for arm, first, last in regions:
-            if first - FUZZ_BEG <= beg and end <= last + FUZZ_END:
-                counts[arm] += count
-                break
+        arm = next((name for name, first, last in regions if first - FUZZ_BEG <= beg and end <= last + FUZZ_END), None)
+        yield arm, beg, count
+
+
+def arm_counts(locus):
+    """Return {'mature', 'star', 'loop'} read counts, as miRDeep2 reports them in result.csv.
+
+    Reads that fit none of the three regions, such as reads spanning two arms,
+    count towards none.
+    """
+    counts = {arm: 0 for arm, _ in ARMS}
+    for arm, _, count in reads_by_arm(locus):
+        if arm is not None:
+            counts[arm] += count
     return counts
 
 
-def five_prime_features(reads):
+def five_prime_features(locus):
     """Return (5' homogeneity, number of distinct reads) for one locus.
 
-    Homogeneity is the share of reads that start at the same position and
-    with the same nucleotide as the most abundant read, out of all reads that
-    start no more than 20 nt after it. Ties for most abundant read go to the
-    first in file order. This reproduces the pypdf-based definition exactly.
+    Homogeneity is the unbiased Simpson index of the 5' start positions of the
+    reads on the mature arm (those arm_counts counts as mature): the
+    probability that two of them, drawn without replacement, start at the same
+    position, sum n_i (n_i - 1) / (N (N - 1)) over start positions i. Star and
+    loop reads do not count, so the value does not depend on which arm the
+    mature is on, and unlike the share of the commonest start or Shannon
+    entropy it is not biased by sequencing depth. A read with a mismatched
+    first base counts at its position. NaN with fewer than 2 mature-arm reads.
+
+    The number of distinct reads counts every read on the precursor.
     """
-    if not reads:
+    if not locus.reads:
         raise ValueError("locus has no aligned reads")
-    positions = [len(aln) - len(aln.lstrip('.')) for aln, _ in reads]
-    df = pd.DataFrame({
-        'nucleotide': [aln[p] for (aln, _), p in zip(reads, positions)],
-        'position': positions,
-        'count': [count for _, count in reads],
-    })
-    top = df.iloc[np.argmax(df['count'].to_numpy())]
-    selected = df[(df['nucleotide'] == top['nucleotide']) & (df['position'] == top['position'])]
-    window = df[(df['position'] - top['position']) <= 20]
-    return selected['count'].sum() / window['count'].sum(), df.shape[0]
+    starts = {}
+    for arm, beg, count in reads_by_arm(locus):
+        if arm == 'mature':
+            starts[beg] = starts.get(beg, 0) + count
+    n = np.array(list(starts.values()), dtype=float)
+    total = n.sum()
+    homogeneity = (n * (n - 1)).sum() / (total * (total - 1)) if total > 1 else math.nan
+    return homogeneity, len(locus.reads)
 
 
 def base_pairs(structure):

CLAUDE.md

@@ -41,7 +41,7 @@ We are revising the code in response to a mock peer review. The work list, with
 ### Open decisions (ask Tejus; do not decide)
 
 - Q1, Q2 and Q4-Q9 in REVISIONS.md are unanswered; B2, B4-B8 and C1/C3/C5 depend on them. Q1 (seed homology: labelling criterion or feature) affects the most items.
-- Feature set: since B1 the pipeline produces the shipped model's features (`five_prime_process_signal_...`, `num_locis`). Since B3 the two overhang features are defined differently from the training data (signed, from `output.mrd`, blank without star reads) until B8 rebuilds it. The entropy features in `training_data/` are computed by functions kept in `folding_properties.py` but not written out.
+- Feature set: since B1 the pipeline produces the shipped model's features (`five_prime_process_signal_...`, `num_locis`). Since B2, 5' homogeneity is the Simpson index of the mature-arm starts (old column name until B5), unlike the training data until B8. Since B3 the two overhang features are defined differently from the training data (signed, from `output.mrd`, blank without star reads) until B8 rebuilds it. The entropy features in `training_data/` are computed by functions kept in `folding_properties.py` but not written out.
 
 ### Known issues not yet in REVISIONS.md
 
@@ -56,7 +56,7 @@ We are revising the code in response to a mock peer review. The work list, with
 - `AutoDeep/AutoDeep`: bash wrapper; put `AutoDeep/` on PATH to get the `AutoDeep` command (there is no setup.py)
 - `AutoDeep/scripts/AutoDeep.sh`: command dispatcher. `infer` runs init_data.py, RNAfold, csv_feature_extraction.py, folding_properties.py, boosted_forest_inference.py
 - `AutoDeep/scripts/features.py`: fixed encodings of the flag columns and labels (Candidate 0, Confident 1, falsepositive 2), shared by all scripts (A1)
-- `AutoDeep/scripts/mrd.py`: parses miRDeep2's `output.mrd` (per locus: `exp` and `obs` strings, `pri_seq`, `pri_struct` and reads; IDs without commas, as init_data.py leaves `result_*.csv`: `strip_commas`, A8) and computes 5' homogeneity, `num_locis`, per-arm read counts by miRDeep2's own rule (`arm_counts`; B1) and the signed 3' overhangs at both duplex ends from the terminal mature-star base pairs (`overhangs`, the duplex rule; B3)
+- `AutoDeep/scripts/mrd.py`: parses miRDeep2's `output.mrd` (per locus: `exp` and `obs` strings, `pri_seq`, `pri_struct` and reads; IDs without commas, as init_data.py leaves `result_*.csv`: `strip_commas`, A8) and computes 5' homogeneity (unbiased Simpson index of the mature-arm reads' 5' starts; B2), `num_locis`, per-arm read counts by miRDeep2's own rule (`arm_counts`, from `reads_by_arm`; B1) and the signed 3' overhangs at both duplex ends from the terminal mature-star base pairs (`overhangs`, the duplex rule; B3)
 - `AutoDeep/scripts/hyperparameters.py`: reads `-hp` files (text "name value", JSON, or a run's config.json) without running code (A3)
 - `AutoDeep/scripts/boosted_forest_inference.py`: selects the model's features by name; `--model PATH` (default: bundled model)
 - `AutoDeep/scripts/boosted_forest_training.py`: `AutoDeep train` (Ray Tune + OptunaSearch). Writes only to `--out_dir` (default `AutoDeepRun/models/<timestamp>/`): miRNA_model.json, config.json, metrics.json, user_data.csv, tuning log. `-t` takes several files; `--seed` (default 42) seeds splits, XGBoost and the Optuna sampler. `load_training_set()` and `tuning_search_alg()` are shared with alt_model_test.py
@@ -66,8 +66,8 @@ We are revising the code in response to a mock peer review. The work list, with
 - `AutoDeep/training_data/`: 3,180 loci in 13 `new_data_*` CSVs with entropy-based 5' features (`normalized_entropy`, `normalized_entropy_s`); not used by the pipeline or shipped model
 - `AutoDeep/model_weights/miRNA_model.json`: shipped model (XGBoost 2.1.1, 983 rounds x 3 classes, depth 5)
 - `AutoDeep/AutoDeep.yml`: portable conda spec and the only requirements list (no requirements.txt); `AutoDeep/AutoDeep.linux-64.lock.yml`: exact Linux env the shipped model was trained in
-- `review/`: REVISIONS.md, reanalysis.py (reviewer's evaluation script), reanalysis_baseline.txt (its output on commit 133bf14), reanalysis_osx-arm64.txt (its output on Tejus's Mac), METHODS_CHANGES.md
-- `tests/`: pytest suite (136 tests); `conftest.py` puts `AutoDeep/scripts` on the import path. `tests/data/mirgenedb_hsa_ten.csv`: ten MirGeneDB 3.0 human precursors with ViennaRNA 2.7.0 folds, chosen in advance for the B3 overhang test
+- `review/`: REVISIONS.md, reanalysis.py (reviewer's evaluation script), reanalysis_baseline.txt (its output on commit 133bf14), reanalysis_osx-arm64.txt (its output on Tejus's Mac), METHODS_CHANGES.md, homogeneity_measures.py (B2's comparison of 5' homogeneity measures on the testing folder; output in homogeneity_measures.txt)
+- `tests/`: pytest suite (139 tests); `conftest.py` puts `AutoDeep/scripts` on the import path. `tests/data/mirgenedb_hsa_ten.csv`: ten MirGeneDB 3.0 human precursors with ViennaRNA 2.7.0 folds, chosen in advance for the B3 overhang test
 - `examples/mirdeep2_tutorial/`: miRDeep2's tutorial output (`mirdeep2_output/`) and AutoDeep's expected `xgboost_judgement.csv`; checked by `tests/test_example.py`. Regenerate the expected file when a Phase B item changes a feature on purpose
 - `.github/workflows/tests.yml`: CI, pytest on ubuntu-latest and macos-latest in an env built from `AutoDeep.yml`, on pushes and PRs to `refactoring` and `main`
 

examples/mirdeep2_tutorial/README.md

@@ -7,7 +7,7 @@ A small, public miRDeep2 run to try AutoDeep on, and the output AutoDeep is expe
 - `mirdeep2_output/`: the two files AutoDeep reads from a miRDeep2 run, `result_16_10_2026_t_19_57_25.csv` and `output.mrd` (from `mirdeep_runs/run_16_10_2026_t_19_57_25/`), unchanged.
 - `expected/xgboost_judgement.csv`: `AutoDeep infer` on that folder with the bundled model (`AutoDeep/model_weights/miRNA_model.json`).
 
-The tutorial yields one novel candidate, `chrII:11534525-11540624_7`, which AutoDeep labels Confident (probabilities: Candidate 0.400, Confident 0.579, false positive 0.021). miRDeep2 also detects 6 known C. elegans miRNAs; AutoDeep only classifies novel candidates.
+The tutorial yields one novel candidate, `chrII:11534525-11540624_7`, which AutoDeep labels Confident (probabilities: Candidate 0.392, Confident 0.587, false positive 0.021). miRDeep2 also detects 6 known C. elegans miRNAs; AutoDeep only classifies novel candidates.
 
 ## Run AutoDeep on it
 

examples/mirdeep2_tutorial/expected/xgboost_judgement.csv

@@ -1,2 +1,2 @@
 provisional_id,five_prime_process_signal_selected_row_count_sum_over_mature_count_sum,num_locis,Folding Energy,3' Overhang Bot,3' Overhang Top,miRDeep2_score,mirDeep_sensitivity,sensitivity_uncertainty,total_read_count,mature_read_count,loop_read_count,star_read_count,significant_randfold_p-value,homologous_seed_in_miRBase,consensus_mature_sequence_length,mature_5'u_or_3'u,consensus_star_sequence_length,consensus_precursor_sequence_length,prediction,Candidate_Prob,Confident_Prob,falsepositive_Prob
-chrII:11534525-11540624_7,0.9876044811497292,397,-24.7,2,2,102167.2,91.0,29.0,200394,200381,0,13,yes,False,22,True,23,61,Confident,0.40007347,0.57929325,0.020633273
+chrII:11534525-11540624_7,0.980503003789957,397,-24.7,2,2,102167.2,91.0,29.0,200394,200381,0,13,yes,False,22,True,23,61,Confident,0.39205605,0.5873195,0.020624492

review/METHODS_CHANGES.md

@@ -2,6 +2,24 @@
 
 Sentences in the manuscript's Methods that must change because the code changed. One entry per change: item ID, old text, new text.
 
+## B2. 5' homogeneity: unbiased Simpson index of the mature-arm 5' starts
+
+- **Old text** (submitted manuscript, Methods, "5' Homogeneity"): "To do this, we considered only the top n sequenced permutations of a potential pre-miRNA mature strand (where n is proportional to the amount of permutations found), then we counted reads associated with the permutation with the highest read count and found the fraction over the chosen permutations that aligned in both position and base pair. The additional feature num loci tabulates the amount of permutations found."
+- **New text:** "To do this, AutoDeep takes the reads that miRDeep2 assigns to the mature arm, those in its mature read count (reads that start at most 2 nt before and end at most 5 nt after the mature sequence), and computes the unbiased Simpson index of their 5' start positions, sum n_i(n_i - 1) / [N(N - 1)], where n_i is the number of reads starting at position i and N the number of mature-arm reads. This is the probability that two reads drawn at random from the mature arm start at the same nucleotide; it is 1 when all of them do. Reads are grouped by start position whatever their first nucleotide. Candidates with fewer than two mature-arm reads have no value. A separate feature counts the distinct reads aligned to the precursor."
+- **Note:** before B2, the code divided the reads that start where the most abundant read starts, with the same first nucleotide, by every read that starts no more than 20 nt after it. When the mature is on the 3' arm, the star and loop reads start upstream and so entered the denominator: in the training data the most abundant read is on the 3' arm in 1,585 of 2,907 loci, and for those the median value was 0.929, against 0.985 for the same share counted on the mature arm only. The "top n permutations" the old text describes were computed but never used (the pypdf code that computed them went in B1). `num_locis` is unchanged. Neither text matched the code; the new text does.
+- **Note:** `training_data_pre/` keeps the old values until B8 rebuilds it, and the shipped model was trained on them. Until D1 retrains it, the model sees values on a different scale: the Simpson index is lower than the dominant read's share for the same locus (about its square: a share of 0.95 gives about 0.90). With only this feature changed, the shipped model changes 88 of the 2,907 training labels (47 Candidate to falsepositive, 25 falsepositive to Candidate). The Results sentence reading a threshold off the trees ("if the 5' ends of reads align with each other less than 91% of the time...") refers to the old definition; take any such threshold from the retrained model (C4).
+
+### For the paper (why the Simpson index; Tejus chose it on 18 Oct 2026)
+
+`review/homogeneity_measures.py` computes all of the numbers below on the mature-arm reads of the 2,907 training loci (output in `review/homogeneity_measures.txt`).
+- **Same information as the alternatives:** the Simpson index, the dominant read's share, Shannon entropy (with or without the Miller-Madow correction) and the mean distance from the commonest start rank the loci almost identically (Spearman >= 0.99). 1,105 loci (38%) have a single 5' start, which all of them score as perfectly homogeneous. Tree models use only the order of values, so the choice changes robustness and readability, not what the model can learn.
+- **No depth bias:** when the reads are thinned at random to 10%, as with a tenfold smaller library, the Simpson index shifts by +0.005 SD on average (+0.018 SD for the 903 loci with fewer than 50 mature-arm reads). The dominant share shifts by +0.053 SD (+0.116), Shannon entropy by -0.109 SD (-0.187). The unbiased estimator corrects exactly for drawing a finite number of reads. This matters for comparing species and libraries of different depth (B7, C1).
+- **Easy to state:** "the probability that two mature-arm reads start at the same nucleotide". Its inverse is the effective number of 5' ends (e.g. 1.1).
+- **Why not normalised entropy:** entropy divided by log2 of the number of distinct starts measures how evenly the reads are spread over the starts that were observed, not how concentrated they are. Two clean starts 1 nt apart score 1, the same as reads scattered over 20 starts, and adding stray reads lowers the value (toy cases in the output).
+- **Why the start position only:** Drosha and Dicer read sequence near their cleavage sites (Drosha: the basal UG and mismatched GHG motifs and the CNNC motif; Dicer: the GYM motif at the cleavage site), so a base change can shift where they cut. A shifted cut shows up as a different start position, which the index counts. A read with a different first base at the same position was still cut there; the first nucleotide matters afterwards, for Argonaute loading. Only 0.39% of mature-arm reads have a mismatched first base. Counting the first base as part of the start would change the index by more than 0.05 in 26 of 2,907 loci (Spearman 0.986). Check the references before citing: Auyeung et al. 2013 (Cell 152:844); Fang and Bartel 2015 (Mol Cell 60:131); Kwon et al. 2019 (Mol Cell 73:505); Lee, Kim and Kim 2023 (Nature 615:323).
+- **Model:** on a scratch copy of `training_data_pre/` with only this feature recomputed, the reviewer's `reanalysis.py` gives XGBoost (all 18 features) an accuracy of 0.864 in random 80/20 splits (unchanged) and 0.818 leave-one-species-out (0.812 before), with the largest per-species changes in B. tabaci (+0.03) and P. xylostella (-0.02). These differences are within the noise of a single run; the gain from B2 is a definition that is correct for 3'-arm miRNAs and robust to depth, not accuracy. The feature's share of the model's gain rises slightly (0.090 to 0.095).
+- **Caveat to state:** reads from a paralogue that align to the precursor with a mismatch at their first base count as starting at the same position. With few reads the index is unbiased but noisy; with fewer than two it is undefined (1 training locus, chrX_23311).
+
 ## B3. Signed 3' overhangs at both ends of the duplex
 
 - **Old text:** TODO (Tejus): paste the sentences that define "3' Overhang Bot" and "3' Overhang Top" (Methods and Table 1).