"""Matched linear-haplotype control for the graph-assisted inference model. Map the identical observed FASTA paths with minimap2, retain paired alternative placements, combine unchanged native IPD candidates, and apply the same scorer. The graph-derived path identifiers describe sequence provenance, not graph use. """ import argparse import collections from contextlib import closing import itertools import json import os from pathlib import Path import shutil import sqlite3 import subprocess import time from build_graph import sha from fragment_support import pair_placements from join_evidence import merge from t1k_evidence import parse from run_context_evidence_disk_v3 import run as evidence from run_context_inference import run as inference def fragments(connection): cursor = connection.execute('SELECT fragment,path,mate,start,end,reverse,score FROM alignments ORDER BY fragment,path,mate') for fragment, records in itertools.groupby(cursor, key=lambda r: r[0]): candidates = {} for path, rows in itertools.groupby(records, key=lambda r: r[1]): mates = collections.defaultdict(list) for _, _, mate, start, end, reverse, score in rows: mates[mate].append(dict(score=score, placements=[dict(path=path, start=start, end=end, reverse=bool(reverse))])) candidates.update(pair_placements(mates[1], mates[2], 1000)) yield dict(read_names=[fragment + '/1', fragment + '/2'], status='supported' if candidates else 'no_concordant_path', candidates=candidates) def alignment_limit_status(connection, secondary_limit): """Check raw emitted records, not deduplicated placements, for truncation.""" if isinstance(secondary_limit, bool) or not isinstance(secondary_limit, int) or secondary_limit < 1: raise ValueError('Secondary limit must be a positive integer') maximum = connection.execute('SELECT MAX(n) FROM (SELECT COUNT(*) AS n FROM alignments GROUP BY fragment,mate)').fetchone()[0] or 0 return maximum, maximum >= secondary_limit + 1 def run(panel, native_evidence, native_reference, baseline, calibration, reads, output, threads, secondary_limit=10000): if isinstance(secondary_limit, bool) or not isinstance(secondary_limit, int) or secondary_limit < 1: raise ValueError('Secondary limit must be a positive integer') if output.exists(): raise FileExistsError(output) if threads < 1 or threads > int(os.environ.get('SLURM_CPUS_PER_TASK', '0')): raise RuntimeError('Requires sufficient Slurm CPUs') import pysam mapper = shutil.which('minimap2') if not mapper: raise FileNotFoundError('minimap2') output.mkdir(parents=True) temporary = output / 'tmp' temporary.mkdir() os.environ['TMPDIR'] = str(temporary.resolve()) os.environ['SQLITE_TMPDIR'] = str(temporary.resolve()) mapping = output / 'linear-mapping' mapping.mkdir() fasta = mapping / 'paths.fa' with fasta.open('wb') as target: for source in sorted(panel.glob('HLA-*.fa')): target.write(source.read_bytes()) read_hashes = {str(p): sha(p) for p in reads} record = dict(status='running', started=time.time(), mapper='minimap2', mapper_sha256=sha(Path(mapper)), mapper_version=subprocess.check_output([mapper, '--version'], text=True).strip(), reads_sha256=read_hashes, reference_sha256=sha(fasta), driver_sha256=sha(Path(__file__)), scope='Linear mapping to identical observed haplotype sequences; no graph alignment', secondary_limit=secondary_limit) command = [mapper, '-ax', 'sr', '--secondary=yes', '-N', str(secondary_limit), '-t', str(threads), str(fasta), *map(str, reads)] record['command'] = command def save(): (mapping / 'manifest.json').write_text(json.dumps(record, indent=2) + '\n') save() try: counts = collections.Counter() with closing(sqlite3.connect(mapping / 'placements.sqlite')) as db, (mapping / 'mapping.log').open('wb') as log: db.execute('PRAGMA temp_store=FILE') db.execute('PRAGMA cache_size=-16384') db.execute('CREATE TABLE alignments(fragment TEXT,path TEXT,mate INTEGER,start INTEGER,end INTEGER,reverse INTEGER,score REAL)') process = subprocess.Popen(command, stdout=subprocess.PIPE, stderr=log) try: with pysam.AlignmentFile(process.stdout, 'r') as source, pysam.AlignmentFile(str(mapping / 'mapping.bam'), 'wb', template=source) as bam: batch = [] for alignment in source: bam.write(alignment) counts['alignment_records'] += 1 if alignment.is_unmapped or alignment.is_supplementary: continue if not alignment.is_paired or alignment.is_read1 == alignment.is_read2: raise ValueError('Linear mapper lost mate identity') name = alignment.query_name if name.endswith(('/1', '/2')): name = name[:-2] batch.append((name, alignment.reference_name, 1 if alignment.is_read1 else 2, alignment.reference_start, alignment.reference_end, int(alignment.is_reverse), float(alignment.get_tag('AS')))) if len(batch) >= 5000: db.executemany('INSERT INTO alignments VALUES (?,?,?,?,?,?,?)', batch) batch.clear() if batch: db.executemany('INSERT INTO alignments VALUES (?,?,?,?,?,?,?)', batch) if process.wait() != 0: raise RuntimeError('Linear mapper failed; see mapping.log') except BaseException: process.kill() process.wait() raise db.commit() db.execute('CREATE INDEX ordered ON alignments(fragment,path,mate)') max_alignments, saturated = alignment_limit_status(db, secondary_limit) record['maximum_alignments_per_mate'] = max_alignments save() if saturated: raise ValueError('Linear alternative cap reached; raise cap before making a matched-control claim') with (mapping / 'fragments.jsonl').open('w') as target: for row in fragments(db): target.write(json.dumps(row, separators=(',', ':')) + '\n') counts['fragments_with_placements'] += 1 counts['supported_fragments'] += bool(row['candidates']) record.update(status='complete', finished=time.time(), counts=dict(counts), output_sha256={n: sha(mapping / n) for n in ['paths.fa', 'mapping.bam', 'fragments.jsonl']}) save() (mapping / 'COMPLETE.json').write_text(json.dumps(record, indent=2) + '\n') # Match fragment identity once across native and linear routes. This uses # the common join row schema; its historical graph_sources field contains # linear path placements in this explicitly labelled control. native_manifest = json.loads((native_evidence / 'COMPLETE.json').read_text()) native_files = list(native_evidence.glob('*_assign.tsv')) if len(native_files) != 1 or native_manifest['status'] != 'complete': raise ValueError('Invalid native export') native_file = native_files[0] if sha(native_file) != native_manifest['output_sha256'][native_file.name]: raise ValueError('Native assignments changed') original_reads = sorted(h for n, h in native_manifest['input_sha256'].items() if n.endswith(('.fq', '.fq.gz', '.fastq', '.fastq.gz'))) if original_reads != sorted(read_hashes.values()): raise ValueError('Native/linear read mismatch') joined = output / 'joined' joined.mkdir() with closing(sqlite3.connect(joined / 'join.sqlite')) as db, native_file.open() as ns, (mapping / 'fragments.jsonl').open() as ls, (joined / 'evidence.jsonl').open('w') as target: for row in merge(db, parse(ns), [map(json.loads, ls)]): target.write(json.dumps(row, separators=(',', ':')) + '\n') (joined / 'COMPLETE.json').write_text(json.dumps(dict(status='complete', native_manifest_sha256=sha(native_evidence / 'COMPLETE.json'), graph_chains=[{str(mapping): sha(mapping / 'COMPLETE.json')}], candidate_source='linear_haplotype_control', output_sha256=sha(joined / 'evidence.jsonl'), driver_sha256=sha(Path(__file__))), indent=2) + '\n') evidence(joined, native_reference, panel, calibration, reads, output / 'evidence') inference(output / 'evidence', baseline, output / 'inference', 'E3') (output / 'COMPLETE.json').write_text(json.dumps(dict(status='complete', experiment='linear-matched', mapping_sha256=sha(mapping / 'COMPLETE.json'), evidence_sha256=sha(output / 'evidence/COMPLETE.json'), inference_sha256=sha(output / 'inference/COMPLETE.json'), driver_sha256=sha(Path(__file__))), indent=2) + '\n') except BaseException as error: (output / 'FAILED.json').write_text(json.dumps(dict(status='failed', error=repr(error)), indent=2) + '\n') raise if __name__ == '__main__': p = argparse.ArgumentParser(description=__doc__) for name in ['panel', 'native-evidence', 'native-reference', 'baseline', 'calibration', 'output']: p.add_argument('--' + name, type=Path, required=True) p.add_argument('--reads', type=Path, nargs=2, required=True) p.add_argument('--threads', type=int, default=4) p.add_argument('--secondary-limit', type=int, default=10000) a = p.parse_args() run(a.panel, a.native_evidence, a.native_reference, a.baseline, a.calibration, a.reads, a.output, a.threads, a.secondary_limit)