prepare_data.py

from pathlib import Path
import sqlite3,json,csv,gzip,hashlib
from collections import Counter
ROOT=Path(__file__).parent; D=ROOT/'data'
con=sqlite3.connect(D/'taxonomy.sqlite')
taxa={r[0]:r for r in con.execute('SELECT taxid,parent_id,name,rank FROM taxa')}
assemblies=dict(con.execute('SELECT accession,taxid FROM assemblies'))
ref={r['assembly_accession'] for r in csv.DictReader(open(D/'refseq_assembly_taxids.csv'))}
common=json.loads((D/'common_names.json').read_text())['names']
metadata={}
source=Path('/home/node/local/git/genbank-annotation-explorer/.cache/coverage/assembly_summary_genbank.txt')
with source.open() as f:
 for line in f:
  if line.startswith('##'): continue
  if line.startswith('#'): header=line[1:].strip().split('\t'); continue
  a=line.split('\t',1)[0]
  if a not in assemblies: continue
  r=dict(zip(header,line.rstrip('\n').split('\t')))
  def num(k):
   try:return int(r.get(k,''))
   except:return None
  annotation=bool(r.get('annotation_provider') not in ('na','',None))
  metadata[a]={'organism':r['organism_name'],'bases':num('genome_size'),'paired':r['gbrs_paired_asm'] if r['gbrs_paired_asm']!='na' else '', 'pair_match':r['paired_asm_comp'],'prior':int(annotation),'ref':int(r['gbrs_paired_asm'] in ref),'provider':r['annotation_provider'] if annotation else '', 'annotation_date':r['annotation_date'] if annotation else '', 'prior_genes':num('protein_coding_gene_count') if annotation else None}
print('metadata',len(metadata),flush=True)
# Keep the real taxonomy topology; omit taxa without displayed descendants.
used=set()
for t in assemblies.values():
 seen=set()
 while t in taxa and t not in seen:
  seen.add(t);used.add(t)
  if t==2759:break
  t=taxa[t][1]
# Sanity: all leaves belong to Eukaryota.
assert 1 not in used
rows=[]
for a,t in sorted(assemblies.items()):
 m=metadata.get(a,{'prior':-1,'ref':0})
 rows.append([a,t,m.get('organism',taxa[t][2]),m.get('bases'),m.get('prior',-1),m.get('ref',0),m.get('paired',''),m.get('provider',''),m.get('prior_genes'),m.get('annotation_date',''),m.get('pair_match','')])
node_rows=[[t,0 if t==2759 else taxa[t][1],taxa[t][2],taxa[t][3],common.get(str(t),'')] for t in sorted(used)]
out={'meta':{'title':'Carbon-A assembly strands — proof of concept','root':2759,'assemblies':len(rows),'taxa':len(set(assemblies.values())),'release':{'assemblies':48167,'taxa':22617,'loci':566340000,'bases':27180000000000},'snapshot':'Explorer taxonomy and covered-assembly mapping, 7 Oct 2026','prior_annotation':'NCBI GenBank assembly summary, cached 7 Oct 2026; named annotation provider required','refseq':'Exact paired GCF accession match to the project’s 2,055-assembly RefSeq comparison corpus; not all of RefSeq','classification_counts':dict(Counter(r[4] for r in rows)),'refseq_matches':sum(r[5] for r in rows),'layout':'NCBI taxonomy; branch lengths are for layout, not evolutionary distance. Each strand is one assembly. Equal strand widths encode assembly counts, not gene counts.'},'node_fields':['taxid','parent','scientific_name','rank','common_name'],'nodes':node_rows,'assembly_fields':['accession','taxid','organism','bases','prior_annotation','refseq_corpus','paired_refseq','annotation_provider','prior_protein_coding_genes','annotation_date','pair_comparison'],'assemblies':rows}
raw=json.dumps(out,separators=(',',':'),ensure_ascii=False).encode()
(D/'tree-data.json').write_bytes(raw)
(D/'tree-data.json.gz').write_bytes(gzip.compress(raw,mtime=0))
print(json.dumps(out['meta'],indent=2)); print('bytes',len(raw),'gzip',(D/'tree-data.json.gz').stat().st_size,'nodes',len(node_rows))