import zipfile,xml.etree.ElementTree as ET,json,hashlib
from pathlib import Path
D=Path(__file__).resolve().parent
book=D/'43587_2026_1234_MOESM15_ESM.xlsx'
analysis=json.loads((D/'public_summary.json').read_text())
ns={'m':'http://schemas.openxmlformats.org/spreadsheetml/2006/main','r':'http://schemas.openxmlformats.org/officeDocument/2006/relationships'}
z=zipfile.ZipFile(book)
ssroot=ET.fromstring(z.read('xl/sharedStrings.xml'))
ss=[''.join(t.text or '' for t in x.findall('.//m:t',ns)) for x in ssroot.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}
sheet=wb.findall('.//m:sheet',ns)[4]
target=rels[sheet.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)
header=rows[0];ix={k:i for i,k in enumerate(header)}
published={}
for r in rows[1:]:
 if len(r)<=max(ix[k] for k in ['Perturbation','Symbol','log2FoldChange','pvalue','padj','is_spliced','Status']): continue
 if r[ix['Perturbation']]!='SRSF1': continue
 sym=r[ix['Symbol']]
 if sym in analysis['gene_ids']:
  lfc=float(r[ix['log2FoldChange']])
  published[sym]={
   'log2FoldChange':lfc,'foldChange':2**lfc,
   'pvalue':float(r[ix['pvalue']]),'adjusted_pvalue':float(r[ix['padj']]),
   'is_spliced':r[ix['is_spliced']]=='1','status':r[ix['Status']]
  }
combined={}
for gene,ours in analysis['summary'].items():
 combined[gene]={
  'our_exploratory_mean_log2_CPM_difference':ours['mean_log2_CPM_difference'],
  'our_exploratory_CPM_plus_half_fold_change':ours['fold_change_CPM_plus_half'],
  'our_donor_lines_increased':ours['donor_lines_increased'],
  'our_donor_lines_decreased':ours['donor_lines_decreased'],
  'authors_supplementary_DESeq2':published.get(gene)
 }
result={
 'source':'Original study supplementary workbook MOESM15, sheet 3e, compared against a targeted secondary CPM analysis of 24 directly mapped raw-count columns',
 'n_independent_donor_lines':6,
 'technical_replicates_per_group':2,
 'dose_ng_ml':2000,'induction_days':3,
 'shared_genes':len(published),
 'genes':combined,
 'interpretation':'This is a methodological cross-check within the same originating study dataset, not independent biological replication. The authors supplementary DESeq2 output strongly corroborates direction and statistical signal for CXCL8 and CDKN1A up, and CDKN2A, SERPINE1, IL1R1, MYC, COL3A1 down. COL1A1 is directionally down but does not pass the reported adjusted-p-value 0.05 threshold (padj ~0.061).',
 'limitations':analysis['limitations']+[
  'The authors DESeq2 table is not independent data; it is the originating study analysis.',
  'The small independent secondary donor-level sign test and the paper DESeq2 model have different statistical assumptions and are not interchangeable.'
 ]
}
p=D/'public_deseq2_crosscheck.json'
p.write_text(json.dumps(result,indent=2,ensure_ascii=False)+'\n')
h=hashlib.sha256(p.read_bytes()).hexdigest()
print('SHARED_GENES',len(published),'OUTPUT_BYTES',p.stat().st_size,'SHA256',h)
for g in ['CXCL8','CDKN1A','CDKN2A','SERPINE1','IL1R1','MYC','COL1A1','COL3A1','IL1B','IL6','SRSF1']:
 x=combined[g]; a=x['authors_supplementary_DESeq2']
 print(g,'ours_FC',round(x['our_exploratory_CPM_plus_half_fold_change'],3),'ours_L+',x['our_donor_lines_increased'],'ours_L-',x['our_donor_lines_decreased'],
       'paper_FC',round(a['foldChange'],3) if a else None,'paper_padj',a['adjusted_pvalue'] if a else None)