From 87109a2a47f5df473214cf6b1b3718921629d13d Mon Sep 17 00:00:00 2001 From: TejusK123 Date: Sun, 18 Oct 2026 23:03:43 -0500 Subject: [PATCH 1/3] B2: tests for 5' homogeneity on the mature arm (fail before the fix) Star and loop reads must not change 5' homogeneity of a 3'-arm miRNA; the value is the unbiased Simpson index of mature-arm 5' starts, whatever the first base, blank with fewer than 2 mature-arm reads. Mark B2 in progress. Co-Authored-By: Claude Opus 5.5 --- review/REVISIONS.md | 3 +- tests/test_five_prime_pipeline.py | 63 +++++++++++++++++++++++++++++++ tests/test_mrd.py | 50 +++++++++++++++++------- 3 files changed, 102 insertions(+), 14 deletions(-) create mode 100644 tests/test_five_prime_pipeline.py diff --git a/review/REVISIONS.md b/review/REVISIONS.md index 2e69baf..b5dbb5d 100644 --- a/review/REVISIONS.md +++ b/review/REVISIONS.md @@ -71,8 +71,9 @@ Do B1 first and B8 last. - [x] **B1. Parse `output.mrd`, not PDFs.** `folding_properties.py` L187-240 extract read alignments from `pdfs_*/.pdf` with pypdf, which breaks across miRDeep2 and PDF-library versions. Fix: a parser for miRDeep2's `output.mrd` (per-locus `exp` string marking flank, mature, loop and star; read lines with counts). Test: on miRDeep2's tutorial run, per-arm read counts from the parser equal the mature, star and loop counts in `result.csv`. -- [ ] **B2. 5'-homogeneity denominator (Q2).** `folding_properties.py` L231, `tdf = selection_df[(selection_df['Position'] - max_position) <= 20]`, keeps every read that starts upstream of the dominant read. For 3p-arm miRNAs, star and loop reads therefore enter the denominator and bias the feature low. Also, L230's `top_loci` (sized by the PDF page count) feeds a feature that is never saved (L244), although the manuscript describes it. +- [~] **B2. 5'-homogeneity denominator (Q2).** `folding_properties.py` L231, `tdf = selection_df[(selection_df['Position'] - max_position) <= 20]`, keeps every read that starts upstream of the dominant read. For 3p-arm miRNAs, star and loop reads therefore enter the denominator and bias the feature low. Also, L230's `top_loci` (sized by the PDF page count) feeds a feature that is never saved (L244), although the manuscript describes it. Test: adding star-arm reads to a synthetic 3p-dominant locus leaves homogeneity unchanged. + Decision (Tejus, 18 Oct 2026; Q2): the unbiased Simpson index of the 5' start positions of the mature-arm reads, grouped by position whatever the first nucleotide; blank with fewer than 2 mature-arm reads. `num_locis` is unchanged. The unsaved `top_loci` went with the pypdf code in B1. The comparison of measures is `review/homogeneity_measures.py`; the reasons are in `review/METHODS_CHANGES.md` (B2). - [x] **B3. Overhang features (Q3).** L41-43 compute "Bot" as `abs(right - left)`, the difference between the unpaired 5' and 3' ends of the whole fold; L69-154 compute "Top" with a padding heuristic. Seven S. frugiperda training rows have negative "Bot" values, so those training features came from different code. Fix: signed 3' overhangs at each duplex end from mature and star coordinates and the base-pair table, as miRDeep2 does (`miRDeep2_core_algorithm.pl`, `star_beg_exp` and `star_end_exp`). Decision (Tejus, 18 Oct 2026): the duplex rule instead of miRDeep2's expected-star rule. At the Drosha end take the outermost base pair between the mature and the star, at the Dicer end the innermost, and count the unpaired nucleotides beyond it on each strand (3'-end strand minus 5'-end strand). The two rules differ on 10 of 148 Pomonella loci with star reads; reasons in `review/METHODS_CHANGES.md` (B3). diff --git a/tests/test_five_prime_pipeline.py b/tests/test_five_prime_pipeline.py new file mode 100644 index 0000000..aca3cdb --- /dev/null +++ b/tests/test_five_prime_pipeline.py @@ -0,0 +1,63 @@ +"""B2: 5' homogeneity counts only the reads on the mature arm. + +Before B2 the denominator held every read starting no more than 20 nt after the +most abundant read, so for a miRNA on the 3' arm the star and loop reads, which +start upstream, lowered it. Runs folding_properties.py on two synthetic loci +with the same 3'-arm mature reads, one of them also with star and loop reads. +""" +import pandas as pd +import pytest + +from test_overhangs_pipeline import COMPUTED, HAIRPIN, TRAINING_CSV, hairpin_sequence, run, substructure + +# Star on the 5' arm, mature on the 3' arm; 1-based positions on HAIRPIN +STAR, MATURE = (13, 34), (51, 72) +LOOP = (STAR[1] + 1, MATURE[0] - 1) +# (start, end, reads): 90 mature reads start at 51 (two 3' ends), 10 start at 52 +MATURE_READS = [(51, 72, 60), (51, 71, 30), (52, 72, 10)] +STAR_AND_LOOP_READS = [(13, 34, 20), (14, 34, 5), (38, 45, 4)] # fewer than the mature, as miRDeep2 calls it +# Probability that two mature-arm reads, drawn without replacement, start at the same position +SIMPSON = (90 * 89 + 10 * 9) / (100 * 99) +HOMOGENEITY = "five_prime_process_signal_selected_row_count_sum_over_mature_count_sum" + + +def exp_line(): + line = ["f"] * len(HAIRPIN) + for (beg, end), mark in ((STAR, "S"), (LOOP, "l"), (MATURE, "M")): + line[beg - 1:end] = mark * (end - beg + 1) + return "".join(line) + + +def write_inputs(tmp_path): + seq = hairpin_sequence() + run_dir, out_dir = tmp_path / "mirdeep2_run", tmp_path / "AutoDeepRun" + run_dir.mkdir() + out_dir.mkdir() + loci = {"syn_3p": MATURE_READS, "syn_3p_with_star_and_loop": MATURE_READS + STAR_AND_LOOP_READS} + mrd, rnafold = [], [] + for pid, reads in loci.items(): + mrd += [f">{pid}", f"exp {exp_line()}"] + if reads is not MATURE_READS: + mrd.append(f"obs {exp_line()}") + mrd += [f"pri_seq {seq}", f"pri_struct {HAIRPIN}\t#MM"] + for n, (beg, end, count) in enumerate(reads): + mrd.append(f"seq_{n}_x{count} {'.' * (beg - 1) + seq[beg - 1:end] + '.' * (len(seq) - end)}\t0") + mrd.append("") + rnafold += [f">{pid}", seq[STAR[0] - 1:MATURE[1]], f"{substructure(STAR[0], MATURE[1])} (-30.00)"] + (run_dir / "output.mrd").write_text("\n".join(mrd) + "\n") + (out_dir / "RNAfold_novel_precursor_miRNAs.txt").write_text("\n".join(rnafold) + "\n") + # The other features: two real loci's values, which 5' homogeneity does not depend on + features = pd.read_csv(TRAINING_CSV).iloc[:len(loci), :-1].drop(columns=COMPUTED) + features["provisional_id"] = list(loci) + features.to_csv(out_dir / "feature_engineered_miRNA_Deep_data_novel_miRNAs.csv", index=False) + return run_dir + + +def test_star_and_loop_reads_leave_homogeneity_unchanged(tmp_path): + run_dir = write_inputs(tmp_path) + run("folding_properties.py", str(run_dir), cwd=tmp_path) + out = pd.read_csv(tmp_path / "AutoDeepRun" / "fully_formatted_data.csv").set_index("provisional_id") + assert out.loc["syn_3p_with_star_and_loop", HOMOGENEITY] == pytest.approx(out.loc["syn_3p", HOMOGENEITY]) + assert out.loc["syn_3p", HOMOGENEITY] == pytest.approx(SIMPSON) + # The number of distinct reads still counts every read on the precursor + assert out.loc["syn_3p_with_star_and_loop", "num_locis"] == len(MATURE_READS + STAR_AND_LOOP_READS) diff --git a/tests/test_mrd.py b/tests/test_mrd.py index dbbeccd..258d10f 100644 --- a/tests/test_mrd.py +++ b/tests/test_mrd.py @@ -1,4 +1,9 @@ -"""B1: 5' homogeneity, num_locis and per-arm read counts from miRDeep2's output.mrd.""" +"""B1: 5' homogeneity, num_locis and per-arm read counts from miRDeep2's output.mrd. + +B2 redefined 5' homogeneity as the unbiased Simpson index of the mature-arm reads' 5' starts. +""" +import math + import pytest from mrd import Locus, arm_counts, five_prime_features, parse_mrd @@ -72,20 +77,19 @@ def test_parse_rejects_ids_equal_without_commas(tmp_path): parse_mrd(str(path)) -def test_homogeneity_and_num_locis(loci): - homogeneity, num_locis = five_prime_features(loci["locus_1"].reads) - # Dominant read starts 'a' at 4: reads 1 and 2 (15). Denominator: every read - # starting at most 20 nt after it, so all but read 6 (10+5+3+2+4 = 24). - assert homogeneity == pytest.approx(15 / 24) - assert num_locis == 6 +def test_homogeneity_counts_starts_whatever_the_first_base(loci): + # All three reads are on the mature arm (5-8, allowing starts from 3). Read 9 + # starts at 3 like read 7 but with a mismatched first base: the same start. + homogeneity, num_locis = five_prime_features(loci["locus_2"]) + assert homogeneity == pytest.approx((7 * 6 + 6 * 5) / (13 * 12)) + assert num_locis == 3 -def test_tie_goes_to_first_read_in_file(loci): - # Reads 7 and 8 both have 6 reads; read 7 comes first, and read 9 starts - # at the same place but with a different (mismatched) nucleotide. - homogeneity, num_locis = five_prime_features(loci["locus_2"].reads) - assert homogeneity == pytest.approx(6 / 13) - assert num_locis == 3 +def test_homogeneity_blank_without_two_mature_reads(loci): + # locus_1's reads are all on the star, the loop or past the precursor + homogeneity, num_locis = five_prime_features(loci["locus_1"]) + assert math.isnan(homogeneity) + assert num_locis == 6 # 5' flank 1-5, mature 6-13, loop 14-25, star 26-33, 3' flank 34-40. With @@ -132,3 +136,23 @@ def test_failed_structure_has_no_arms(): def test_missing_exp_raises(): with pytest.raises(ValueError): arm_counts(Locus(None, [read(6, 13, 10)])) + with pytest.raises(ValueError): + five_prime_features(Locus(None, [read(6, 13, 10)])) + + +def test_homogeneity_is_simpson_index_of_mature_arm_starts(): + # Mature arm: starts 6 (8 reads, two different 3' ends) and 7 (2 reads). The + # read starting at 3 is outside miRDeep2's mature window; star and loop reads + # do not count. Unbiased Simpson index: sum n(n - 1) / (N(N - 1)). + locus = Locus(EXP, [read(6, 13, 6), read(6, 12, 2), read(7, 13, 2), read(3, 13, 100), + read(26, 33, 50), read(16, 23, 9)]) + homogeneity, num_locis = five_prime_features(locus) + assert homogeneity == pytest.approx((8 * 7 + 2 * 1) / (10 * 9)) + assert num_locis == 6 + + +def test_homogeneity_of_one_mature_read_is_blank(): + homogeneity, _ = five_prime_features(Locus(EXP, [read(6, 13, 1), read(26, 33, 50)])) + assert math.isnan(homogeneity) + homogeneity, _ = five_prime_features(Locus(EXP, [read(6, 13, 2)])) + assert homogeneity == 1 From fdd031a0941d80310788de851638d05d29ef774b Mon Sep 17 00:00:00 2001 From: TejusK123 Date: Sun, 18 Oct 2026 23:03:43 -0500 Subject: [PATCH 2/3] B2: 5' homogeneity as the unbiased Simpson index of mature-arm 5' starts reads_by_arm() assigns reads to arms by miRDeep2's rule for both arm_counts() and five_prime_features(). The Pomonella and FAW tests compare homogeneity with output.mrd until B8 rebuilds the training data; the tutorial's expected output is regenerated (homogeneity 0.988 -> 0.981, label unchanged). Co-Authored-By: Claude Opus 5.5 --- AutoDeep/scripts/folding_properties.py | 4 +- AutoDeep/scripts/mrd.py | 64 +++++++++++-------- examples/mirdeep2_tutorial/README.md | 2 +- .../expected/xgboost_judgement.csv | 2 +- tests/test_ids_with_commas.py | 17 +++-- tests/test_mrd_parity.py | 35 +++++----- 6 files changed, 71 insertions(+), 53 deletions(-) diff --git a/AutoDeep/scripts/folding_properties.py b/AutoDeep/scripts/folding_properties.py index 58c7d49..065a842 100644 --- a/AutoDeep/scripts/folding_properties.py +++ b/AutoDeep/scripts/folding_properties.py @@ -352,7 +352,7 @@ 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 @@ -360,6 +360,8 @@ def calculate_homogeneity_score(five_prime_df, miRNA_id=None): 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 diff --git a/AutoDeep/scripts/mrd.py b/AutoDeep/scripts/mrd.py index a6c13e8..d787dbd 100644 --- a/AutoDeep/scripts/mrd.py +++ b/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): diff --git a/examples/mirdeep2_tutorial/README.md b/examples/mirdeep2_tutorial/README.md index 9af14f2..ed1148d 100644 --- a/examples/mirdeep2_tutorial/README.md +++ b/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 diff --git a/examples/mirdeep2_tutorial/expected/xgboost_judgement.csv b/examples/mirdeep2_tutorial/expected/xgboost_judgement.csv index 19f93f6..1563d43 100644 --- a/examples/mirdeep2_tutorial/expected/xgboost_judgement.csv +++ b/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 diff --git a/tests/test_ids_with_commas.py b/tests/test_ids_with_commas.py index b6d8f57..f96f2c8 100644 --- a/tests/test_ids_with_commas.py +++ b/tests/test_ids_with_commas.py @@ -3,7 +3,9 @@ The S. frugiperda (FAW) run's contig names contain commas. init_data.py removes them from result_*.csv (this copy of the run was already rewritten that way), while output.mrd keeps them. The pipeline must still find every locus and give -the 5' features of the training rows made from this run. +the number of distinct reads of the training rows made from this run, and the +5' homogeneity computed from output.mrd (B2 redefined it; the training rows +keep the old values until B8). Needs the local run, so it is skipped unless AUTODEEP_TEST_RUN points at a directory holding mirDeep_FAW/ (result_*.csv, output.mrd). @@ -17,6 +19,8 @@ import pandas as pd import pytest +from mrd import five_prime_features, parse_mrd + REPO = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) PKG = os.path.join(REPO, "AutoDeep") TRAINING_CSV = os.path.join(PKG, "training_data_pre", "original_data_FAW.csv") @@ -25,10 +29,10 @@ pytestmark = pytest.mark.skipif(not RUN, reason="set AUTODEEP_TEST_RUN to the local testing folder") -FIVE_PRIME = ["five_prime_process_signal_selected_row_count_sum_over_mature_count_sum", "num_locis"] +HOMOGENEITY = "five_prime_process_signal_selected_row_count_sum_over_mature_count_sum" -def test_faw_run_gives_training_5p_features(tmp_path): +def test_faw_run_gives_every_locus_and_its_5p_features(tmp_path): # init_data.py rewrites result_*.csv in place, so run on a copy src = os.path.join(RUN, "mirDeep_FAW") run_dir = tmp_path / "mirDeep_FAW" @@ -47,5 +51,8 @@ def test_faw_run_gives_training_5p_features(tmp_path): training = pd.read_csv(TRAINING_CSV).set_index("provisional_id") assert len(training) == 145 assert sorted(out.index) == sorted(training.index) - pd.testing.assert_frame_equal(out.loc[training.index, FIVE_PRIME], training[FIVE_PRIME], - check_dtype=False, check_exact=False, atol=1e-9, rtol=0) + assert (out.loc[training.index, "num_locis"] == training["num_locis"]).all() + loci = parse_mrd(os.path.join(src, "output.mrd")) + expected = pd.Series({pid: five_prime_features(loci[pid])[0] for pid in training.index}) + pd.testing.assert_series_equal(out.loc[training.index, HOMOGENEITY], expected, check_names=False, + check_exact=False, atol=1e-9, rtol=0) diff --git a/tests/test_mrd_parity.py b/tests/test_mrd_parity.py index afca296..bb6049f 100644 --- a/tests/test_mrd_parity.py +++ b/tests/test_mrd_parity.py @@ -1,7 +1,8 @@ """B1: features from output.mrd match the pypdf-era training data on a real run. -B3 redefined the 3' overhangs, so those two features differ from the training -data until B8 rebuilds it; the pipeline must write the values from output.mrd. +B2 redefined 5' homogeneity and B3 the 3' overhangs, so those three features +differ from the training data until B8 rebuilds it; the pipeline must write the +values from output.mrd. Needs a local miRDeep2 run, so it is skipped unless AUTODEEP_TEST_RUN points at a directory holding mirdeep_Pomonella/ (result_*.csv, output.mrd). @@ -26,13 +27,14 @@ HOMOGENEITY = "five_prime_process_signal_selected_row_count_sum_over_mature_count_sum" OVERHANGS = ["3' Overhang Bot", "3' Overhang Top"] +REDEFINED = [HOMOGENEITY] + OVERHANGS @pytest.fixture(scope="module") -def mrd_overhangs(): - """{id: (Drosha end, Dicer end)} from the run's output.mrd (B3).""" +def mrd_values(): + """{id: (5' homogeneity, Drosha end, Dicer end)} from the run's output.mrd (B2, B3).""" loci = parse_mrd(os.path.join(RUN, "mirdeep_Pomonella", "output.mrd")) - return {pid: overhangs(locus) for pid, locus in loci.items()} + return {pid: (five_prime_features(locus)[0],) + overhangs(locus) for pid, locus in loci.items()} @pytest.fixture(scope="module") @@ -40,12 +42,10 @@ def training(): return pd.read_csv(TRAINING_CSV).set_index("provisional_id") -def test_mrd_features_match_training_data(training): +def test_mrd_num_locis_match_training_data(training): loci = parse_mrd(os.path.join(RUN, "mirdeep_Pomonella", "output.mrd")) for pid, row in training.iterrows(): - homogeneity, num_locis = five_prime_features(loci[pid].reads) - assert homogeneity == pytest.approx(row[HOMOGENEITY], abs=1e-9), pid - assert num_locis == row["num_locis"], pid + assert five_prime_features(loci[pid])[1] == row["num_locis"], pid @pytest.fixture(scope="module") @@ -64,7 +64,7 @@ def pipeline_run(tmp_path_factory): def test_pipeline_features_match_training_data(pipeline_run, training): model = xgboost.XGBClassifier() model.load_model(os.path.join(PKG, "model_weights", "miRNA_model.json")) - names = [name for name in model.get_booster().feature_names if name not in OVERHANGS] + names = [name for name in model.get_booster().feature_names if name not in REDEFINED] out = pd.read_csv(os.path.join(pipeline_run, "fully_formatted_data.csv")).set_index("provisional_id") assert set(names) <= set(out.columns) @@ -73,24 +73,25 @@ def test_pipeline_features_match_training_data(pipeline_run, training): check_dtype=False, check_exact=False, atol=1e-9, rtol=0) -def test_pipeline_overhangs_come_from_output_mrd(pipeline_run, training, mrd_overhangs): +def test_pipeline_redefined_features_come_from_output_mrd(pipeline_run, training, mrd_values): out = pd.read_csv(os.path.join(pipeline_run, "fully_formatted_data.csv")).set_index("provisional_id") - expected = pd.DataFrame([mrd_overhangs[pid] for pid in training.index], index=training.index, columns=OVERHANGS) - # Blank (no star reads) must stay blank - pd.testing.assert_frame_equal(out.loc[training.index, OVERHANGS], expected, check_dtype=False) + expected = pd.DataFrame([mrd_values[pid] for pid in training.index], index=training.index, columns=REDEFINED) + # Blank (no star reads, or fewer than 2 mature-arm reads) must stay blank + pd.testing.assert_frame_equal(out.loc[training.index, REDEFINED], expected, check_dtype=False, + check_exact=False, atol=1e-9, rtol=0) -def test_pipeline_predictions_match_training_features(pipeline_run, training, mrd_overhangs): +def test_pipeline_predictions_match_training_features(pipeline_run, training, mrd_values): # The run's stored xgboost_judgement.csv came from an older model, so compare # with the shipped model's predictions on the pypdf-era training features, - # with the overhangs from output.mrd (B3) + # with 5' homogeneity and the overhangs from output.mrd (B2, B3) from features import decode_labels, encode_features model = xgboost.XGBClassifier() model.load_model(os.path.join(PKG, "model_weights", "miRNA_model.json")) names = model.get_booster().feature_names features = training.copy() - features[OVERHANGS] = [mrd_overhangs[pid] for pid in features.index] + features[REDEFINED] = [mrd_values[pid] for pid in features.index] expected = decode_labels(model.predict(encode_features(features)[names].to_numpy())) new = pd.read_csv(os.path.join(pipeline_run, "xgboost_judgement.csv")).set_index("provisional_id") From d165070d5b941a023de3363cd27242bd0582f995 Mon Sep 17 00:00:00 2001 From: TejusK123 Date: Sun, 18 Oct 2026 23:03:43 -0500 Subject: [PATCH 3/3] B2: Methods text, the comparison of homogeneity measures, CLAUDE.md notes Co-Authored-By: Claude Opus 5.5 --- CLAUDE.md | 8 +- review/METHODS_CHANGES.md | 18 +++++ review/homogeneity_measures.py | 136 ++++++++++++++++++++++++++++++++ review/homogeneity_measures.txt | 43 ++++++++++ 4 files changed, 201 insertions(+), 4 deletions(-) create mode 100644 review/homogeneity_measures.py create mode 100644 review/homogeneity_measures.txt diff --git a/CLAUDE.md b/CLAUDE.md index 48b5659..066f5b9 100644 --- a/CLAUDE.md +++ b/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//`): 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` diff --git a/review/METHODS_CHANGES.md b/review/METHODS_CHANGES.md index fe7590f..3acf4c8 100644 --- a/review/METHODS_CHANGES.md +++ b/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). diff --git a/review/homogeneity_measures.py b/review/homogeneity_measures.py new file mode 100644 index 0000000..14f4308 --- /dev/null +++ b/review/homogeneity_measures.py @@ -0,0 +1,136 @@ +"""Compare ways of measuring 5' homogeneity on the training loci (B2, Q2). + +Usage: python review/homogeneity_measures.py > review/homogeneity_measures.txt + + holds the miRDeep2 run behind each training CSV (the local +testing folder in CLAUDE.local.md). Every measure is computed on the reads that +miRDeep2 assigns to the mature arm. Reported: toy cases, agreement between the +measures, how well each separates the curated labels on its own, how much each +shifts when the reads are thinned at random (as with lower sequencing depth), +and how much counting the first nucleotide as part of the start changes the +Simpson index. B2 uses the unbiased Simpson index of the start positions. +""" +import glob +import os +import sys + +import numpy as np +import pandas as pd + +REPO = os.path.dirname(os.path.dirname(os.path.abspath(__file__))) +sys.path.insert(0, os.path.join(REPO, 'AutoDeep', 'scripts')) +from mrd import parse_mrd, reads_by_arm # noqa: E402 + +RUNS = {'Andreii': 'mirDeep_Andreii', 'Arabidopsis': 'mirDeep_Arabidopsis', 'Bemisia': 'mirDeep_Bemisicaa', + 'Drosophila': 'mirDeep_Drosophila', 'FAW': 'mirDeep_FAW', 'Human': 'mirDeep_Human', + 'Pomonella': 'mirdeep_Pomonella', 'Xylostella': 'mirdeep_Xylostella'} +LABEL = 'Candidate/Confident/no-mirna' +MEASURES = ['simpson', 'p_max', 'within1', 'H_bits', 'H_mm', 'H_norm', 'mad_nt'] +HIGHER_IS_HOMOGENEOUS = {'simpson': 1, 'p_max': 1, 'within1': 1, 'H_bits': -1, 'H_mm': -1, 'H_norm': -1, 'mad_nt': -1} +SEED = 42 + + +def simpson(groups): + """Unbiased Simpson index: the probability that two reads, drawn without replacement, fall in the same group.""" + n = np.array(list(groups.values()), dtype=float) + N = n.sum() + return (n * (n - 1)).sum() / (N * (N - 1)) if N > 1 else np.nan + + +def measures(starts): + """starts: {5' start: reads}.""" + pos = np.array(sorted(starts)) + n = np.array([starts[p] for p in pos], dtype=float) + N, K = n.sum(), len(n) + if N == 0: + return {k: np.nan for k in MEASURES} + p = n / N + H = -(p * np.log2(p)).sum() + 0.0 # + 0.0 turns -0.0 into 0.0 + mode = pos[np.argmax(n)] + return { + 'simpson': simpson(starts), # P(two reads share a start) + 'p_max': p.max(), # share at the commonest start + 'within1': p[np.abs(pos - mode) <= 1].sum(), # share within 1 nt of it + 'H_bits': H, # Shannon entropy + 'H_mm': H + (K - 1) / (2 * N * np.log(2)), # with Miller-Madow correction + 'H_norm': H / np.log2(K) if K > 1 else 0.0, # entropy / log2(distinct starts) + 'mad_nt': (p * np.abs(pos - mode)).sum(), # mean distance from the commonest start + } + + +def group(keys, counts): + out = {} + for k, c in zip(keys, counts): + if c: + out[k] = out.get(k, 0) + c + return out + + +def auc(pos, neg): + pos, neg = np.asarray(pos.dropna()), np.asarray(neg.dropna()) + return ((pos[:, None] > neg[None, :]).sum() + 0.5 * (pos[:, None] == neg[None, :]).sum()) / (len(pos) * len(neg)) + + +def main(testing): + print("Toy cases (reads at each 5' start):") + toys = {'1000 at one start': {0: 1000}, '1000 + 10 stray reads': {0: 1000, **{i: 1 for i in range(3, 13)}}, + '500 + 500, 1 nt apart': {0: 500, 1: 500}, 'the same + 50 stray': {0: 500, 1: 500, 9: 50}, + '90% + 10% at +1 nt': {0: 900, 1: 100}, '90% + 10% scattered': {0: 900, **{i: 10 for i in range(3, 13)}}, + 'scattered over 20 starts': {i: 10 for i in range(20)}} + print(pd.DataFrame({k: measures(v) for k, v in toys.items()}).T[MEASURES].round(3).to_string()) + + rng = np.random.default_rng(SEED) + recs = [] + for species, run in RUNS.items(): + loci = {pid.replace(',', ''): locus for pid, locus in parse_mrd(os.path.join(testing, run, 'output.mrd')).items()} + for f in sorted(glob.glob(os.path.join(REPO, 'AutoDeep', 'training_data_pre', f'original_data_{species}*.csv'))): + for pid, label in pd.read_csv(f)[['provisional_id', LABEL]].itertuples(index=False): + locus = loci[pid] + mature = [(beg, aln[beg - 1], count) for (arm, beg, count), (aln, _) in zip(reads_by_arm(locus), locus.reads) + if arm == 'mature'] + begs = [b for b, _, _ in mature] + counts = np.array([c for _, _, c in mature]) + rec = {'species': species, 'label': label, 'N': counts.sum(), **measures(group(begs, counts))} + rec['simpson_with_first_base'] = simpson(group([(b, nt.lower()) for b, nt, _ in mature], counts)) + rec['first_base_mismatched'] = sum(c for _, nt, c in mature if nt.isupper()) + for frac in (0.25, 0.1): + draws = [measures(group(begs, rng.binomial(counts, frac))) for _ in range(5)] + for k in MEASURES: + rec[f'{k}@{frac}'] = np.nanmean([d[k] for d in draws]) + recs.append(rec) + df = pd.DataFrame(recs) + + print(f"\n{len(df)} training loci; mature-arm reads per locus: median {df.N.median():.0f}, " + f"fewer than 50: {(df.N < 50).sum()}, fewer than 2: {(df.N < 2).sum()}; " + f"a single start: {(df.p_max == 1).sum()}") + print('\nSpearman correlation between the measures:') + print(df[MEASURES].corr(method='spearman').round(2).to_string()) + + print('\nAUC of each measure alone (oriented so that higher = more homogeneous):') + for k in MEASURES: + s = df[k] * HIGHER_IS_HOMOGENEOUS[k] + pairs = [('Confident', 'falsepositive'), ('Candidate', 'falsepositive'), ('Confident', 'Candidate')] + print(f' {k:8s} ' + ' '.join(f'{a} vs {b} {auc(s[df.label == a], s[df.label == b]):.3f}' for a, b in pairs)) + + print('\nDepth: mean shift when reads are thinned to 25% and 10% (5 draws each), in SDs of the measure ' + 'across loci, for all loci and for loci with fewer than 50 mature-arm reads; Spearman between full depth and 10%') + low = df.N < 50 + for k in MEASURES: + sd = df[k].std() + shift = {frac: (df[f'{k}@{frac}'] - df[k]) / sd for frac in (0.25, 0.1)} + rho = df[[k, f'{k}@0.1']].corr(method='spearman').iloc[0, 1] + print(f' {k:8s} 25%: {shift[0.25].mean():+.3f} 10%: {shift[0.1].mean():+.3f} ' + f'10%, <50 reads: {shift[0.1][low].mean():+.3f} Spearman {rho:.3f}') + + print("\nFirst nucleotide (output.mrd prints a mismatched base in upper case):") + print(f" mature-arm reads with a mismatched first base: {df.first_base_mismatched.sum() / df.N.sum():.2%}") + diff = df.simpson_with_first_base - df.simpson + print(f" Simpson index with the first base as part of the start, minus position only: differs by more than " + f"0.05 in {(diff.abs() > 0.05).sum()} loci; Spearman " + f"{df[['simpson', 'simpson_with_first_base']].corr(method='spearman').iloc[0, 1]:.3f}") + + +if __name__ == '__main__': + if len(sys.argv) != 2: + sys.exit(__doc__) + main(sys.argv[1]) diff --git a/review/homogeneity_measures.txt b/review/homogeneity_measures.txt new file mode 100644 index 0000000..a325a6c --- /dev/null +++ b/review/homogeneity_measures.txt @@ -0,0 +1,43 @@ +Toy cases (reads at each 5' start): + simpson p_max within1 H_bits H_mm H_norm mad_nt +1000 at one start 1.000 1.000 1.000 0.000 0.000 0.000 0.000 +1000 + 10 stray reads 0.980 0.990 0.990 0.113 0.120 0.033 0.074 +500 + 500, 1 nt apart 0.499 0.500 1.000 1.000 1.001 1.000 0.500 +the same + 50 stray 0.455 0.476 0.952 1.229 1.230 0.775 0.905 +90% + 10% at +1 nt 0.820 0.900 1.000 0.469 0.470 0.469 0.100 +90% + 10% scattered 0.811 0.900 0.900 0.801 0.808 0.232 0.750 +scattered over 20 starts 0.045 0.050 0.100 4.322 4.390 1.000 9.500 + +2907 training loci; mature-arm reads per locus: median 161, fewer than 50: 903, fewer than 2: 1; a single start: 1105 + +Spearman correlation between the measures: + simpson p_max within1 H_bits H_mm H_norm mad_nt +simpson 1.00 1.00 0.69 -1.00 -1.00 -0.99 -0.99 +p_max 1.00 1.00 0.69 -1.00 -1.00 -0.99 -0.99 +within1 0.69 0.69 1.00 -0.71 -0.71 -0.63 -0.73 +H_bits -1.00 -1.00 -0.71 1.00 1.00 0.98 0.99 +H_mm -1.00 -1.00 -0.71 1.00 1.00 0.98 0.99 +H_norm -0.99 -0.99 -0.63 0.98 0.98 1.00 0.98 +mad_nt -0.99 -0.99 -0.73 0.99 0.99 0.98 1.00 + +AUC of each measure alone (oriented so that higher = more homogeneous): + simpson Confident vs falsepositive 0.580 Candidate vs falsepositive 0.612 Confident vs Candidate 0.455 + p_max Confident vs falsepositive 0.577 Candidate vs falsepositive 0.609 Confident vs Candidate 0.455 + within1 Confident vs falsepositive 0.541 Candidate vs falsepositive 0.617 Confident vs Candidate 0.390 + H_bits Confident vs falsepositive 0.581 Candidate vs falsepositive 0.615 Confident vs Candidate 0.452 + H_mm Confident vs falsepositive 0.582 Candidate vs falsepositive 0.616 Confident vs Candidate 0.453 + H_norm Confident vs falsepositive 0.584 Candidate vs falsepositive 0.598 Confident vs Candidate 0.473 + mad_nt Confident vs falsepositive 0.585 Candidate vs falsepositive 0.618 Confident vs Candidate 0.453 + +Depth: mean shift when reads are thinned to 25% and 10% (5 draws each), in SDs of the measure across loci, for all loci and for loci with fewer than 50 mature-arm reads; Spearman between full depth and 10% + simpson 25%: -0.001 10%: +0.005 10%, <50 reads: +0.018 Spearman 0.955 + p_max 25%: +0.016 10%: +0.053 10%, <50 reads: +0.116 Spearman 0.949 + within1 25%: +0.008 10%: +0.036 10%, <50 reads: +0.091 Spearman 0.923 + H_bits 25%: -0.045 10%: -0.109 10%, <50 reads: -0.187 Spearman 0.948 + H_mm 25%: -0.027 10%: -0.079 10%, <50 reads: -0.166 Spearman 0.947 + H_norm 25%: +0.005 10%: -0.035 10%, <50 reads: -0.200 Spearman 0.934 + mad_nt 25%: -0.012 10%: -0.045 10%, <50 reads: -0.108 Spearman 0.948 + +First nucleotide (output.mrd prints a mismatched base in upper case): + mature-arm reads with a mismatched first base: 0.39% + Simpson index with the first base as part of the start, minus position only: differs by more than 0.05 in 26 loci; Spearman 0.986