Retain original SiHa RNA for reproducible viral coding regression

examples/retain_siha_fixture.py

@@ -0,0 +1,97 @@
+"""Retain the existing SiHa regression slice, without downloading or realigning."""
+
+import argparse
+import gzip
+from hashlib import sha256
+import json
+from pathlib import Path
+from tempfile import TemporaryDirectory
+
+import pysam
+
+from isovar import reconstruct_viral_sequences
+
+
+TARGET = 'isovar_viral_a22b9eba3846fcf4496e065f'
+SOURCE = 'SRR12650977:SiHa-1:bounded-original-RNA'
+ROOT = Path(__file__).resolve().parents[1]
+FASTA = ROOT / 'tests/data/viral/hpv16-GDC.fa.gz'
+
+
+def digest(value):
+    return sha256(value).hexdigest()
+
+
+def retain(bam_path, replay_path, output):
+    bam_path, replay_path, output = Path(bam_path), Path(replay_path), Path(output)
+    receipt_bytes = replay_path.read_bytes()
+    receipt = json.loads(receipt_bytes)
+    bam_sha = digest(bam_path.read_bytes())
+    if receipt['input']['sha256'] != bam_sha:
+        raise ValueError('Original BAM does not match the historical replay receipt')
+    names = {row['fragment'][1] for row in receipt['support']['witnesses']
+             if any(TARGET in segment['overlapping_sequences'] for segment in row['segments'])
+             and not row['competing_placement_outside_targets']}
+    output.mkdir(parents=True, exist_ok=False)
+    with pysam.AlignmentFile(str(bam_path)) as source:
+        header = source.header.to_dict()
+        sam = [record.to_string() for record in source if record.query_name in names]
+        (output / 'source-header.sam.gz').write_bytes(gzip.compress(str(source.header).encode(), mtime=0))
+        header['SQ'] = [row for row in header['SQ'] if row['SN'] == 'HPV16']
+        with pysam.AlignmentFile(str(output / 'original-rna.bam'), 'wb', header=header) as destination:
+            for text in sam:
+                # Parsing fails if a selected mate/alternative uses any other
+                # contig. No record is silently dropped or its fields changed.
+                record = pysam.AlignedSegment.fromstring(text, destination.header)
+                if record.to_string() != text:
+                    raise ValueError('A selected original SAM record changed')
+                destination.write(record)
+    (output / 'full-input-replay.json.gz').write_bytes(gzip.compress(receipt_bytes, mtime=0))
+    with TemporaryDirectory() as tmp:
+        tmp = Path(tmp)
+        fasta = tmp / 'reference.fa'
+        fasta.write_bytes(gzip.decompress(FASTA.read_bytes()))
+        with pysam.AlignmentFile(str(output / 'original-rna.bam')) as bam:
+            result = reconstruct_viral_sequences(
+                bam, genome_fasta=fasta, sample_id='SiHa-1', source=SOURCE,
+                contigs=['HPV16'], reference_profile=None, max_discovery_reads=500,
+                max_paths=32, reference_cache=tmp / 'reference-cache')
+    target, = [row for row in result['sequences'] if row['id'] == TARGET]
+    coding, = [row for row in target['coding']
+               if any(ref['protein_id'] == 'NP_041326.1' for ref in row['reference_orfs'])]
+    manifest = dict(
+        source=dict(sample='GSM4785628', run='SRR12650977', project='PRJNA663801', series='GSE158033',
+                    sample_name='SiHa-1', library_prefix_pairs=1000000,
+                    url='https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM4785628'),
+        source_bam_sha256=bam_sha, source_result_sha256=digest(receipt_bytes),
+        source_fasta=receipt['reference']['files']['genome_fasta'],
+        selection='whole query-name groups with a segment overlapping the historical target; no recorded placement outside HPV16',
+        selected_query_names=sorted(names), selected_fragments=len(names), selected_records=len(sam),
+        sam_record_sha256=[digest(text.encode()) for text in sam], historical_target_id=TARGET,
+        replay=dict(source=SOURCE, contigs=['HPV16'], reference_profile=None,
+                    max_discovery_reads=500, max_paths=32, min_overlap=30),
+        reference_fasta=dict(path='../hpv16-GDC.fa.gz', sha256=digest(FASTA.read_bytes()),
+                             sequence_sha256=result['coding_reference']['contigs'][0]['sequence_sha256']),
+        orfdata=dict(version=result['coding_reference']['version'],
+                     catalogue_sha256=result['coding_reference']['contigs'][0]['reference']['catalogue_sha256'],
+                     source=result['coding_reference']['contigs'][0]['reference']['source']),
+        expected=dict(paths=len(result['sequences']), target_fragments=result['support']['target_fragments'],
+                      target_rna_sha256=digest(target['sequence'].encode()),
+                      target_rna_length=len(target['sequence']), target_protein_sha256=digest(coding['amino_acids'].encode()),
+                      target_protein_length=len(coding['amino_acids']),
+                      compatible_fragments=target['compatible_fragments'],
+                      supporting_fragments=target['supporting_fragments'],
+                      supporting_segments=target['supporting_segments']),
+        files={path.name: digest(path.read_bytes()) for path in sorted(output.iterdir())})
+    (output / 'manifest.json').write_text(json.dumps(manifest, indent=2) + '\n')
+    return manifest
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument('--bam', required=True)
+    parser.add_argument('--replay', required=True)
+    parser.add_argument('--output-dir', required=True)
+    args = parser.parse_args()
+    manifest = retain(args.bam, args.replay, args.output_dir)
+    print(json.dumps(manifest['expected'], indent=2))

isovar/__init__.py

@@ -10,7 +10,7 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
-__version__ = "1.64.0"
+__version__ = "1.64.1"
 
 
 from .allele_interpretations import reconcile_allele_interpretations

notes/siha-original-rna-fixture.md

@@ -0,0 +1,19 @@
+# Original SiHa RNA fixture (#531)
+
+Retain a bounded original-read selection from the already acquired public
+SRR12650977 competitive alignment. Select complete query-name groups whose
+recorded placements are all on HPV16 and support the previously reported local
+path. Preserve every selected SAM record, mate, flag, quality and alignment tag.
+Remove only unused header sequence entries; retain the complete source header
+and exact source BAM/result/reference hashes separately.
+
+Use the existing sequence-pinned HPV16 FASTA for a small custom-reference replay.
+Do not claim that this slice repeats whole-library competition: retain the
+separate historical full-input receipt (26 paths, 71 target fragments) and its
+human-competitive placements. Pin the ORFdata catalogue and source CDS identities.
+Assert the existing local RNA/coding prediction, bounded support and ambiguity;
+do not infer protein expression, a full viral transcript or viral absence.
+
+Provide a deterministic extraction recipe, hashes, public-source attribution and
+redistribution policy, plus one offline API/CLI regression. No reconstruction
+algorithm or new biological sequence engineering is part of this change.

notes/viral-coding-validation.md

@@ -17,6 +17,13 @@ transcript or translation in the sample.
 
 ## Saved real RNA replay
 
+The bounded original-read regression now lives in
+[`tests/data/viral/siha`](../tests/data/viral/siha/README.md). Its 23-fragment
+slice retains the same selected coding prediction, with an offline CLI test,
+original SAM hashes, reference/catalogue pins and a deterministic extraction
+recipe. The unchanged full-input receipt is stored separately; its 71-fragment
+counts below are not counts from the bounded fixture.
+
 The existing public SiHa-1 RNA experiment uses GSM4785628 / SRR12650977,
 project PRJNA663801. Its one-million-pair prefix and competitive human/viral
 alignment are described in [the recovery experiment](https://github.com/openvax/isovar/issues/514).

tests/data/viral/siha/README.md

@@ -0,0 +1,63 @@
+# Original SiHa-1 RNA regression slice
+
+These are original public RNA observations from **GSM4785628 / SRR12650977**,
+GSE158033 / PRJNA663801, deposited by Zhi-Ming Zheng's NCI laboratory. The
+[GEO record](https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM4785628)
+identifies total-RNA sequencing of the SiHa cervical-cancer cell line.
+
+`original-rna.bam` retains 46 original alignment records from 23 query-name
+groups. The earlier experiment used a one-million-pair prefix and competitive
+alignment to GDC GRCh38.d1.vd1. Its exact source BAM, full FASTA and replay
+hashes are in `manifest.json`; `source-header.sam.gz` preserves the complete
+alignment header and minimap2 command/version. No read sequence, quality,
+mate flag, tag or selected alignment was changed. Only unused `@SQ` header
+entries were removed. Canonical SAM hashes verify every retained record.
+
+Selection uses the saved full-input witnesses: include every query-name group
+with at least one segment overlapping the existing 344-nt RNA path and no
+recorded placement outside HPV16, then retain **all** original BAM records for
+those names. The selection is intentionally conditioned on this one case. It
+does not estimate viral abundance or test whole-library specificity/absence.
+Three compatible fragments with competing human placements are outside this
+slice; their evidence remains in the separate full-input receipt.
+
+The bounded replay yields seven RNA paths from 23 fragments. The selected
+344-nt path reproduces the previously reported 98-aa E7 coding prediction and
+its L28F difference from NP_041326.1. Its coding bases are tiled by original
+RNA; no single segment spans the CDS. RNA strand and transcript completeness
+remain unknown. A complete coding prediction is not evidence of protein
+expression or an intact full-length viral transcript.
+
+`full-input-replay.json.gz` is the **unchanged historical receipt** for 26 paths
+and 71 target fragments, made with ORFdata 0.29.0. It must not be compared as
+though these were counts from the 23-fragment slice. The new replay separately
+pins the tested ORFdata 0.39.1 catalogue, HPV16 GenBank source and reference
+sequence hashes. A changed catalogue requires an explicit fixture review.
+
+The exact HPV16 FASTA already lives at `../hpv16-GDC.fa.gz`. The small replay
+uses `custom` reference validation because unused human contigs are absent.
+The original BAM header has no M5 checksum, and the replay preserves that
+unverified-BAM flag; the manifest separately pins the historical alignment
+reference. This is a subset replay, not a new competitive alignment.
+
+From the repository root:
+
+```sh
+gzip -cd tests/data/viral/hpv16-GDC.fa.gz > /tmp/siha-reference.fa
+isovar viral --bam tests/data/viral/siha/original-rna.bam \
+  --genome-fasta /tmp/siha-reference.fa --sample-id SiHa-1 \
+  --alignment-reference custom --contig HPV16 \
+  --max-discovery-reads 500 --max-paths 32 --output-dir /tmp/siha-replay
+```
+
+To reproduce the extraction from the saved original inputs, run
+`python examples/retain_siha_fixture.py --bam SOURCE.bam --replay FULL-REPLAY.json
+--output-dir NEW-DIRECTORY`. It validates the original BAM checksum, preserves
+all selected SAM fields, reruns the existing reconstruction and records hashes.
+No download or realignment is performed.
+
+The [NCBI molecular-data policy](https://www.ncbi.nlm.nih.gov/home/about/policies/#data)
+places no NCBI restrictions on molecular-data use or redistribution, while
+retaining the submitter-rights caveat. This public molecular-data subset is
+distributed with source attribution under that policy; it is not relicensed
+as project-authored sequence data. The project's code remains Apache-2.0.

tests/data/viral/siha/manifest.json

@@ -0,0 +1,134 @@
+{
+  "source": {
+    "sample": "GSM4785628",
+    "run": "SRR12650977",
+    "project": "PRJNA663801",
+    "series": "GSE158033",
+    "sample_name": "SiHa-1",
+    "library_prefix_pairs": 1000000,
+    "url": "https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM4785628"
+  },
+  "source_bam_sha256": "139748c3cd839c5edcf5358a693f52e7c0b59ecb9f0815bc90f260e87197208c",
+  "source_result_sha256": "cbcaa31f643cf926d506c3d72ce6a249aface8529282cdf2003f3cc0289413bd",
+  "source_fasta": {
+    "path": "/Users/iskander/code/isovar/.cache/gdc-viral-reconstruction-20261015/.cache/GRCh38.d1.vd1.fa",
+    "sha256": "fbde121e898f00cc064eaafe5acc58940d7fe069ad14d34acd40c5d675e6b69b"
+  },
+  "selection": "whole query-name groups with a segment overlapping the historical target; no recorded placement outside HPV16",
+  "selected_query_names": [
+    "SRR12650977.105184",
+    "SRR12650977.112911",
+    "SRR12650977.13194",
+    "SRR12650977.163131",
+    "SRR12650977.264995",
+    "SRR12650977.37489",
+    "SRR12650977.380150",
+    "SRR12650977.518757",
+    "SRR12650977.554077",
+    "SRR12650977.573820",
+    "SRR12650977.611919",
+    "SRR12650977.662935",
+    "SRR12650977.680139",
+    "SRR12650977.719197",
+    "SRR12650977.729885",
+    "SRR12650977.776851",
+    "SRR12650977.831234",
+    "SRR12650977.861040",
+    "SRR12650977.905283",
+    "SRR12650977.907102",
+    "SRR12650977.926848",
+    "SRR12650977.933521",
+    "SRR12650977.940631"
+  ],
+  "selected_fragments": 23,
+  "selected_records": 46,
+  "sam_record_sha256": [
+    "c39a2662d5953928bc199e98b2a3d3a8bb6efadfeea6a3625d18d42bcf593f92",
+    "54307c50a2df00f1cbf054315636afd39f816e4bdf95ba1628b339770d64ddaa",
+    "d1ad9d246351777ed5372cb25c5b7620e38a6689da1cfc7ed0556a4c01ee7a2c",
+    "25e6f0e464225bd4f6d8719f9ab2dc758ed53bb6619e86a08489053890c9a255",
+    "88eb42896b1718c0ff8b41cb6502fc9daa53a6f34c5776429fd2814debe4c3bb",
+    "65eb9d337f1c857a8d22cab6d9128784ded6d987ab59f83834fed0528a66effd",
+    "b0967d52492f57f380ea1e4178b26d1306ec38078935e9c5aea82bdec2640970",
+    "1e49c02ac1b1f545cbff277618336ee0f97da94ac4988836687a9b81a4ecdb8f",
+    "7e95209b07fb8fae02f7a952a2ae9be60f647dbc9c6a86dd29c6b0bc58a085c9",
+    "05ec523ea3c7c09ef43ae3a62148967b23a07915dc50e91f91df2c96da26b66b",
+    "69ae9be13a408ce919fb096a10cbac404aa8f53a62d086531a86bb233620f57e",
+    "a6392c3cbd11d39da6ce166bdc3ff69aab2ec67b8b87fa84ec0526eadf5fa4dd",
+    "ef46960d4cd94cd4acaa9d8e32cd35564ed93451eb00b277112a25c61f30bdde",
+    "cca7327e618eba32d860e938b2321818fe37eb6467e4ecd703d3f904a7ead22d",
+    "237fd5d479c1d6a53c6b28b636d2e2b3c1578e6c087081ede94b81a472dc419a",
+    "c4c59ad46a4671b8054b638317f7e5b586bcda02167a17944545a3f3e5a3c5fc",
+    "9bff0a9016c1eb67f2bd7764953429d45a00032d2aa5835bd230334cdd1bd656",
+    "1fb3c6ff7d04ee5bd1a66a732deab5eaf33df50c9c7edc7a5bc216ad90144d32",
+    "dd1b111ae18b4196d9201131c2a82378ec2133f247316b1f049303a918ba5ac3",
+    "3a9923489179f73f81f443235a3cc8edb607e2670a28ef3066eb58698f5da049",
+    "dcda6ff668fd22d4668532ead8f6eac5f01a607b6251f3b6429040b0ac294c38",
+    "7ce226081eea7c352cbec707a1fe01234c848327c9276c5477a2f7b8439c15c2",
+    "96ff1a638949344335b03dfdae7bd80637ad5984921a37dbeae2ce0fcdbeb251",
+    "6a3c7c891f097012bdfb49630f4733a65d7bc5c70b6be3202ded209e63ec4c47",
+    "e24e74bad46ce3782718ce5c124b8b6ea61086b2caaa52a07d58d1cf4a1e3555",
+    "5a5ce8de5a1e6750487ff81472b6f5abc61a47a671ce96697d3fe6082e4ac633",
+    "b106eb4d5142e474d2c6d74b57b6a35368065e63df9c05f791e040fa9230d038",
+    "c16c6ad4122f5a076b9e1f48e26ab92e15fa29920b2cfb41f669835b7a3ea275",
+    "d42f71f2212b14341ce961018e49dc8edebaa992970a397c36999ea92215b543",
+    "d39f86ec364e31ca6a9f87c8876d0e322051bd11d9e11b5d508dd0dc3b0ba7b6",
+    "3dd492725c237e14f79032cd5019a485f80b1060308261b7981f8a04f14bb1a7",
+    "9636bc0d54974e9def4126531d97fc204b92bec0a2cd7b6f93df639e41ead6cb",
+    "1dcf9a5706c9a2a5ff1254f1994ffdcde9cdc754a59c65562b0209bb8c4d9492",
+    "7d110f48de93e5bd53a8c29e5277ee18ceb0d05a06df06e4f6d4df6e2112bdae",
+    "4d8844472b37a7e407b0c07e2ed028ae97bba610f8455349c4e58dba89a5aded",
+    "1e8f8c78104e7b79c079dec1fdf93b72ecc878b0441568a4978c20832e746ba2",
+    "62a95edf90d1a2e79e4f316f31e7273738e9cc4d1051b6f3adcc6441a01ec803",
+    "915e0206b9ba00b5e5ca5ae00e75152e69ba176ec793f8d939e01e7f6daedd63",
+    "cebd3f64a488f0a2087541e49afda1345bca162175741d8528ec4674abbcc2db",
+    "adc28951bc0b19c0269ef2e5dadada54ae3b450a50c860ef73f77e34f51d2f46",
+    "56df33cb182add9eb712f01f53a977d5b5eacc5c93a9f65d5e5c1048cb41e0db",
+    "8493a48063add8d31d0c34037c03cbd207d6403f3a50cf8e2463681b82b6f63f",
+    "6fad366c30d3835a0f9e8f1af1fdb26f809157ba90b3f0cbc2ef64f51c4df6b9",
+    "c0d5f6cac5a4e4aad4d4dec4cdded5140418b537e97043c6c4009cebe954585d",
+    "6633c6dc8b8d42ef90a30a12c9dc667113c42f8abd75a43e26d398581ef6eb8d",
+    "e2ddfed9d83a403dbe200c2622a19b86e34a83fcf07488e7231dd869041cef35"
+  ],
+  "historical_target_id": "isovar_viral_a22b9eba3846fcf4496e065f",
+  "replay": {
+    "source": "SRR12650977:SiHa-1:bounded-original-RNA",
+    "contigs": [
+      "HPV16"
+    ],
+    "reference_profile": null,
+    "max_discovery_reads": 500,
+    "max_paths": 32,
+    "min_overlap": 30
+  },
+  "reference_fasta": {
+    "path": "../hpv16-GDC.fa.gz",
+    "sha256": "3359f936fb20d595b1ef3b91ee79909f4e5b0ef2dc67988e99fd32a5c9089609",
+    "sequence_sha256": "a9e63c3d3754c2b09942e16eac9f8c9092245b61647861e12397b7413b56e365"
+  },
+  "orfdata": {
+    "version": "0.39.1",
+    "catalogue_sha256": "aca0c72948d4cda896d14b32203981bce4759e0c70212c05dca6f3d0e16ed4b4",
+    "source": {
+      "accession": "NC_001526.4",
+      "sha256": "282be442f226cc809fefbda616ac136b8d20dfde23fad09558f3f653a3de96a4",
+      "url": "https://www.ncbi.nlm.nih.gov/sviewer/viewer.fcgi?id=NC_001526.4&db=nuccore&report=genbank&retmode=text"
+    }
+  },
+  "expected": {
+    "paths": 7,
+    "target_fragments": 23,
+    "target_rna_sha256": "6d0ea9cd27f96aa0a60221aeea6de5e58694c242995044fb1f0a59e9ad187af8",
+    "target_rna_length": 344,
+    "target_protein_sha256": "8c9403c9751757f2bef8e48f03645658d0015519e2998dd0fdd012b81fe8a78b",
+    "target_protein_length": 98,
+    "compatible_fragments": 8,
+    "supporting_fragments": 23,
+    "supporting_segments": 33
+  },
+  "files": {
+    "full-input-replay.json.gz": "ded43f3a27676c14efcf0655f6a3533525e92d34256fd668f7a8d3d3dec60b76",
+    "original-rna.bam": "a0414a94c9a1cc77183492eb2cd8991167872f879b13e483f2e986a5d3869fae",
+    "source-header.sam.gz": "a27d48b1da3e3f299f3139bf2fc14f0cf551841271a3c974b8e1c346ede91299"
+  }
+}