"""Portable reconstruction adapted from accepted matrix_correction_v4.py. Inputs: model/TK1151_uniform_20260930.gtf, gene_identity_map.tsv, sample_order.tsv, gtf/RUN.gtf and abundance/RUN.gene_abundance.tsv. The original source data and per-run StringTie outputs are never changed. """ import argparse import collections import csv import gzip import hashlib import math import re from array import array from pathlib import Path ATTR = re.compile(r'(\w+) "([^"]*)"') METRICS = ('Coverage', 'FPKM', 'TPM') MODEL_SHA = 'b63dabd24e9d3cd96fb574a1adfe0d0f64c385da0eb3d272d012f299694b2245' def lines(path): return gzip.open(path, 'rt', encoding='utf-8') if path.suffix == '.gz' else path.open(encoding='utf-8') def sha(path): digest = hashlib.sha256() with (gzip.open(path, 'rb') if path.suffix == '.gz' else path.open('rb')) as handle: for block in iter(lambda: handle.read(1024 * 1024), b''): digest.update(block) return digest.hexdigest() def table(path): with path.open(encoding='utf-8', newline='') as handle: return list(csv.DictReader(handle, delimiter='\t')) def reconstruct_sample(sid, inputs, guide, tx_gene, exons, tx_lengths, union_lengths): transcripts = {} with lines(inputs / 'gtf' / (sid + '.gtf')) as handle: for line in handle: if line.startswith('#'): continue row = line.rstrip().split('\t') if row[2] != 'transcript': continue attributes = dict(ATTR.findall(row[8])) tid, gid = attributes['transcript_id'], attributes['gene_id'] assert tx_gene[tid] == gid values = tuple(float(attributes[key]) for key in ('cov', 'FPKM', 'TPM')) if tid in transcripts: previous = transcripts[tid] assert all(x == 0 or y == 0 or x == y for x, y in zip(previous, values)) values = tuple(max(x, y) for x, y in zip(previous, values)) transcripts[tid] = values assert set(transcripts) == set(tx_gene) contributions = collections.defaultdict(list) for tid, values in transcripts.items(): contributions[tx_gene[tid]].append((values, tid)) original = {} for row in table(inputs / 'abundance' / (sid + '.gene_abundance.tsv')): gid = row['Gene ID'] values = tuple(float(row[key]) for key in METRICS) original[gid] = tuple(max(x, y) for x, y in zip(original.get(gid, (0, 0, 0)), values)) assert set(original) == set(guide) corrected = {} for gid, items in contributions.items(): fpkm = math.fsum(values[1] for values, tid in items) tpm = math.fsum(values[2] for values, tid in items) significant_zero = original[gid] == (0, 0, 0) and max(fpkm, tpm) > 2e-6 * (1 + len(items)) if significant_zero: assert guide[gid]['strand'] == '.' assert len(items) == 1 and len(exons[items[0][1]]) == 1 assert tx_lengths[items[0][1]] == union_lengths[gid] coverage = items[0][0][0] else: coverage = original[gid][0] corrected[gid] = (coverage, fpkm, tpm) assert all(math.isfinite(value) and value >= 0 for value in corrected[gid]) assert set(corrected) == set(guide) assert abs(math.fsum(values[2] for values in corrected.values()) - 1000000) < 1 return corrected def main(): parser = argparse.ArgumentParser() parser.add_argument('--inputs', type=Path, required=True) parser.add_argument('--output-dir', type=Path, required=True) parser.add_argument('--model', type=Path) parser.add_argument('--gene-map', type=Path) parser.add_argument('--sample-order', type=Path) parser.add_argument('--samples', nargs='+', help='Optional subset for validation; not a full-release export.') parser.add_argument('--reference-dir', type=Path, help='Compare every generated subset cell to accepted full matrices.') args = parser.parse_args() model = args.model or args.inputs / 'model/TK1151_uniform_20260930.gtf' guide_path = args.gene_map or args.inputs / 'model/gene_identity_map.tsv' order_path = args.sample_order or args.inputs / 'samples.tsv' assert sha(model) == MODEL_SHA, 'Use the exact released fixed model.' guide = {row['uniform_gene_id']: row for row in table(guide_path)} genes = list(guide) order = [row['sample_id'] for row in table(order_path)] samples = args.samples or order assert len(set(order)) == len(order) == 109 and set(samples) <= set(order) assert len(set(samples)) == len(samples) tx_gene, exons, gene_exons = {}, collections.defaultdict(set), collections.defaultdict(set) with lines(model) as handle: for line in handle: if line.startswith('#'): continue row = line.rstrip().split('\t') attributes = dict(ATTR.findall(row[8])) tid, gid = attributes['transcript_id'], attributes['gene_id'] tx_gene[tid] = gid if row[2] == 'exon': interval = (int(row[3]), int(row[4])) exons[tid].add(interval) gene_exons[gid].add(interval) assert len(genes) == 98473 and len(tx_gene) == 257009 tx_lengths = {tid: sum(end - start + 1 for start, end in spans) for tid, spans in exons.items()} union_lengths = {} for gid, spans in gene_exons.items(): merged = [] for start, end in sorted(spans): if merged and start <= merged[-1][1] + 1: merged[-1][1] = max(merged[-1][1], end) else: merged.append([start, end]) union_lengths[gid] = sum(end - start + 1 for start, end in merged) columns = {metric: [] for metric in METRICS} for sid in samples: corrected = reconstruct_sample(sid, args.inputs, guide, tx_gene, exons, tx_lengths, union_lengths) for index, metric in enumerate(METRICS): columns[metric].append(array('d', (corrected[gid][index] for gid in genes))) print('Reconstructed', sid, len(corrected), 'loci', flush=True) args.output_dir.mkdir(parents=True, exist_ok=True) for metric in METRICS: filename = 'TK1151_uniform_' + metric.lower() + '_matrix.tsv' output = args.output_dir / filename if output.exists(): raise FileExistsError('Use an empty output directory; refusing to overwrite ' + str(output)) with output.open('w', encoding='utf-8', newline='') as handle: writer = csv.writer(handle, delimiter='\t') writer.writerow(['uniform_gene_id'] + samples) for index, gid in enumerate(genes): writer.writerow([gid] + [format(column[index], '.10g') for column in columns[metric]]) if args.reference_dir: with (args.reference_dir / filename).open(encoding='utf-8', newline='') as handle, output.open(encoding='utf-8', newline='') as generated: reference = csv.reader(handle, delimiter='\t') subset = csv.reader(generated, delimiter='\t') header = next(reference) subset_header = next(subset) indices = [header.index(sid) for sid in samples] assert subset_header[1:] == samples count = 0 for row in reference: actual = next(subset) assert actual == [row[0]] + [row[index] for index in indices], (metric, row[0]) count += 1 assert next(subset, None) is None and count == 98473 print('All reference subset cells match', metric, count * len(samples), flush=True) print(filename, sha(output), flush=True) if __name__ == '__main__': main()