diff --git a/CHANGELOG.md b/CHANGELOG.md index f7f6505..575431b 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -1,5 +1,14 @@ # Change Log +## [v10.11.6](https://github.com/openvax/varcode/tree/v10.11.6) (2026-10-05) + +- Exercise real Sid RNA evidence through `VariantCollection.effects()` for + NTF3 compound codons and MAP2 complex haplotypes (#464). Pin independent + NCBI coding/protein references, retain complete available sequencing + templates, and verify individual predictions, phase uncertainty and joint + effect provenance. Keep corrected, deletion-only and historical MAP2 + alleles distinct; RNA support does not replace their DNA predictions. + ## [v10.11.5](https://github.com/openvax/varcode/tree/v10.11.5) (2026-10-05) - Test RNA read phasing and molecular phase resolution on real Sid RNA diff --git a/tests/README.md b/tests/README.md index 618f0fd..a2979ce 100644 --- a/tests/README.md +++ b/tests/README.md @@ -227,7 +227,7 @@ cache); later runs reuse them offline. in [#528](https://github.com/openvax/varcode/issues/528). It uses the published `openvax-v2` bundle, with manifest SHA256 `b42dbce529cb5724153f035fca3783a41e673b31610a30851b2ef96d4ba43734`. -The test exports three named members into a pytest temporary directory, +The phasing checks export three named members into a pytest temporary directory, preserving their original alignment records and indexed BAM format: | Source label | Target/member suffix | Records | Templates | Supplementary records | @@ -294,3 +294,103 @@ python -m pytest -q tests/test_osteosarc_rna_phasing.py The first run can download the verified 28 MB shared bundle; later runs reuse it offline. No sibling checkout or BAM copies in `tests/data` are required. + +## Real Sid RNA protein effects + +The same module passes `MolecularPhaseResolver(RNAReadPhasingSource(...))` to +`VariantCollection.effects()` for NTF3 and complex MAP2 cases, following +[#464](https://github.com/openvax/varcode/issues/464). It checks every individual +DNA prediction against the collection without RNA phasing, including other +transcripts, and verifies the selected transcript's complete protein against +independent expectations. NTF3 runs with `fast`, `protein_diff` and the opt-in +`transcript_model`; MAP2 runs with the two established point-edit annotators. +Joint legacy effects expose their edited cDNA, and the transcript model's +realized candidates expose theirs; the tests check their CDS, affected codon +and stop codon. Joint membership, protein and phase provenance survive JSON. +The experimental model currently loses its realized candidates/hypotheses +and their cDNA on JSON round-trip; this reproduced limitation is tracked in +[#574](https://github.com/openvax/varcode/issues/574). Its full CDS/codon checks +use the in-memory products, and its round-trip checks cover the preserved +protein, membership and phase source. The audit also reproduced a support-cache +key defect ([#575](https://github.com/openvax/varcode/issues/575)); returned +counts are correct here, but repeated calls rescan the BAM. + +`tests/data/sid_rna_effect_references.json` is a separate offline reference +asset (SHA256 `d9e21cc96450d4a14e7b41d0900df48792c9f64b103b9ce01ab5cc841508d29e`). +On 2026-10-05, NCBI EFetch supplied these versioned primary records: + +| Gene | RefSeq mRNA / protein | GenBank CDS, one-based inclusive | Ensembl 81 transcript | +|---|---|---|---| +| NTF3 | [NM_002527.5](https://www.ncbi.nlm.nih.gov/nuccore/NM_002527.5) / NP_002518.1 | 84–857 | ENST00000331010 | +| MAP2 | [NM_002374.4](https://www.ncbi.nlm.nih.gov/nuccore/NM_002374.4) / NP_002365.3 | 321–5804 | ENST00000360351 | + +The acquired CDS and `/translation` protein strings match those Ensembl +transcripts exactly. Independently retrieved GRCh38 chromosome windows +`NC_000012.12:5494300–5494500` and `NC_000002.12:209694750–209694830` each match +one interval of the corresponding RefSeq CDS, at zero-based offsets 85 and +2579. The asset retains these strings, URLs and response/sequence hashes. +To repeat acquisition, fetch each recorded URL, extract its CDS interval and +`/translation`, and compare the sequence hashes. GenBank descriptive metadata +can change while the accession's sequence remains the same. + +The expected edits were applied directly to these RefSeq CDS strings and +translated independently with [NCBI's standard genetic code](https://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi#SG1), +without importing Varcode, PyEnsembl or using upstream peptide annotations. +The asset records fixed CDS edit offsets, codons, literal protein replacements, +complete mutant CDS/protein hashes, lengths and stops. Tests read those frozen +expectations; they do not compute a second prediction from Varcode. + +- NTF3's zero-based CDS offsets 166 and 167 give `AAG>AGG` (`p.K56R`) and + `AAG>AAT` (`p.K56N`) individually. Together they give `AAG>AGT` (`p.K56S`), + with local peptide `DVSENY` instead of `DVKENY`. Both short-read and ONT cis + evidence produce this joint effect while retaining both individual effects. +- The real ONT trans companion `12:5494466 G>A` changes `CCG>CCA` at proline + 84. Its silent effect and the separate `p.K56R` prediction remain present. + Increasing the compound pair's minimum to 10 leaves its nine observations + insufficient: phase is unknown, and the two individual predictions remain. + +The MAP2 reads use three additional `openvax-v2` members under +`isovar/osteosarc/figure_comparisons/corpus/context-evidence.json.gz#/map2/sources/`: + +| Member suffix | Source | Regional member records | Complete available records / templates | Q20 compound templates | +|---|---|---:|---:|---:| +| `11/original_records` | `IPISRC044_tumor_T1_ucla_rna.md` | 123 | 125 / 78 | 4 | +| `12/original_records` | `IPISRC044_tumor_T2_ucla_rna.md` | 25 | 26 / 18 | 2 | +| `3/original_records` | `IPISRC044_T2_sclrs_ONT_dedup` | 32 | 33 / 32 | 1 | + +These are original regional acquisitions retained in the neutral bundle. The +fixture selects their RG/QNAMEs, then copies **all** available records of those +templates from the bundle source into an indexed temporary BAM. No alignment +fields change; record digest multiplicities must equal the source's complete +template subset. The ONT input retains one supplementary record, which the +default evidence filter excludes. The two bulk inputs have no supplementary +records. Source products are RNA products, checked by their manifest URLs. + +The corrected GRCh38 MAP2 allele is explicitly decomposed into +`2:209694769 C>A`, `2:209694770 T>G` and deletion of +`2:209694773–209694800 GCTACTGTGTGTTCAATAAGTACACAGT`. The existing +`register_haplotype()` API tests this hypothesis's original sequence between +five-base genomic flanks: `AAGACAGGGCCCAT`. Independent aligned-position, +query-sequence and quality inspection found four T1 bulk, two T2 bulk and one +T2 ONT supporting RG/QNAME templates. The bulk witnesses encode a `28N` gap; +the ONT witness encodes `22D2M6D`. Exact context allows the same RNA sequence +to support the hypothesis without assigning a genomic deletion mechanism. +Default MAPQ/Q20, pairing, flag, five-base edge and two-fragment filters remain +in force; ONT uses the existing 5% error allowance. Missing or different context +stays unknown and supplies no reference/trans evidence. + +The bulk evidence yields a joint 904-aa protein: reference residues 1–866 +followed by `RAHCHHLFKTVRIYQGRVVPFTKALMIKFEEIWPQTFH`, then `TGA`. The single +ONT compound witness leaves phase unknown and yields no joint prediction. +Separate predictions for the two substitutions and deletion remain in all cases. +Additional checks retain the corrected complex allele, the 28-base deletion +alone (904 aa, with `LAH` instead of `RAH` at the junction), and the historical +22-base deletion (906 aa, with `STQSH` at the junction) as distinct DNA predictions. +Only the corrected compound sequence matches these RNA witnesses. Zero matches +for the other hypotheses do not remove their effects or establish their absence. + +These proteins are reference-frame predictions for observed local sequence +combinations. The tests do not reconstruct full RNA molecules, establish an +expressed isoform, validate tumor haplotypes, estimate allele fractions from the +bounded bundle, or pool processed libraries as independent molecules. No runtime +dependency, RNA assembly implementation or global phasing policy changes. diff --git a/tests/data/README.md b/tests/data/README.md index b138640..61a85a8 100644 --- a/tests/data/README.md +++ b/tests/data/README.md @@ -13,6 +13,7 @@ No BAM/CRAM files or shared RNA bundles belong here. See |---|---| | `osteosarc_variants.json`, `osteosarc_snapshot_2026-09-18t.zip` | [Snapshot identity, receipts, archive hash and offline export recipe](../README.md#regenerate-from-the-verified-snapshot). `../collect_osteosarc_variants.py` exports the reviewed snapshot with the pinned Osteosarc version. | | `osteosarc_observed_junctions.json` | Each record retains its source page, ONT read name, reference identity and source-sequence SHA256. [Shared read verification](../README.md#shared-test-data-openvax-v2) uses the external bundle; only the bounded junction windows are stored here. | +| `sid_rna_effect_references.json` | Versioned NCBI RefSeq CDS/protein records and GRCh38 genomic windows, acquired 2026-10-05 independently of Varcode/PyEnsembl. URLs, accession versions, response and sequence hashes, CDS edits and manually audited protein expectations are recorded in the file. [RNA effect verification](../README.md#real-sid-rna-protein-effects) explains the audit. | | `osteosarc_esvee_somatic.vcf`, `osteosarc_esvee_shared_svs.vcf` | Published Sid esvee records with FORMAT/sample columns removed. The subset selection and reference releases are documented in [the fusion tests](../test_osteosarc_fusions.py) and [shared SV tests](../test_osteosarc_shared_svs.py). | | `real_callers/*.vcf` | [Pinned GATK source and synthetic Strelka2/VEP construction notes](real_callers/README.md). | | `spec_examples/*.vcf` | [VCF specification sections and source](spec_examples/README.md). | @@ -25,8 +26,8 @@ No BAM/CRAM files or shared RNA bundles belong here. See The original external acquisition recipes for the legacy inputs were not recorded. Their repository history is the available provenance; no new upstream -identity or scientific validation is inferred here. All fixture payloads in -this release are byte-for-byte those in +identity or scientific validation is inferred here. The legacy fixture payloads +are byte-for-byte those in [commit 6f29cf8](https://github.com/openvax/varcode/tree/6f29cf8c5201189ae2803ad319aacca8053f7ddf/tests/data). To recover those exact payloads, use the raw files under that commit's `tests/data/` tree and verify them against `SHA256SUMS`. The new inventory notes diff --git a/tests/data/SHA256SUMS b/tests/data/SHA256SUMS index 1c9b261..e71574c 100644 --- a/tests/data/SHA256SUMS +++ b/tests/data/SHA256SUMS @@ -1,4 +1,4 @@ -b91532598e8ef0034ac5f9db42ad23b80d4881aac57e278c1c70db57c77c25f1 README.md +12f256e28c8420d5f0e6a22a89c9fb714f4f96f85c0474d6978f1df7b1ffa9f7 README.md e569c37ac74f13111be67d230f0bf03891a419af31afb04c7e3a28212cba8bac dbnsfp_validation_set.csv a3d5571678f7a4717bec59325fe77131d54d814b27f99c494d944d78a07de501 different-samples.1.vcf b2cb6090f568cfacdb0df4a9c280161923c255712fa368ac433863a623d126a5 different-samples.2.vcf @@ -26,6 +26,7 @@ eabc38358591c380ddb676e2d53a564cb01a33e428b67186a009f5b3cad4ce93 real_callers/s 91c8c7dd7e9915a73fd0e709dde8d5cc74d4f8fa3af14df0dc63dfd400771674 real_callers/vep_annotated_csq.vcf a3d5571678f7a4717bec59325fe77131d54d814b27f99c494d944d78a07de501 same-samples.1.vcf a3d5571678f7a4717bec59325fe77131d54d814b27f99c494d944d78a07de501 same-samples.2.vcf +d9e21cc96450d4a14e7b41d0900df48792c9f64b103b9ce01ab5cc841508d29e sid_rna_effect_references.json 94bf15f0cca60a2d2e14526cc1da7a3933aaedea18f57ba4ada75462ad8943b6 simple.1.vcf 6910e2b188142b0a9702aac57cc88d29987c47e9d4adc1db0bc24b3543d86b47 simple.2.vcf 5d988950e60a08f40fb406d7d40a4e2bf2c4df402c73395753db39933d8409ea somatic_hg19_14muts.space_in_sample_name.vcf diff --git a/tests/data/sid_rna_effect_references.json b/tests/data/sid_rna_effect_references.json new file mode 100644 index 0000000..ca5dbfd --- /dev/null +++ b/tests/data/sid_rna_effect_references.json @@ -0,0 +1,250 @@ +{ + "assembly": "GRCh38", + "ensembl_release": 81, + "acquired": "2026-10-05", + "translation_table": 1, + "translation_source": "https://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi#SG1", + "references": { + "NTF3": { + "transcript_id": "ENST00000331010", + "mrna_accession": "NM_002527.5", + "protein_accession": "NP_002518.1", + "genbank_url": "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=NM_002527.5&rettype=gb&retmode=text", + "genbank_response_sha256": "c45e0ca5bf63b225d512e6a987097c39c5ab42a3efe3cffcc0d62a98607eb684", + "cds_interval_1based": [ + 84, + 857 + ], + "cds": "ATGTCCATCTTGTTTTATGTGATATTTCTCGCTTATCTCCGTGGCATCCAAGGTAACAACATGGATCAAAGGAGTTTGCCAGAAGACTCGCTCAATTCCCTCATTATTAAGCTGATCCAGGCAGATATTTTGAAAAACAAGCTCTCCAAGCAGATGGTGGACGTTAAGGAAAATTACCAGAGCACCCTGCCCAAAGCTGAGGCTCCCCGAGAGCCGGAGCGGGGAGGGCCCGCCAAGTCAGCATTCCAGCCGGTGATTGCAATGGACACCGAACTGCTGCGACAACAGAGACGCTACAACTCACCGCGGGTCCTGCTGAGCGACAGCACCCCCTTGGAGCCCCCGCCCTTGTATCTCATGGAGGATTACGTGGGCAGCCCCGTGGTGGCGAACAGAACATCACGGCGGAAACGGTACGCGGAGCATAAGAGTCACCGAGGGGAGTACTCGGTATGTGACAGTGAGAGTCTGTGGGTGACCGACAAGTCATCGGCCATCGACATTCGGGGACACCAGGTCACGGTGCTGGGGGAGATCAAAACGGGCAACTCTCCCGTCAAACAATATTTTTATGAAACGCGATGTAAGGAAGCCAGGCCGGTCAAAAACGGTTGCAGGGGTATTGATGATAAACACTGGAACTCTCAGTGCAAAACATCCCAAACCTACGTCCGAGCACTGACTTCAGAGAACAATAAACTCGTGGGCTGGCGGTGGATACGGATAGACACGTCCTGTGTGTGTGCCTTGTCGAGAAAAATCGGAAGAACATGA", + "cds_sha256": "0b2e583378a782d557a331dd93a4fdfec04d36ba6c66b9b4164fd3af3d309598", + "protein": "MSILFYVIFLAYLRGIQGNNMDQRSLPEDSLNSLIIKLIQADILKNKLSKQMVDVKENYQSTLPKAEAPREPERGGPAKSAFQPVIAMDTELLRQQRRYNSPRVLLSDSTPLEPPPLYLMEDYVGSPVVANRTSRRKRYAEHKSHRGEYSVCDSESLWVTDKSSAIDIRGHQVTVLGEIKTGNSPVKQYFYETRCKEARPVKNGCRGIDDKHWNSQCKTSQTYVRALTSENNKLVGWRWIRIDTSCVCALSRKIGRT", + "protein_sha256": "51d8e5a0dcb4df0c94896dd5fab8ca5deeb4523aaee8ba6ec6948763443b91e2", + "genomic_window": { + "accession": "NC_000012.12", + "start": 5494300, + "end": 5494500, + "cds_offset": 85, + "url": "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=NC_000012.12&rettype=fasta&retmode=text&seq_start=5494300&seq_stop=5494500", + "response_sha256": "288a13644ba261b6bd8adf7f6a94875f8213bf20e8ff92ae23c2c047b79a8726", + "sequence": "ACTCGCTCAATTCCCTCATTATTAAGCTGATCCAGGCAGATATTTTGAAAAACAAGCTCTCCAAGCAGATGGTGGACGTTAAGGAAAATTACCAGAGCACCCTGCCCAAAGCTGAGGCTCCCCGAGAGCCGGAGCGGGGAGGGCCCGCCAAGTCAGCATTCCAGCCGGTGATTGCAATGGACACCGAACTGCTGCGACAAC", + "sequence_sha256": "ae0b62fd71e7116f361ff6cc80c4ad9724a4fee34ca1a683cae3e8a970667cf7" + } + }, + "MAP2": { + "transcript_id": "ENST00000360351", + "mrna_accession": "NM_002374.4", + "protein_accession": "NP_002365.3", + "genbank_url": "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=NM_002374.4&rettype=gb&retmode=text", + "genbank_response_sha256": "36bc9959cdefa8a34f946bdb6a4024237f860dd603367f9d5faa69cd98bc76a9", + "cds_interval_1based": [ + 321, + 5804 + ], + "cds": "ATGGCAGATGAACGGAAAGATGAAGCAAAGGCACCTCACTGGACCTCAGCACCGCTAACAGAGGCATCTGCACACTCACATCCACCTGAGATTAAGGATCAAGGCGGAGCAGGGGAAGGACTTGTCCGAAGCGCCAATGGATTCCCATACAGGGAGGATGAAGAGGGTGCCTTTGGAGAGCATGGGTCACAGGGCACCTATTCAAATACCAAAGAGAATGGGATCAACGGAGAGCTGACCTCAGCTGACAGAGAAACAGCAGAGGAGGTGTCTGCAAGGATAGTTCAAGTAGTCACTGCTGAGGCTGTAGCAGTCCTGAAAGGTGAACAAGAGAAAGAAGCTCAACATAAAGACCAGACTGCAGCTCTGCCTTTAGCAGCTGAAGAAACAGCTAATCTGCCTCCTTCTCCACCCCCATCACCTGCCTCAGAACAGACTGTCACAGTGGAGGAAGATTTACTTACAGCCTCGAAGATGGAGTTCCACGATCAACAGGAATTGACTCCCTCTACAGCTGAGCCTTCAGACCAGAAGGAAAAGGAGTCAGAGAAGCAAAGTAAGCCTGGTGAAGACCTTAAACATGCTGCCTTAGTTTCTCAGCCAGAGACAACTAAAACTTACCCTGATAAAAAGGACATGCAAGGCACGGAAGAAGAAAAAGCACCCCTAGCTTTGTTTGGGCACACTCTTGTTGCCAGCCTGGAAGACATGAAACAGAAGACAGAACCAAGCCTTGTAGTACCTGGCATTGACCTCCCTAAAGAGCCTCCAACTCCAAAAGAACAAAAGGACTGGTTCATCGAAATGCCAACGGAAGCAAAAAAGGATGAGTGGGGTTTAGTTGCCCCCATATCTCCTGGCCCTCTGACTCCCATGAGGGAAAAAGATGTATTTGATGATATCCCAAAATGGGAAGGGAAACAGTTTGATTCTCCCATGCCAAGTCCCTTTCAAGGGGGAAGCTTCACTCTTCCTTTAGATGTCATGAAGAATGAAATAGTTACAGAAACATCGCCCTTTGCCCCTGCCTTTTTACAGCCAGATGACAAAAAATCTCTGCAACAAACCAGTGGCCCAGCTACTGCCAAAGATAGTTTTAAAATTGAAGAGCCCCATGAGGCTAAACCTGACAAAATGGCAGAAGCACCACCCTCAGAGGCAATGACCTTACCCAAAGATGCTCACATTCCAGTTGTAGAAGAACATGTTATGGGGAAAGTTTTAGAGGAAGAAAAGGAGGCCATAAATCAAGAGACTGTGCAGCAAAGGGATACTTTCACCCCCAGTGGACAGGAACCTATACTTACTGAAAAGGAAACTGAGCTGAAGCTTGAAGAAAAAACCACCATTTCTGACAAAGAAGCTGTGCCAAAAGAGAGTAAACCCCCAAAACCTGCAGATGAAGAAATAGGCATAATTCAGACCTCCACAGAGCACACTTTCTCAGAACAGAAAGACCAAGAGCCTACCACAGATATGTTGAAACAGGACTCGTTCCCTGTAAGTTTGGAGCAAGCAGTTACAGATTCAGCCATGACCTCTAAAACACTGGAGAAAGCCATGACCGAACCATCTGCATTAATTGAAAAGAGCTCAATTCAGGAACTTTTTGAAATGAGAGTTGATGACAAAGATAAGATTGAAGGAGTTGGAGCTGCAACATCAGCTGAGCTTGATATGCCATTTTATGAAGATAAATCAGGAATGTCCAAGTACTTTGAAACATCTGCCTTGAAAGAAGAAGCAACAAAAAGCATTGAGCCAGGCAGTGATTACTATGAACTGAGTGACACTAGAGAAAGTGTCCATGAGTCTATTGATACCATGTCTCCCATGCATAAAAATGGTGACAAGGAGTTTCAAACAGGAAAAGAATCCCAGCCCAGTCCTCCAGCACAAGAAGCAGGGTACAGCACTCTCGCACAGAGTTATCCATCAGATTTACCTGAAGAACCCAGTTCTCCTCAAGAAAGAATGTTCACTATTGATCCAAAAGTGTATGGAGAGAAAAGGGACCTCCACAGTAAGAATAAGGATGATTTGACCCTTAGCAGGAGTTTAGGACTTGGTGGTAGGTCTGCAATAGAACAAAGAAGCATGTCAATCAATTTGCCGATGTCTTGCCTAGATTCCATAGCCCTTGGATTTAACTTTGGTCGGGGACATGATCTTTCTCCTCTGGCTTCCGATATTCTAACCAACACTAGTGGAAGTATGGATGAAGGGGATGATTACCTTCCAGCCACCACACCTGCACTGGAGAAAGCCCCTTGCTTCCCTGTAGAAAGCAAAGAGGAAGAACAGATAGAGAAAGTAAAAGCTACTGGAGAAGAAAGTACTCAAGCGGAGATATCATGTGAGTCTCCTTTCCTAGCCAAAGATTTTTACAAAAATGGTACTGTCATGGCACCTGACCTTCCTGAAATGCTAGATCTGGCAGGCACAAGGTCAAGATTGGCTTCTGTGAGTGCAGATGCTGAGGTTGCCAGGAGGAAATCAGTCCCATCAGAGACTGTGGTTGAGGATAGTCGTACTGGCTTGCCCCCGGTAACTGATGAAAACCATGTCATTGTAAAAACGGACAGTCAGCTCGAAGACCTGGGCTACTGTGTGTTCAATAAGTACACAGTCCCATTGCCATCACCTGTTCAAGACAGTGAGAATTTATCAGGGGAGAGTGGTACCTTTTACGAAGGCACTGATGATAAAGTTCGAAGAGATTTGGCCACAGACCTTTCACTGATTGAAGTGAAACTGGCAGCAGCCGGAAGAGTCAAAGATGAGTTCAGTGTTGACAAAGAAGCATCCGCGCATATCTCTGGTGACAAATCAGGACTGAGTAAGGAGTTTGACCAAGAGAAGAAAGCTAATGATAGGTTGGATACTGTACTAGAAAAGAGTGAAGAACATGCTGATTCAAAAGAACATGCCAAGAAAACTGAAGAGGCTGGTGATGAAATAGAAACATTCGGATTAGGAGTAACCTATGAGCAAGCTTTGGCCAAAGATTTGTCAATACCAACAGATGCATCCTCTGAGAAAGCAGAGAAGGGTCTTAGTTCAGTGCCAGAGATAGCTGAGGTAGAACCATCCAAAAAGGTGGAACAAGGTCTGGATTTTGCTGTCCAGGGTCAACTAGATGTTAAAATTAGTGACTTTGGACAGATGGCTTCAGGGCTAAACATAGATGATAGAAGGGCAACAGAGCTAAAACTTGAGGCTACACAGGACATGACCCCCTCATCCAAAGCACCGCAGGAGGCAGATGCATTTATGGGTGTTGAGTCTGGCCACATGAAAGAAGGCACTAAAGTTAGTGAGACAGAAGTCAAAGAGAAGGTGGCCAAGCCTGACTTGGTGCACCAGGAGGCTGTAGACAAGGAGGAGTCCTATGAATCTAGTGGTGAGCATGAAAGTCTCACCATGGAGTCCTTGAAAGCTGATGAGGGCAAGAAGGAAACATCTCCAGAATCATCTCTAATTCAAGATGAGATTGCCGTCAAATTGTCAGTGGAAATACCTTGCCCACCTGCTGTTTCAGAGGCTGATTTAGCCACAGATGAGAGAGCTGATGTCCAGATGGAATTTATTCAGGGGCCAAAAGAAGAAAGCAAAGAGACCCCAGATATATCCATCACGCCTTCTGATGTTGCAGAGCCATTGCATGAAACGATCGTATCTGAACCAGCAGAGATTCAGAGTGAGGAAGAAGAGATAGAAGCCCAGGGAGAATATGATAAACTGCTCTTCCGCTCAGACACCCTTCAGATAACTGACCTGGGTGTCTCAGGTGCCAGGGAGGAATTTGTGGAGACCTGCCCAAGTGAACACAAAGGAGTGATTGAGTCTGTTGTGACCATCGAGGATGATTTCATCACTGTAGTGCAAACCACAACTGATGAAGGGGAGTCAGGGTCCCACAGCGTGCGTTTTGCAGCCCTAGAGCAGCCTGAGGTGGAAAGGAGACCATCTCCTCATGATGAAGAAGAGTTTGAAGTAGAAGAGGCAGCTGAAGCCCAGGCAGAACCCAAAGATGGTTCCCCAGAGGCTCCAGCTTCCCCTGAGAGAGAAGAGGTTGCACTTTCTGAATATAAGACAGAAACCTATGACGATTACAAAGATGAGACCACCATTGACGACTCCATCATGGACGCTGACAGCCTCTGGGTGGACACTCAAGATGATGATAGGAGCATCATGACAGAACAGTTAGAAACTATTCCTAAAGAGGAGAAAGCTGAAAAGGAAGCTCGGAGATCATCTCTTGAGAAACATAGAAAAGAAAAGCCTTTTAAAACCGGGAGAGGCAGAATTTCCACTCCTGAAAGAAAAGTAGCTAAAAAGGAACCTAGCACAGTCTCCAGAGATGAAGTGAGAAGGAAAAAAGCAGTTTATAAGAAGGCTGAACTTGCTAAAAAAACAGAAGTTCAGGCCCACTCTCCCTCCAGGAAATTCATTTTAAAACCTGCTATCAAATATACTAGACCAACTCATCTCTCCTGTGTTAAGCGGAAAACCACAGCAGCAGGTGGGGAATCAGCTCTGGCTCCCAGTGTATTTAAACAGGCAAAGGACAAAGTCTCTGACGGAGTAACCAAGAGCCCAGAAAAGCGCTCTTCTCTCCCAAGACCTTCCTCCATTCTCCCTCCTCGGCGAGGTGTGTCAGGAGACAGAGATGAGAATTCCTTCTCTCTCAACAGTTCTATCTCTTCTTCAGCACGGCGGACCACCAGGTCAGAGCCAATTCGCAGAGCAGGGAAGAGTGGTACCTCAACACCCACTACCCCTGGGTCTACTGCCATCACTCCTGGCACCCCACCAAGTTATTCTTCACGCACACCAGGCACTCCTGGAACCCCTAGCTATCCCAGGACCCCTCACACACCAGGAACCCCCAAGTCTGCCATCTTGGTGCCGAGTGAGAAGAAGGTCGCCATCATACGTACTCCTCCAAAATCTCCTGCGACTCCCAAGCAGCTTCGGCTTATTAACCAACCACTGCCAGACCTGAAGAATGTCAAATCCAAAATCGGATCAACAGACAACATCAAATACCAGCCTAAAGGGGGGCAGGTACAAATTGTTACCAAGAAAATAGACCTAAGCCATGTGACATCCAAATGTGGCTCTCTGAAGAACATCCGCCACAGGCCAGGTGGCGGACGTGTGAAAATTGAGAGTGTAAAACTAGATTTCAAAGAAAAGGCCCAAGCTAAAGTTGGTTCTCTTGATAATGCTCATCATGTACCTGGAGGTGGTAATGTCAAGATTGACAGCCAAAAGTTGAACTTCAGAGAGCATGCTAAAGCCCGTGTGGACCATGGGGCTGAGATCATTACACAGTCCCCAGGCAGATCCAGCGTGGCATCACCCCGACGACTCAGCAATGTCTCCTCGTCTGGAAGCATCAACCTGCTCGAATCTCCTCAGCTTGCCACTTTGGCTGAGGATGTCACTGCTGCACTCGCTAAGCAGGGCTTGTGA", + "cds_sha256": "c62e3f6a2a5cc95f511e8933b70b3a2c3e23b74b157ed49d5ecaa856d11ececa", + "protein": "MADERKDEAKAPHWTSAPLTEASAHSHPPEIKDQGGAGEGLVRSANGFPYREDEEGAFGEHGSQGTYSNTKENGINGELTSADRETAEEVSARIVQVVTAEAVAVLKGEQEKEAQHKDQTAALPLAAEETANLPPSPPPSPASEQTVTVEEDLLTASKMEFHDQQELTPSTAEPSDQKEKESEKQSKPGEDLKHAALVSQPETTKTYPDKKDMQGTEEEKAPLALFGHTLVASLEDMKQKTEPSLVVPGIDLPKEPPTPKEQKDWFIEMPTEAKKDEWGLVAPISPGPLTPMREKDVFDDIPKWEGKQFDSPMPSPFQGGSFTLPLDVMKNEIVTETSPFAPAFLQPDDKKSLQQTSGPATAKDSFKIEEPHEAKPDKMAEAPPSEAMTLPKDAHIPVVEEHVMGKVLEEEKEAINQETVQQRDTFTPSGQEPILTEKETELKLEEKTTISDKEAVPKESKPPKPADEEIGIIQTSTEHTFSEQKDQEPTTDMLKQDSFPVSLEQAVTDSAMTSKTLEKAMTEPSALIEKSSIQELFEMRVDDKDKIEGVGAATSAELDMPFYEDKSGMSKYFETSALKEEATKSIEPGSDYYELSDTRESVHESIDTMSPMHKNGDKEFQTGKESQPSPPAQEAGYSTLAQSYPSDLPEEPSSPQERMFTIDPKVYGEKRDLHSKNKDDLTLSRSLGLGGRSAIEQRSMSINLPMSCLDSIALGFNFGRGHDLSPLASDILTNTSGSMDEGDDYLPATTPALEKAPCFPVESKEEEQIEKVKATGEESTQAEISCESPFLAKDFYKNGTVMAPDLPEMLDLAGTRSRLASVSADAEVARRKSVPSETVVEDSRTGLPPVTDENHVIVKTDSQLEDLGYCVFNKYTVPLPSPVQDSENLSGESGTFYEGTDDKVRRDLATDLSLIEVKLAAAGRVKDEFSVDKEASAHISGDKSGLSKEFDQEKKANDRLDTVLEKSEEHADSKEHAKKTEEAGDEIETFGLGVTYEQALAKDLSIPTDASSEKAEKGLSSVPEIAEVEPSKKVEQGLDFAVQGQLDVKISDFGQMASGLNIDDRRATELKLEATQDMTPSSKAPQEADAFMGVESGHMKEGTKVSETEVKEKVAKPDLVHQEAVDKEESYESSGEHESLTMESLKADEGKKETSPESSLIQDEIAVKLSVEIPCPPAVSEADLATDERADVQMEFIQGPKEESKETPDISITPSDVAEPLHETIVSEPAEIQSEEEEIEAQGEYDKLLFRSDTLQITDLGVSGAREEFVETCPSEHKGVIESVVTIEDDFITVVQTTTDEGESGSHSVRFAALEQPEVERRPSPHDEEEFEVEEAAEAQAEPKDGSPEAPASPEREEVALSEYKTETYDDYKDETTIDDSIMDADSLWVDTQDDDRSIMTEQLETIPKEEKAEKEARRSSLEKHRKEKPFKTGRGRISTPERKVAKKEPSTVSRDEVRRKKAVYKKAELAKKTEVQAHSPSRKFILKPAIKYTRPTHLSCVKRKTTAAGGESALAPSVFKQAKDKVSDGVTKSPEKRSSLPRPSSILPPRRGVSGDRDENSFSLNSSISSSARRTTRSEPIRRAGKSGTSTPTTPGSTAITPGTPPSYSSRTPGTPGTPSYPRTPHTPGTPKSAILVPSEKKVAIIRTPPKSPATPKQLRLINQPLPDLKNVKSKIGSTDNIKYQPKGGQVQIVTKKIDLSHVTSKCGSLKNIRHRPGGGRVKIESVKLDFKEKAQAKVGSLDNAHHVPGGGNVKIDSQKLNFREHAKARVDHGAEIITQSPGRSSVASPRRLSNVSSSGSINLLESPQLATLAEDVTAALAKQGL", + "protein_sha256": "1f134e0d2480079a47dae31ef1ad0d7d4dcb7694b679c186a886ad8ad4862e2f", + "genomic_window": { + "accession": "NC_000002.12", + "start": 209694750, + "end": 209694830, + "cds_offset": 2579, + "url": "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=NC_000002.12&rettype=fasta&retmode=text&seq_start=209694750&seq_stop=209694830", + "response_sha256": "a22367c53bf15f7e9351878979f40d637374a53d3c4078168f38ead881af1bdb", + "sequence": "GGACAGTCAGCTCGAAGACCTGGGCTACTGTGTGTTCAATAAGTACACAGTCCCATTGCCATCACCTGTTCAAGACAGTGA", + "sequence_sha256": "1645239ae61beccb1949e33230965042ff79ccce6f46c62e3be5d8382832be85" + } + } + }, + "expectations": { + "NTF3_first": { + "gene": "NTF3", + "cds_edits": [ + [ + 166, + "A", + "G" + ] + ], + "aa_pos": 55, + "ref_codon": "AAG", + "alt_codon": "AGG", + "protein_prefix_length": 55, + "protein_replacement": "R", + "protein_suffix_start": 56, + "protein_length": 257, + "protein_sha256": "f51f6deefc369a40ff1096f56e9cad1cdb714cfafd5291a2572b5da9a8797020", + "mutant_cds_sha256": "e746c47198f01e668c4aacfd915cdc3f339142f2ba9e0db0c7eabe8bf036a3c5", + "stop_codon": "TGA", + "stop_codon_offset": 771 + }, + "NTF3_second": { + "gene": "NTF3", + "cds_edits": [ + [ + 167, + "G", + "T" + ] + ], + "aa_pos": 55, + "ref_codon": "AAG", + "alt_codon": "AAT", + "protein_prefix_length": 55, + "protein_replacement": "N", + "protein_suffix_start": 56, + "protein_length": 257, + "protein_sha256": "b016bc597c727b1f4f2d851bc75dba01f8dd378b525afd5082914a52a3b8439f", + "mutant_cds_sha256": "ebbfa6029f0042df74f1ab65cd8fc02efbe8b8f4d8521d71bfe5b1dbfeb9952d", + "stop_codon": "TGA", + "stop_codon_offset": 771 + }, + "NTF3_compound": { + "gene": "NTF3", + "cds_edits": [ + [ + 166, + "AG", + "GT" + ] + ], + "aa_pos": 55, + "ref_codon": "AAG", + "alt_codon": "AGT", + "protein_prefix_length": 55, + "protein_replacement": "S", + "protein_suffix_start": 56, + "protein_length": 257, + "protein_sha256": "7de2ab6041ae304fa0eb320572fffd7e8e78b04e7eb6fd2ef3c85fde993af9a7", + "mutant_cds_sha256": "25704c5891ce21c53e6f8d5b3e599cf122126e1d727afef22cd048b300e7f308", + "stop_codon": "TGA", + "stop_codon_offset": 771 + }, + "NTF3_companion": { + "gene": "NTF3", + "cds_edits": [ + [ + 251, + "G", + "A" + ] + ], + "aa_pos": 83, + "ref_codon": "CCG", + "alt_codon": "CCA", + "protein_prefix_length": 0, + "protein_replacement": "", + "protein_suffix_start": 0, + "protein_length": 257, + "protein_sha256": "51d8e5a0dcb4df0c94896dd5fab8ca5deeb4523aaee8ba6ec6948763443b91e2", + "mutant_cds_sha256": "12efcb1f2fbf01e70bbb34e1028989f4f7568d91d5e63cfb5d8230f7d1397e50", + "stop_codon": "TGA", + "stop_codon_offset": 771 + }, + "MAP2_first": { + "gene": "MAP2", + "cds_edits": [ + [ + 2598, + "C", + "A" + ] + ], + "aa_pos": 866, + "ref_codon": "CTG", + "alt_codon": "ATG", + "protein_prefix_length": 866, + "protein_replacement": "M", + "protein_suffix_start": 867, + "protein_length": 1827, + "protein_sha256": "826410140bc91d4cc7d47852463d60e832da3777ea754a858bbb2ef5cd38703b", + "mutant_cds_sha256": "7bba8f272c2e5f15fcd995c3819ffe3d06de7fd3a2b4a405ca4716a3e74a4229", + "stop_codon": "TGA", + "stop_codon_offset": 5481 + }, + "MAP2_second": { + "gene": "MAP2", + "cds_edits": [ + [ + 2599, + "T", + "G" + ] + ], + "aa_pos": 866, + "ref_codon": "CTG", + "alt_codon": "CGG", + "protein_prefix_length": 866, + "protein_replacement": "R", + "protein_suffix_start": 867, + "protein_length": 1827, + "protein_sha256": "c243e8dea5562c37105d3b1804bd85971ac7581c60081b06e29d0dcde2996e5c", + "mutant_cds_sha256": "2a6c56d9ba0be398c9935dddfd621128035c42b3bf47e7ee683bc59748d9d882", + "stop_codon": "TGA", + "stop_codon_offset": 5481 + }, + "MAP2_deletion": { + "gene": "MAP2", + "cds_edits": [ + [ + 2602, + "GCTACTGTGTGTTCAATAAGTACACAGT", + "" + ] + ], + "aa_pos": 867, + "ref_codon": "GGC", + "alt_codon": "GCC", + "protein_prefix_length": 867, + "protein_replacement": "AHCHHLFKTVRIYQGRVVPFTKALMIKFEEIWPQTFH", + "protein_suffix_start": 1827, + "protein_length": 904, + "protein_sha256": "600dd17aafb9e8c5cfe55683e187d3c158105315e9d65567906055c830573f85", + "mutant_cds_sha256": "14b554e97c9f6834acb281583d574011e1ed6beb89579bf3ec53fe07989f339a", + "stop_codon": "TGA", + "stop_codon_offset": 2712 + }, + "MAP2_compound": { + "gene": "MAP2", + "cds_edits": [ + [ + 2598, + "CTGGGCTACTGTGTGTTCAATAAGTACACAGT", + "AGGG" + ] + ], + "aa_pos": 866, + "ref_codon": "CTG", + "alt_codon": "AGG", + "protein_prefix_length": 866, + "protein_replacement": "RAHCHHLFKTVRIYQGRVVPFTKALMIKFEEIWPQTFH", + "protein_suffix_start": 1827, + "protein_length": 904, + "protein_sha256": "a3385b6bfd321c5d98145c3981549d7a8a07a6cfbc7775536ccda16a23ae205c", + "mutant_cds_sha256": "48f5fb1d9dedbb10a5c0c59bde0b789c866bb051805cc512a7f4c20a77231c22", + "stop_codon": "TGA", + "stop_codon_offset": 2712 + }, + "MAP2_historical": { + "gene": "MAP2", + "cds_edits": [ + [ + 2598, + "CTGGGCTACTGTGTGTTCAATA", + "" + ] + ], + "aa_pos": 866, + "ref_codon": "CTG", + "alt_codon": "AGT", + "protein_prefix_length": 866, + "protein_replacement": "STQSHCHHLFKTVRIYQGRVVPFTKALMIKFEEIWPQTFH", + "protein_suffix_start": 1827, + "protein_length": 906, + "protein_sha256": "9e2b1813ce2fe5c5577e6281bbfce1590e49edb5444c0531a12401ba3f7d527d", + "mutant_cds_sha256": "7fd5fed764b4ddf3380aca50ccca79a2e8d68b6318f8a8c7516aeb25a965bc04", + "stop_codon": "TGA", + "stop_codon_offset": 2718 + } + } +} diff --git a/tests/test_osteosarc_rna_phasing.py b/tests/test_osteosarc_rna_phasing.py index 3846825..adfcbfd 100644 --- a/tests/test_osteosarc_rna_phasing.py +++ b/tests/test_osteosarc_rna_phasing.py @@ -13,7 +13,11 @@ import pytest from pyensembl import cached_release -from varcode import MolecularPhaseResolver, RNAReadPhasingSource, Variant +from varcode import ( + EffectCollection, MolecularPhaseResolver, RNAReadPhasingSource, Variant, + VariantCollection, +) +from varcode.effects.effect_classes import HaplotypeEffect, Silent, Substitution SHORT_SOURCE = "2024.06.11.bostongene.align.tcga.protocol.dr32Aligned.sorted" @@ -21,8 +25,15 @@ SHORT_NTF3 = SHORT_SOURCE + ".NTF3-chr12-5494381-compound" LONG_NTF3 = LONG_SOURCE + ".NTF3-chr12-5494381-compound" LONG_EXOC4 = LONG_SOURCE + ".EXOC4-chr7-133274996" -MEMBERS = (SHORT_NTF3, LONG_NTF3, LONG_EXOC4) +MAP2_PREFIX = "isovar/osteosarc/figure_comparisons/corpus/context-evidence.json.gz#/map2/sources/" +SHORT_MAP2_T1 = MAP2_PREFIX + "11/original_records" +SHORT_MAP2_T2 = MAP2_PREFIX + "12/original_records" +LONG_MAP2_T2 = MAP2_PREFIX + "3/original_records" +MAP2_MEMBERS = (SHORT_MAP2_T1, SHORT_MAP2_T2, LONG_MAP2_T2) +MEMBERS = (SHORT_NTF3, LONG_NTF3, LONG_EXOC4) + MAP2_MEMBERS MANIFEST_SHA256 = "b42dbce529cb5724153f035fca3783a41e673b31610a30851b2ef96d4ba43734" +REFERENCE_PATH = Path(__file__).parent / "data" / "sid_rna_effect_references.json" +REFERENCE_SHA256 = "d9e21cc96450d4a14e7b41d0900df48792c9f64b103b9ce01ab5cc841508d29e" NTF3 = ("12", 5494381, "A", "G") NTF3_COMPOUND_PARTNER = ("12", 5494382, "G", "T") @@ -30,12 +41,21 @@ NTF3_COMPANION = ("12", 5494466, "G", "A") EXOC4 = ("7", 133274996, "G", "T") EXOC4_COMPANION = ("7", 133895694, "G", "A") +MAP2_PARTS = ( + ("2", 209694769, "C", "A"), + ("2", 209694770, "T", "G"), + ("2", 209694773, "GCTACTGTGTGTTCAATAAGTACACAGT", ""), +) +MAP2_COMPOUND = ("2", 209694768, "CCTGGGCTACTGTGTGTTCAATAAGTACACAGT", "CAGGG") +MAP2_HISTORICAL = ("2", 209694768, "CCTGGGCTACTGTGTGTTCAATA", "C") @pytest.fixture(scope="module") def shared_phasing_reads(tmp_path_factory): osteosarc = pytest.importorskip("osteosarc", reason="install .[test-data]") - pytest.importorskip("pysam", reason="install .[rna]") + pysam = pytest.importorskip("pysam", reason="install .[rna]") + from osteosarc.records import read_records, read_template + root = Path(osteosarc.fetch_bundle("openvax-v2")) manifest_bytes = (root / "manifest.json").read_bytes() assert hashlib.sha256(manifest_bytes).hexdigest() == MANIFEST_SHA256 @@ -44,15 +64,33 @@ def shared_phasing_reads(tmp_path_factory): compound = recipe["targets"]["NTF3-chr12-5494381-compound"] assert (compound["assembly"], compound["contig"], compound["position"], compound["ref"], compound["alt"]) == ("GRCh38", "chr12", 5494381, "AG", "GT") - paths = osteosarc.export_bundle( - "openvax-v2", tmp_path_factory.mktemp("sid-rna-phasing"), members=MEMBERS) - return paths, root, json.loads(manifest_bytes) + destination = tmp_path_factory.mktemp("sid-rna-phasing") + paths = osteosarc.export_bundle("openvax-v2", destination, members=MEMBERS) + manifest = json.loads(manifest_bytes) + # The legacy MAP2 members are regional acquisitions. Include every other + # original record of their selected templates that is available in the + # bundle, including distant mates and supplementary placements. + for index, member in enumerate(MAP2_MEMBERS): + selected = {record.template for record in read_records(paths[member])} + source = manifest["members"][member]["source"] + complete_path = destination / ("complete-map2-%d.bam" % index) + with pysam.AlignmentFile(root / manifest["sources"][source]["bam"], "rb") as bam: + with pysam.AlignmentFile(complete_path, "wb", template=bam) as output: + for read in bam: + if read_template(read) in selected: + output.write(read) + pysam.index(str(complete_path)) + paths[member] = complete_path + return paths, root, manifest @pytest.mark.parametrize("member,records,templates,supplementary", [ (SHORT_NTF3, 12, 6, 0), (LONG_NTF3, 15, 15, 0), (LONG_EXOC4, 40, 38, 2), + (SHORT_MAP2_T1, 125, 78, 0), + (SHORT_MAP2_T2, 26, 18, 0), + (LONG_MAP2_T2, 33, 32, 1), ]) def test_shared_phasing_inputs_keep_complete_templates( shared_phasing_reads, member, records, templates, supplementary): @@ -61,13 +99,18 @@ def test_shared_phasing_inputs_keep_complete_templates( paths, root, manifest = shared_phasing_reads exported = list(read_records(paths[member])) identities = Counter(record.digest for record in exported) - assert identities == Counter(manifest["members"][member]["records"]) + member_identities = Counter(manifest["members"][member]["records"]) + assert not member_identities - identities + if member not in MAP2_MEMBERS: + assert identities == member_identities selected = {record.template for record in exported} source = manifest["members"][member]["source"] + source_records = list(read_records(root / manifest["sources"][source]["bam"])) + assert selected == {record.template for record in source_records + if record.digest in member_identities} # Compare to every record of each selected RG/QNAME in the bundle's source, # rather than refetching a locus and losing mates or distant split pieces. - complete = Counter(record.digest for record in read_records( - root / manifest["sources"][source]["bam"]) if record.template in selected) + complete = Counter(record.digest for record in source_records if record.template in selected) assert identities == complete assert len(exported) == records assert len(selected) == templates @@ -107,6 +150,226 @@ def test_real_sid_rna_fragment_evidence( source.close() +@pytest.fixture(scope="module") +def effect_references(): + """Offline RefSeq references and separately audited protein expectations.""" + payload = REFERENCE_PATH.read_bytes() + assert hashlib.sha256(payload).hexdigest() == REFERENCE_SHA256 + return json.loads(payload) + + +@pytest.mark.parametrize("gene", ["NTF3", "MAP2"]) +def test_primary_effect_references_match_the_selected_transcripts(effect_references, gene): + reference = effect_references["references"][gene] + transcript = cached_release(81).transcript_by_id(reference["transcript_id"]) + assert transcript.coding_sequence == reference["cds"] + assert transcript.protein_sequence == reference["protein"] + assert hashlib.sha256(transcript.coding_sequence.encode()).hexdigest() == reference["cds_sha256"] + assert hashlib.sha256(transcript.protein_sequence.encode()).hexdigest() == reference["protein_sha256"] + window = reference["genomic_window"] + # Independently acquired chromosome bases map to this exact CDS interval. + start = min(transcript.start_codon_spliced_offsets) + assert transcript.spliced_offset(window["start"]) - start == window["cds_offset"] + assert transcript.sequence[transcript.spliced_offset(window["start"]): + transcript.spliced_offset(window["end"]) + 1] == window["sequence"] + + +def transcript_effects(effects, transcript_id, *, joint=False): + return [effect for effect in effects + if getattr(getattr(effect, "transcript", None), "id", None) == transcript_id + and (len(getattr(effect, "variants", ())) > 1) is joint] + + +def individual_signatures(effects): + """Snapshot all individual predictions, including other transcript models.""" + return sorted((str(effect.variant), getattr(getattr(effect, "transcript", None), "id", ""), + type(effect).__name__, effect.short_description, effect.mutant_protein_sequence) + for effect in effects if len(getattr(effect, "variants", ())) < 2) + + +def assert_primary_protein(effect, case, references): + expected = references["expectations"][case] + reference = references["references"][expected["gene"]] + protein = (reference["protein"][:expected["protein_prefix_length"]] + + expected["protein_replacement"] + + reference["protein"][expected["protein_suffix_start"]:]) + assert len(protein) == expected["protein_length"] + assert hashlib.sha256(protein.encode()).hexdigest() == expected["protein_sha256"] + if case == "NTF3_companion": + # Silent effects intentionally expose no mutant protein sequence. + assert isinstance(effect, Silent) + assert effect.aa_pos == 83 and effect.aa_ref == "P" + assert effect.modifies_protein_sequence is False + return + assert effect.mutant_protein_sequence == protein + + +def assert_primary_cdna(effect, case, references): + """Verify an available joint sequence product against the RefSeq edit audit.""" + expected = references["expectations"][case] + reference = references["references"][expected["gene"]] + if effect.mutant_transcript is not None: + cdnas = [effect.mutant_transcript.cdna_sequence] + else: + cdnas = [outcome.mutant.cdna_sequence + for candidate in effect.realized_candidates for outcome in candidate.outcomes] + assert cdnas + start = min(effect.transcript.start_codon_spliced_offsets) + delta = sum(len(alt) - len(ref) for _, ref, alt in expected["cds_edits"]) + for cdna in cdnas: + cds = cdna[start:start + len(reference["cds"]) + delta] + assert hashlib.sha256(cds.encode()).hexdigest() == expected["mutant_cds_sha256"] + assert cds[expected["aa_pos"] * 3:expected["aa_pos"] * 3 + 3] == expected["alt_codon"] + stop = expected["stop_codon_offset"] + assert cds[stop:stop + 3] == expected["stop_codon"] + + +@pytest.mark.parametrize("annotator", ["fast", "protein_diff", "transcript_model"]) +@pytest.mark.parametrize("member", [SHORT_NTF3, LONG_NTF3]) +def test_real_ntf3_cis_reads_predict_the_compound_codon( + shared_phasing_reads, effect_references, annotator, member): + paths, _, _ = shared_phasing_reads + variants = tuple(Variant(*allele, genome=cached_release(81)) + for allele in (NTF3, NTF3_COMPOUND_PARTNER)) + collection = VariantCollection(variants) + baseline = individual_signatures(collection.effects(annotator=annotator)) + source = RNAReadPhasingSource( + str(paths[member]), variants=variants, + phasing_error_rate=0.01 if member == SHORT_NTF3 else 0.05) + try: + resolver = MolecularPhaseResolver(source) + assert resolver.in_cis(*variants) is True + effects = collection.effects(annotator=annotator, phase_resolver=resolver) + assert individual_signatures(effects) == baseline + transcript_id = effect_references["references"]["NTF3"]["transcript_id"] + joint, = transcript_effects(effects, transcript_id, joint=True) + assert isinstance(joint, Substitution if annotator == "transcript_model" else HaplotypeEffect) + assert joint.variants == variants + assert joint.phase_source == "rna_reads" + assert joint.annotator == effects.annotator == annotator + assert joint.annotator_version == effects.annotator_version + assert_primary_protein(joint, "NTF3_compound", effect_references) + assert_primary_cdna(joint, "NTF3_compound", effect_references) + individuals = {effect.variant: effect for effect in transcript_effects(effects, transcript_id)} + for variant, case in zip(variants, ("NTF3_first", "NTF3_second")): + assert_primary_protein(individuals[variant], case, effect_references) + restored = EffectCollection.from_json(effects.to_json()) + restored_joint, = transcript_effects(restored, transcript_id, joint=True) + assert restored_joint.variants == variants + assert restored_joint.phase_source == joint.phase_source + assert_primary_protein(restored_joint, "NTF3_compound", effect_references) + finally: + source.close() + + +@pytest.mark.parametrize("annotator", ["fast", "protein_diff", "transcript_model"]) +@pytest.mark.parametrize("alleles,cases,minimum,phase", [ + ((NTF3, NTF3_COMPANION), ("NTF3_first", "NTF3_companion"), 2, False), + ((NTF3, NTF3_COMPOUND_PARTNER), ("NTF3_first", "NTF3_second"), 10, None), +]) +def test_real_ntf3_trans_and_unknown_reads_keep_separate_predictions( + shared_phasing_reads, effect_references, annotator, alleles, cases, minimum, phase): + paths, _, _ = shared_phasing_reads + variants = tuple(Variant(*allele, genome=cached_release(81)) for allele in alleles) + collection = VariantCollection(variants) + baseline = individual_signatures(collection.effects(annotator=annotator)) + source = RNAReadPhasingSource( + str(paths[LONG_NTF3]), variants=variants, min_alt_reads=minimum, phasing_error_rate=0.05) + try: + resolver = MolecularPhaseResolver(source) + assert resolver.in_cis(*variants) is phase # Unknown must remain distinct from trans. + assert all(not resolver.phased_partners(variant) for variant in variants) + effects = collection.effects(annotator=annotator, phase_resolver=resolver) + assert individual_signatures(effects) == baseline + assert not any(len(getattr(effect, "variants", ())) > 1 for effect in effects) + transcript_id = effect_references["references"]["NTF3"]["transcript_id"] + individuals = {effect.variant: effect for effect in transcript_effects(effects, transcript_id)} + assert set(individuals) == set(variants) + for variant, case in zip(variants, cases): + assert_primary_protein(individuals[variant], case, effect_references) + finally: + source.close() + + +@pytest.mark.parametrize("annotator", ["fast", "protein_diff"]) +@pytest.mark.parametrize("member,support,phase", [ + (SHORT_MAP2_T1, 4, True), (SHORT_MAP2_T2, 2, True), (LONG_MAP2_T2, 1, None), +]) +def test_real_map2_anchored_reads_predict_the_complete_complex_allele( + shared_phasing_reads, effect_references, annotator, member, support, phase): + paths, _, _ = shared_phasing_reads + variants = tuple(Variant(*allele, genome=cached_release(81)) for allele in MAP2_PARTS) + collection = VariantCollection(variants) + baseline = individual_signatures(collection.effects(annotator=annotator)) + source = RNAReadPhasingSource( + str(paths[member]), variants=variants, + phasing_error_rate=0.05 if member == LONG_MAP2_T2 else 0.01) + try: + source.register_haplotype(variants) + assert tuple(source.supports_variant(variant) for variant in variants) == (support,) * 3 + resolver = MolecularPhaseResolver(source) + for left, right in ((0, 1), (0, 2), (1, 2)): + assert resolver.in_cis(variants[left], variants[right]) is phase + effects = collection.effects(annotator=annotator, phase_resolver=resolver) + assert individual_signatures(effects) == baseline + transcript_id = effect_references["references"]["MAP2"]["transcript_id"] + individuals = {effect.variant: effect for effect in transcript_effects(effects, transcript_id)} + assert set(individuals) == set(variants) + for variant, case in zip(variants, ("MAP2_first", "MAP2_second", "MAP2_deletion")): + assert_primary_protein(individuals[variant], case, effect_references) + if phase is None: + assert not any(len(getattr(effect, "variants", ())) > 1 for effect in effects) + assert not any(resolver.has_evidence(variant) for variant in variants) + else: + joint, = transcript_effects(effects, transcript_id, joint=True) + assert joint.variants == variants + assert joint.phase_source == "rna_reads" + assert joint.annotator == effects.annotator == annotator + assert_primary_protein(joint, "MAP2_compound", effect_references) + assert_primary_cdna(joint, "MAP2_compound", effect_references) + restored_joint, = transcript_effects( + EffectCollection.from_json(effects.to_json()), transcript_id, joint=True) + assert restored_joint.variants == variants + assert restored_joint.phase_source == joint.phase_source + assert_primary_protein(restored_joint, "MAP2_compound", effect_references) + assert_primary_cdna(restored_joint, "MAP2_compound", effect_references) + finally: + source.close() + + +@pytest.mark.parametrize("annotator", ["fast", "protein_diff"]) +@pytest.mark.parametrize("member,support", [ + (SHORT_MAP2_T1, 4), (SHORT_MAP2_T2, 2), (LONG_MAP2_T2, 1), +]) +def test_real_map2_rna_does_not_conflate_competing_dna_alleles( + shared_phasing_reads, effect_references, annotator, member, support): + paths, _, _ = shared_phasing_reads + variants = tuple(Variant(*allele, genome=cached_release(81)) + for allele in (MAP2_COMPOUND, MAP2_PARTS[-1], MAP2_HISTORICAL)) + collection = VariantCollection(variants) + baseline = individual_signatures(collection.effects(annotator=annotator)) + source = RNAReadPhasingSource( + str(paths[member]), variants=variants, + phasing_error_rate=0.05 if member == LONG_MAP2_T2 else 0.01) + try: + for variant in variants: + source.register_haplotype([variant]) + assert tuple(source.supports_variant(variant) for variant in variants) == (support, 0, 0) + resolver = MolecularPhaseResolver(source) + assert resolver.in_cis(variants[0], variants[1]) is None + assert resolver.in_cis(variants[0], variants[2]) is None + effects = collection.effects(annotator=annotator, phase_resolver=resolver) + assert individual_signatures(effects) == baseline + assert not any(len(getattr(effect, "variants", ())) > 1 for effect in effects) + transcript_id = effect_references["references"]["MAP2"]["transcript_id"] + individuals = {effect.variant: effect for effect in transcript_effects(effects, transcript_id)} + assert set(individuals) == set(variants) + for variant, case in zip(variants, ("MAP2_compound", "MAP2_deletion", "MAP2_historical")): + assert_primary_protein(individuals[variant], case, effect_references) + finally: + source.close() + + def test_real_sid_cis_call_requires_enough_fragments(shared_phasing_reads): paths, _, _ = shared_phasing_reads variants = tuple(Variant(*allele, genome=cached_release(81)) diff --git a/varcode/version.py b/varcode/version.py index 6c37e06..900edd1 100644 --- a/varcode/version.py +++ b/varcode/version.py @@ -1 +1 @@ -__version__ = "10.11.5" +__version__ = "10.11.6"