#!/usr/bin/env python3
"""Reproduce Predome's public-data case records. Python standard library only.

This computes descriptive summaries of retained public data. It does not call an
LLM, repeat a laboratory experiment or establish biological/clinical validity.
Pass --data-dir when running the downloaded reproduction bundle standalone.
Pass --refresh-pubchem to explicitly fetch the official reference properties.
"""
from pathlib import Path
import argparse, csv, datetime, gzip, hashlib, io, json, math, shutil, tarfile, time, urllib.request, zipfile

ROOT = Path(__file__).resolve().parents[1]
PUBCHEM = 'https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/cid/3385/property/MolecularFormula,MolecularWeight,CanonicalSMILES/JSON'
ARCHIVE_HASH = 'e4c78d89055b94cf18730e098aa8326ad5a7ae8468b0bb3392102d4085eb1459'
def utc(): return datetime.datetime.now(datetime.timezone.utc).isoformat(timespec='milliseconds').replace('+00:00','Z')
def digest(value):
    data = value if isinstance(value,bytes) else json.dumps(value,sort_keys=True,separators=(',',':')).encode()
    return hashlib.sha256(data).hexdigest()
def write(path,value): path.write_text(json.dumps(value,indent=2,ensure_ascii=False)+'\n')

def main():
    parser=argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--data-dir',type=Path)
    parser.add_argument('--refresh-pubchem',action='store_true')
    args=parser.parse_args()
    dest=args.data_dir.resolve() if args.data_dir else ROOT/'public/case-data'
    dest.mkdir(parents=True,exist_ok=True)
    started=utc();clock=time.perf_counter();trace=[]
    def stage(name,action):
        at=utc();begin=time.perf_counter();output=action()
        trace.append({'name':name,'engine':'deterministic Python','status':'completed','started_at':at,'completed_at':utc(),'elapsed_ms':round((time.perf_counter()-begin)*1000,3),'output':output})
        return output
    archive=(dest/'GSE90944_RAW.tar').read_bytes()
    if digest(archive)!=ARCHIVE_HASH: raise ValueError('Original GEO archive fingerprint changed; review before recomputing.')
    selected=json.loads((dest/'gse90944-selected-genes.json').read_text())
    tables={}
    def load_source():
        with tarfile.open(fileobj=io.BytesIO(archive)) as tar:
            for member in tar.getmembers():
                if not member.isfile() or not member.name.endswith('.count.txt.gz'):continue
                accession=member.name.split('_')[0]
                text=gzip.decompress(tar.extractfile(member).read()).decode()
                rows=[line.split('\t') for line in text.strip().splitlines()]
                if any(len(row)!=2 for row in rows):raise ValueError('Malformed original count table')
                values={gene:int(count) for gene,count in rows}
                if len(values)!=len(rows) or any(v<0 for v in values.values()):raise ValueError('Duplicate or negative source counts')
                tables[accession]=values
        if set(tables)!={s['sample_id'] for s in selected['samples']}:raise ValueError('Unexpected source sample set')
        return {'archive_sha256':digest(archive),'tables':len(tables),'table_rows':{k:len(v) for k,v in tables.items()},'source_date':'2017','source_kind':'submitter-provided processed HTSeq counts; not FASTQ'}
    stage('Source verification',load_source)
    samples=selected['samples'];gene_rows=[];measurements=[]
    def validate():
        matched=0
        for sample in samples:
            ident=sample['sample_id'];library=sum(v for gene,v in tables[ident].items() if not gene.startswith('__'))
            if library!=sample['library_size']:raise ValueError('Library denominator differs from complete source table')
            metadata=(dest/f'{ident}.soft.txt').read_text()
            expected='FN Treatment' if sample['group']=='exposed' else 'NC,'
            if expected not in metadata:raise ValueError('Sample group differs from source metadata')
            for gene in selected['genes']:
                row=next(x for x in gene['measurements'] if x['sample_id']==ident)
                count=tables[ident][gene['symbol']]
                if count!=row['raw_count'] or abs(row['cpm']-count/library*1e6)>1e-8:raise ValueError('Selected count or CPM differs from source')
                matched+=1
        return {'matched_counts':matched,'matched_library_denominators':len(samples),'selected_genes':len(selected['genes']),'missing_values':0,'negative_counts':0,'biological_independence':'not established by source metadata'}
    stage('Numeric and provenance checks',validate)
    def calculate():
        for gene in selected['genes']:
            points=[]
            for sample in samples:
                count=tables[sample['sample_id']][gene['symbol']]
                point={'gene':gene['symbol'],'sample_id':sample['sample_id'],'sample_label':sample['label'],'group':sample['group'],'raw_count':count,'library_size':sample['library_size'],'cpm':round(count/sample['library_size']*1e6,9)}
                points.append(point);measurements.append(point)
            mean=lambda group:sum(p['cpm'] for p in points if p['group']==group)/3
            control=mean('control');exposed=mean('exposed')
            if control<=0 or exposed<=0:raise ValueError('A zero mean requires an explicitly revised ratio method')
            gene_rows.append({'symbol':gene['symbol'],'pathway':gene['pathway'],'control_mean_cpm':round(control,9),'exposed_mean_cpm':round(exposed,9),'ratio':round(exposed/control,9),'log2_ratio':round(math.log2(exposed/control),9),'samples':points})
        return {'formula':'CPM=count/sum(non-summary gene counts)*1000000; ratio=mean(exposed CPM)/mean(control CPM)','samples':len(samples),'counts':len(measurements),'derived_values':len(measurements),'genes':len(gene_rows),'no_pseudocount':True,'inferential_test':'not performed','result_sha256':digest(gene_rows)}
    stage('Descriptive calculation',calculate)
    # Preserve the retained input JSON byte-for-byte; adapt only this published
    # report's navigation links to the static website's real asset locations.
    published_dataset={**selected['dataset'],
      'download_csv':'/case-data/gse90944-selected-genes.csv',
      'download_json':'/case-data/gse90944-selected-genes.json',
      'download_archive':'/case-data/GSE90944_RAW.tar'}
    report={'case_id':'host-microbe','record_kind':'computational_reproduction','new_wet_lab_experiment':False,'llm_used':False,'question':'How do ten prespecified transcripts compare descriptively between F. nucleatum-exposed and medium-control HT-29 samples in GSE90944?','dataset':published_dataset,'samples':samples,'results':gene_rows,'trace':trace,'input_sha256':{'original_archive':digest(archive),'selected_json':digest((dest/'gse90944-selected-genes.json').read_bytes())},'result_sha256':digest(gene_rows),'run':{'started_at':started,'completed_at':utc(),'elapsed_ms':round((time.perf_counter()-clock)*1000,3)},'interpretation':'Transcript abundance differences in this selected public-data subset are descriptive. They do not establish autophagic activity, a causal mechanism, chemotherapy resistance or clinical efficacy.','next_validation':{'status':'proposed_not_executed','tasks':['Review sample independence and batch metadata.','Repeat a prespecified analysis on an independent dataset.','Use an appropriate experimental assay before claiming protein function or drug response.']}}
    write(dest/'host-microbe-report.json',report)
    with (dest/'host-microbe-measurements.csv').open('w',newline='') as f:
        writer=csv.DictWriter(f,fieldnames=list(measurements[0]));writer.writeheader();writer.writerows(measurements)
    with (dest/'host-microbe-summary.csv').open('w',newline='') as f:
        keys=['symbol','control_mean_cpm','exposed_mean_cpm','ratio','log2_ratio']
        writer=csv.DictWriter(f,fieldnames=keys);writer.writeheader();writer.writerows({k:row[k] for k in keys} for row in gene_rows)
    molecular_path=dest/'pubchem-3385.json';pubchem_error=None
    if args.refresh_pubchem:
        try:
            request=urllib.request.Request(PUBCHEM,headers={'User-Agent':'Predome-public-reference/1.0','Accept':'application/json'})
            with urllib.request.urlopen(request,timeout=20) as response:
                raw=response.read(100000)
                if response.status!=200:raise ValueError('PubChem returned a non-success response')
            parsed=json.loads(raw);prop=parsed['PropertyTable']['Properties'][0]
            if prop['CID']!=3385 or not isinstance(prop.get('MolecularFormula'),str) or float(prop['MolecularWeight'])<=0:raise ValueError('Unexpected PubChem response schema')
            molecular_path.write_bytes(raw)
            write(dest/'pubchem-retrieval.json',{'source_url':PUBCHEM,'collected_at':utc(),'sha256':digest(raw)})
        except Exception as error:pubchem_error=type(error).__name__+': '+str(error)
    molecular={'case_id':'molecular-reference','record_kind':'official_reference_audit','question':'Can a 5-fluorouracil reference profile preserve the compound identity, properties and source?','source_url':'https://pubchem.ncbi.nlm.nih.gov/compound/3385','api_url':PUBCHEM,'status':'source_not_retrieved','llm_used':False,'drug_effect_tested':False,'properties':None,'retrieval_error':pubchem_error,'notice':'A compound reference profile is not a new molecule, a drug-response prediction or evidence of efficacy.'}
    if molecular_path.exists():
        raw=molecular_path.read_bytes();props=json.loads(raw)['PropertyTable']['Properties'][0]
        if props['CID']!=3385:raise ValueError('Wrong PubChem compound')
        molecular.update({'status':'reference_verified','properties':props,'sha256':digest(raw),'retrieval':json.loads((dest/'pubchem-retrieval.json').read_text()),'checks':{'matching_cid':True,'molecular_weight_positive':float(props['MolecularWeight'])>0,'formula_present':bool(props['MolecularFormula'])}})
    write(dest/'molecular-reference-report.json',molecular)
    registry_path=ROOT/'public/data/registry-snapshot.json'
    if args.data_dir:registry_path=dest/'registry-snapshot.json'
    registry=json.loads(registry_path.read_text());records=registry['registry']
    clinical={'case_id':'clinical-question','record_kind':'registry_metadata_audit','question':'What is present, and still missing, before a registry reference becomes a forecastable clinical question?','source_collected_at':registry['collected_at'],'source_url':registry['source_url'],'source_sha256':registry['source_version'],'snapshot_sha256':digest(registry_path.read_bytes()),'audit_executed_at':utc(),'records':records,'summary':{'records':len(records),'recruiting':sum(r['status']=='RECRUITING' for r in records),'phase_missing':sum(not r['phase'] for r in records),'phase_not_applicable':sum('NA' in r['phase'] for r in records),'phase1_mentions':sum('PHASE1' in r['phase'] for r in records),'phase2_mentions':sum('PHASE2' in r['phase'] for r in records)},'model_submissions':[],'registered_forecasts':0,'required_before_forecast':['Confirm no relevant result is already public.','Review the exact endpoint, population, comparator and protocol.','Set the evidence cutoff and an authoritative result source.','Register resolution, unresolved and void rules before any forecast.'],'notice':'Metadata completeness is checked. Recruitment status, enrollment and update date are not clinical outcomes, upcoming readout dates or evidence of treatment efficacy.'}
    write(dest/'clinical-question-report.json',clinical)
    shutil.copy2(registry_path,dest/'registry-snapshot.json') if registry_path.resolve()!=(dest/'registry-snapshot.json').resolve() else None
    bundle={'schema_version':1,'hostMicrobe':report,'molecular':molecular,'clinical':clinical}
    write(dest/'case-records.json',bundle)
    if not args.data_dir:write(ROOT/'src/lib/case-data.json',bundle)
    here=Path(__file__).resolve()
    if here!=(dest/'build-case-evidence.py').resolve():shutil.copy2(here,dest/'build-case-evidence.py')
    (dest/'REPRODUCE.txt').write_text('Predome public-data reproduction\n\nRun: python3 build-case-evidence.py --data-dir .\nPython 3 standard library only. No model key is required.\nThe default run is offline; --refresh-pubchem explicitly requests the official reference again.\nRun timestamps and elapsed time change; source bytes and result_sha256 should remain the same.\nThese are 2017 public processed RNA-seq counts from the original investigators, not a Predome wet-lab experiment.\nThe displayed ten-gene selection is not a genome-wide screen; no significance, causality or efficacy is established.\n')
    members=['GSE90944_RAW.tar','GSE90944.soft.txt','gse90944-selected-genes.json','gse90944-selected-genes.csv','registry-snapshot.json','build-case-evidence.py','REPRODUCE.txt']+[f"{s['sample_id']}.soft.txt" for s in samples]
    if molecular_path.exists():members+=['pubchem-3385.json','pubchem-retrieval.json']
    with zipfile.ZipFile(dest/'reproduce-public-cases.zip','w',zipfile.ZIP_DEFLATED) as z:
        for name in members:z.write(dest/name,name)
    print(json.dumps({'case_records':3,'matched_original_counts':len(measurements),'library_denominators':len(samples),'result_sha256':report['result_sha256'],'run':report['run'],'pubchem_status':molecular['status'],'registry_references':len(records)}))

if __name__=='__main__':main()