import gzip,csv,math,statistics,json,hashlib,zipfile,xml.etree.ElementTree as ET
from collections import defaultdict
from pathlib import Path

D=Path(__file__).resolve().parent
H=D/'GSE325735_human_raw_counts.txt.gz'
X=D/'43587_2026_1234_MOESM15_ESM.xlsx'
OUT=D/'public_summary.json'
NS={'m':'http://schemas.openxmlformats.org/spreadsheetml/2006/main','r':'http://schemas.openxmlformats.org/officeDocument/2006/relationships'}

def sha(p):
 h=hashlib.sha256()
 with open(p,'rb') as f:
  for b in iter(lambda:f.read(1<<20),b''):h.update(b)
 return h.hexdigest()

def read_sheet(path,sheet_index=1):
 z=zipfile.ZipFile(path); ss=[]
 if 'xl/sharedStrings.xml' in z.namelist():
  rt=ET.fromstring(z.read('xl/sharedStrings.xml'))
  ss=[''.join(t.text or '' for t in x.findall('.//m:t',NS)) for x in rt.findall('m:si',NS)]
 wb=ET.fromstring(z.read('xl/workbook.xml'))
 rel=ET.fromstring(z.read('xl/_rels/workbook.xml.rels'))
 rels={x.attrib['Id']:x.attrib['Target'] for x in rel}
 sh=wb.findall('.//m:sheet',NS)[sheet_index-1]
 target=rels[sh.attrib['{'+NS['r']+'}id']]
 if not target.startswith('/'): target='xl/'+target
 root=ET.fromstring(z.read(target)); rows=[]
 for row in root.findall('.//m:sheetData/m:row',NS):
  vals=[]
  for c in row.findall('m:c',NS):
   v=c.find('m:v',NS)
   if v is None: val=''
   elif c.attrib.get('t')=='s': val=ss[int(v.text)]
   else: val=v.text
   vals.append(val)
  rows.append(vals)
 return sh.attrib['name'],rows

# Supplementary Table 15, sheet 3a: direct SampleID to treatment/donor metadata.
sheet,rows=read_sheet(X,1)
header=rows[0]
ix={name:i for i,name in enumerate(header)}
required=['SampleID','Replicate','Gene','Line','Passage','Line_age','Dox','Line_coded']
if not all(k in ix for k in required): raise SystemExit('Unexpected supplementary table schema')
with gzip.open(H,'rt') as f:
 rd=csv.reader(f,delimiter='\t')
 count_header=next(rd)
count_aliases=count_header[1:]
if len(set(count_aliases))!=len(count_aliases): raise SystemExit('Duplicate count matrix aliases')
alias_index={a:i for i,a in enumerate(count_aliases)}
selected=[]
for row in rows[1:]:
 if len(row)<=max(ix.values()): continue
 if row[ix['Dox']]!='2000' or row[ix['Gene']] not in ('SRSF1','BFP'): continue
 sid=row[ix['SampleID']]
 if sid in alias_index:
  selected.append({k:row[ix[k]] for k in required})
# The expected comparison has 6 donors x 2 replicates x 2 groups = 24 samples.
if len(selected)!=24: raise SystemExit('Expected 24 directly mapped SRSF1/BFP samples, found '+str(len(selected)))
groups=defaultdict(lambda:defaultdict(list))
for r in selected:
 key=(r['Line_coded'],r['Line'],r['Passage'],r['Dox'])
 groups[key][r['Gene']].append(r)
if len(groups)!=6: raise SystemExit('Expected 6 independent donor strata, found '+str(len(groups)))
for k,v in groups.items():
 if len(v.get('SRSF1',[]))!=2 or len(v.get('BFP',[]))!=2:
  raise SystemExit('Unbalanced replicate group: '+repr((k,{g:len(x) for g,x in v.items()})))

# A priori targeted panel: SRSF1 manipulation, inflammatory/oncogenic context,
# senescence/cell-cycle, ECM and stress response. Expression is RNA-level only.
genes={
'SRSF1':'ENSG00000136450','IL1R1':'ENSG00000115594','IL1B':'ENSG00000125538',
'IL6':'ENSG00000136244','CXCL8':'ENSG00000169429','MYC':'ENSG00000136997',
'KRAS':'ENSG00000133703','MAPK1':'ENSG00000100030','MAPK3':'ENSG00000102882',
'NFKBIA':'ENSG00000100906','CDKN1A':'ENSG00000124762','CDKN2A':'ENSG00000147889',
'SERPINE1':'ENSG00000106366','MKI67':'ENSG00000148773','CCND1':'ENSG00000110092',
'TP53':'ENSG00000141510','COL1A1':'ENSG00000108821','COL3A1':'ENSG00000168542'}
id_to_symbol={v:k for k,v in genes.items()}
tot=[0]*len(count_aliases); counts={}
with gzip.open(H,'rt') as f:
 rd=csv.reader(f,delimiter='\t'); next(rd)
 for row in rd:
  if len(row)!=len(count_aliases)+1: continue
  vals=[int(x or 0) for x in row[1:]]
  for i,v in enumerate(vals): tot[i]+=v
  gid=row[0].split('.')[0]
  if gid in id_to_symbol: counts[id_to_symbol[gid]]=vals
missing=sorted(set(genes)-set(counts))
if missing: raise SystemExit('Missing gene IDs: '+','.join(missing))
logcpm={g:[math.log2((v/t*1e6)+0.5) if t else float('nan') for v,t in zip(vals,tot)] for g,vals in counts.items()}

effects_by_gene=defaultdict(list); donor_rows=[]
for (coded,line,passage,dox),v in sorted(groups.items()):
 rowout={'donor_code':coded,'donor_line':line,'passage':int(passage),'dox_ng_ml':int(dox),
         'SRSF1_sample_ids':[x['SampleID'] for x in v['SRSF1']],
         'BFP_control_sample_ids':[x['SampleID'] for x in v['BFP']]}
 diffs={}
 for g,vec in logcpm.items():
  s=statistics.mean(vec[alias_index[x['SampleID']]] for x in v['SRSF1'])
  c=statistics.mean(vec[alias_index[x['SampleID']]] for x in v['BFP'])
  diffs[g]=s-c
  effects_by_gene[g].append((coded,diffs[g]))
 rowout['log2_CPM_plus_half_differences']=diffs
 donor_rows.append(rowout)

summary={}
for g,items in effects_by_gene.items():
 vals=[x[1] for x in items]
 pos=sum(x>0 for x in vals); neg=sum(x<0 for x in vals); nz=pos+neg
 k=min(pos,neg)
 p=1.0 if nz==0 else min(1.0,2*sum(math.comb(nz,j) for j in range(k+1))/(2**nz))
 mean=statistics.mean(vals)
 summary[g]={
  'n_independent_donor_lines':len(vals),
  'mean_log2_CPM_difference':mean,
  'median_log2_CPM_difference':statistics.median(vals),
  'fold_change_CPM_plus_half':2**mean,
  'donor_lines_increased':pos,'donor_lines_decreased':neg,
  'exploratory_two_sided_exact_sign_test_p':p,
  'directional_summary_only':True
 }
result={
 'dataset':'GSE325735',
 'contrast':'Dox 2000 ng/mL, SRSF1 transgene versus mTagBFP2/BFP control, matched by donor line and passage',
 'mapping_source':'Direct sample-ID mapping from Nature Aging Supplementary Table 15, sheet 3a (MOESM15), whose SampleID values match the human raw-count matrix column aliases.',
 'mapping_qc':{
  'supplementary_table_sample_rows_Dox2000_SRSF1_or_BFP':len([r for r in rows[1:] if len(r)>ix['Dox'] and r[ix['Dox']]=='2000' and r[ix['Gene']] in ('SRSF1','BFP')]),
  'directly_mapped_comparison_samples':len(selected),
  'comparison_groups':len(groups),
  'independent_donor_line_count':len(groups),
  'replicates_per_group':2,
  'all_24_comparison_samples_directly_mapped':True
 },
 'method':'Raw gene counts; log2(CPM+0.5); two technical replicates averaged within each group; paired SRSF1-minus-BFP contrasts computed per independent donor line. No mapping by GEO accession order.',
 'human_raw_counts_sha256':sha(H),
 'supplementary_table15_sha256':sha(X),
 'gene_ids':genes,
 
 'summary':summary,
 'limitations':[
  'Secondary reanalysis of the originating study, not independent biological replication.',
  'Only six donor-derived fibroblast lines, one Dox dose and one 3-day induction schedule are analyzed.',
  'Technical replicates are not treated as independent biological replicates.',
  'Exploratory sign tests have low power at n=6; no p-value is presented as confirmatory evidence.',
  'Selected marker panel only; no full differential-expression, splice-isoform, protein activity, tumor-incidence, organ-toxicity or clinical analysis.',
  'Bulk RNA expression cannot establish human benefit or safety.',
  'No therapeutic use, gene therapy, or self-experimentation is supported.'
 ]
}
OUT.write_text(json.dumps(result,indent=2,ensure_ascii=False)+'\n')
print('VALIDATED_SAMPLE_MAPPING',len(selected),'/',len([r for r in rows[1:] if len(r)>ix['Dox'] and r[ix['Dox']]=='2000' and r[ix['Gene']] in ('SRSF1','BFP')]))
print('DONOR_LINES',sorted(k[0] for k in groups))
for g in ['SRSF1','IL1R1','IL1B','IL6','CXCL8','MYC','KRAS','MAPK1','MAPK3','NFKBIA','CDKN1A','CDKN2A','SERPINE1','MKI67','CCND1','TP53','COL1A1','COL3A1']:
 x=summary[g]
 print(g,'FC',round(x['fold_change_CPM_plus_half'],3),'L+',x['donor_lines_increased'],'L-',x['donor_lines_decreased'],'signP',round(x['exploratory_two_sided_exact_sign_test_p'],4))
print('JSON_BYTES',OUT.stat().st_size,'SHA256',sha(OUT))