GENOMICS-CDDBRG/Genome_Visualization_Tool
0
1"""2app.py - Gulumba Genomics Suite v1.03 4A single Streamlit application combining, on a shared, tested scientific core5(genomics_core.py), the functionality previously spread across five separate6tools: GenomeScope Pro (composition/codon-bias viewer), the Prokka annotation7workbench, the Pan-Genome analysis app, and the Comparative Genomics app.8 9Each tool is its own tab with its own "Run" button and its own inputs, exactly10as the separate apps worked before, so nothing about how you use any one tool11changes. What changes underneath:12 13 - All GC%, codon usage, ENC, CAI, and assembly-statistics calculations are14 computed once, in genomics_core.py, and reused by every tool. The15 original tools each had their own, occasionally disagreeing, copy of16 these formulas (e.g. GC% over unambiguous bases only vs. over full17 length including N). That disagreement is resolved by having one18 implementation.19 - Comparative Genomics gained an alignment-free mode (MinHash/Mash distance20 and ANI, plus a neighbour-joining tree) for genomes too large for the21 original pairwise Needleman-Wunsch approach, which does not scale past a22 few Mb. The original pairwise mode is kept for what it was always good23 at: genes, plasmids, amplicons, and other short, closely related24 sequences.25 - Pan-Genome analysis gained Heaps' law openness estimation and Tettelin-26 style gene-accumulation curves (with permutation resampling), so you can27 report whether the pan-genome is open or closed, not just the core/28 accessory/unique counts.29 - Every figure is publication-styled (600 dpi raster + vector PDF/SVG,30 colour-blind-safe palettes, Type-42 embedded fonts) and every heuristic31 call (origin/terminus prediction, self-derived CAI weights, alignment-32 free ANI) carries an honest confidence/provenance label in the figure or33 table itself, not just in a caption you might not export.34 35Requirements36 Core (always required): streamlit, pandas, numpy, scipy, matplotlib37 GenBank/EMBL parsing: biopython38 Annotation workbench: prokka (external tool on PATH)39 Pan-genome tool: pyrodigal (pip), cd-hit (external tool on PATH)40 41Run with: streamlit run app.py42"""43 44import io45import os46import shutil47import tempfile48import zipfile49from pathlib import Path50from typing import Dict, List, Optional51 52import numpy as np53import pandas as pd54import matplotlib.pyplot as plt55import streamlit as st56 57import genomics_core as gc58import genomics_figures as gf59import genomics_parsing as gp60import genomics_external as ge61 62st.set_page_config(page_title="Gulumba Genomics Suite v1.0", page_icon="\U0001F9EC", layout="wide")63gc.set_publication_style()64 65WORK_ROOT = Path(tempfile.gettempdir()) / "unified_genomics_suite"66WORK_ROOT.mkdir(parents=True, exist_ok=True)67 68 69# ============================================================================70# Shared helpers71# ============================================================================72 73def fig_to_png_bytes(fig, dpi=600) -> bytes:74 buf = io.BytesIO()75 fig.savefig(buf, format="png", dpi=dpi, bbox_inches="tight")76 buf.seek(0)77 return buf.getvalue()78 79 80def fig_to_vector_bytes(fig, fmt="pdf") -> bytes:81 buf = io.BytesIO()82 fig.savefig(buf, format=fmt, bbox_inches="tight")83 buf.seek(0)84 return buf.getvalue()85 86 87def build_zip(figures: Dict[str, bytes], tables: Dict[str, pd.DataFrame], readme: str = "") -> bytes:88 mem = io.BytesIO()89 with zipfile.ZipFile(mem, "w", compression=zipfile.ZIP_DEFLATED) as zf:90 for name, data in figures.items():91 zf.writestr(f"figures/{name}", data)92 for name, df in tables.items():93 zf.writestr(f"tables/{name}", df.to_csv(index=False))94 if readme:95 zf.writestr("README.txt", readme)96 mem.seek(0)97 return mem.getvalue()98 99 100METHODS_NOTE = (101 "Methods note: GC/AT content is reported over unambiguous (A/C/G/T) bases only, "102 "so ambiguity codes and gaps do not bias the percentage. The effective number of "103 "codons follows Wright (1990). RSCU and the codon adaptation index follow Sharp "104 "and Li (1986, 1987). Origin/terminus prediction from cumulative GC skew follows "105 "Lobry (1996) and Grigoriev (1998), and is reported with a signal-to-noise "106 "confidence flag rather than as a certainty - always confirm with dedicated tools "107 "for a final assembly. Alignment-free genome distance uses MinHash/Mash sketching "108 "(Ondov et al. 2016); ANI is estimated from the Mash Jaccard-containment distance, "109 "not from a base-level alignment, and should be read as an estimate at high "110 "similarity (>90% ANI) where the underlying approximation is most accurate."111)112 113 114def genetic_code_selector(key_prefix: str) -> gc.CodonContext:115 options = {116 "11 Bacterial, Archaeal and Plant Plastid": 11,117 "1 Standard": 1,118 "2 Vertebrate Mitochondrial": 2,119 "4 Mold, Protozoan and Coelenterate Mitochondrial": 4,120 "5 Invertebrate Mitochondrial": 5,121 "9 Echinoderm and Flatworm Mitochondrial": 9,122 }123 label = st.selectbox("Genetic code table", list(options.keys()), key=f"{key_prefix}_gcode")124 return gc.build_codon_context(options[label])125 126 127def theme_selector(key_prefix: str) -> str:128 return st.selectbox("Figure colour theme", ["Ocean", "Ember", "Forest", "Royal"], key=f"{key_prefix}_theme")129 130 131def collect_cds_sequences(records: List[gp.GenomeRecord], tables: Dict[str, pd.DataFrame],132 min_len_bp: int = 90) -> List[str]:133 """Extract CDS sequences from parsed feature tables (any input format)."""134 seqs = []135 for rec in records:136 fdf = tables.get(rec.id)137 if fdf is None or fdf.empty:138 continue139 cds = fdf[fdf["feature_type"].str.lower().isin(["cds"])]140 for _, row in cds.iterrows():141 start, end = int(row["start"]), int(row["end"])142 seq = rec.seq[start:end]143 if int(row.get("strand", 0)) == -1:144 seq = gp.reverse_complement(seq)145 seq = gc.valid_cds(seq)146 if len(seq) >= min_len_bp:147 seqs.append(seq)148 return seqs149 150 151# ============================================================================152# Page header153# ============================================================================154 155st.markdown("## \U0001F9EC Unified Genomics Suite")156st.caption(157 "Gulumba Genome composition and codon usage \u00b7 Prokka annotation workbench \u00b7 "158 "Pan-genome analysis \u00b7 Comparative genomics. One shared, tested scientific "159 "core; publication-ready figures and tables throughout."160)161 162tab_names = [163 "1. Genome Composition & Codon Usage",164 "2. Annotation Workbench (Prokka)",165 "3. Pan-Genome Analysis",166 "4. Comparative Genomics",167]168tab1, tab2, tab3, tab4 = st.tabs(tab_names)169 170 171# ============================================================================172# TOOL 1 - Gulumba Genome composition & codon usage173# (replaces GenomeScope Pro 2.1 + genome_analytics.py)174# ============================================================================175 176with tab1:177 st.markdown("#### Single-genome structure, GC content, and codon usage bias")178 st.caption(179 "Accepts GenBank, EMBL, FASTA, GFF/GTF+FASTA, BED+FASTA, an annotation "180 "table+FASTA, or a ZIP bundle of any of these."181 )182 183 c1, c2 = st.columns(2)184 with c1:185 primary_file = st.file_uploader(186 "Genome file (or ZIP bundle)", key="t1_primary",187 type=["gb", "gbk", "genbank", "embl", "fa", "fasta", "fna", "fas",188 "gff", "gff3", "gtf", "bed", "csv", "tsv", "zip", "gz"],189 )190 with c2:191 secondary_file = st.file_uploader(192 "Matching FASTA (only needed for GFF/BED/table input)", key="t1_secondary",193 type=["fa", "fasta", "fna", "fas", "gz"],194 )195 196 ctx1 = genetic_code_selector("t1")197 theme1 = theme_selector("t1")198 window_size = st.number_input("Window size for GC/skew (bp)", min_value=100, max_value=200000, value=1000, step=100, key="t1_window")199 min_coding_bp = st.number_input("Minimum CDS length for codon analysis (bp)", min_value=30, max_value=3000, value=90, step=30, key="t1_mincds")200 201 run1 = st.button("Run genome composition & codon usage analysis", type="primary", key="t1_run")202 203 if run1:204 if primary_file is None:205 st.error("Please upload a genome file.")206 else:207 try:208 with st.spinner("Parsing input..."):209 secondary_bytes = secondary_file.getvalue() if secondary_file else None210 secondary_name = secondary_file.name if secondary_file else None211 parsed = gp.parse_any(primary_file.name, primary_file.getvalue(), secondary_name, secondary_bytes)212 213 if not parsed.records:214 st.error("No sequence records were parsed from the uploaded file.")215 else:216 main_rec = max(parsed.records, key=lambda r: len(r.seq))217 feature_df = parsed.tables.get(main_rec.id, gp.empty_feature_table())218 219 with st.spinner("Computing composition, skew, and codon metrics..."):220 window_df = gc.window_table(main_rec.seq, window=int(window_size))221 ori_ter = gc.predict_ori_ter(window_df)222 nuc_counts = gc.nucleotide_counts(main_rec.seq)223 dinuc_df = gc.dinucleotide_oe(main_rec.seq)224 225 lengths = [len(r.seq) for r in parsed.records]226 assembly = gc.assembly_statistics(lengths)227 228 cds_seqs = collect_cds_sequences(parsed.records, parsed.tables, int(min_coding_bp))229 counts = gc.aggregate_codon_counts(cds_seqs, ctx1)230 stops = gc.stop_codon_counts(cds_seqs, ctx1)231 weights, weight_source, weight_note = gc.build_cai_weights(ctx1, all_seqs=cds_seqs)232 usage = gc.codon_usage_table(counts, ctx1, stops)233 per_gene = gc.per_gene_codon_table(cds_seqs, ctx1, weights, min_codons=int(min_coding_bp) // 3)234 aa_df = gc.amino_acid_usage(cds_seqs, ctx1)235 neutrality = (gc.neutrality_regression(per_gene["GC3"].values, per_gene["GC12"].values)236 if not per_gene.empty else {})237 238 with st.spinner("Rendering figures..."):239 composition_records = pd.DataFrame([{240 "record_id": r.id, "length_bp": len(r.seq),241 "gc_percent": gc.gc_percent(r.seq),242 "entropy_bits": gc.shannon_entropy(r.seq),243 **{f"count_{b}": gc.nucleotide_counts(r.seq).get(b, 0) for b in "ATGCN"},244 } for r in parsed.records])245 figs = {246 "01_composition_overview": gf.fig_composition_overview(composition_records, theme1),247 "02_skew_and_origin": gf.fig_skew_and_ori(window_df, ori_ter, main_rec.id, theme1),248 "03_dinucleotide_oe": gf.fig_dinucleotide_oe(dinuc_df, theme1),249 "04_assembly_curve": gf.fig_assembly_curve(lengths, assembly, theme1),250 "05_rscu_heatmap": gf.fig_rscu_heatmap(usage),251 "06_codon_bar": gf.fig_codon_bar(usage, theme1),252 "07_amino_acid_usage": gf.fig_amino_acid_usage(aa_df, theme1),253 "08_enc_plot": gf.fig_enc_plot(per_gene, theme1),254 "09_neutrality_plot": gf.fig_neutrality_plot(per_gene, neutrality, theme1),255 "10_pr2_plot": gf.fig_pr2_plot(per_gene, theme1),256 }257 if not feature_df.empty:258 figs["11_feature_map"] = gf.fig_feature_map(feature_df, len(main_rec.seq), main_rec.id)259 260 st.success(f"Analysis complete: {main_rec.id} ({len(main_rec.seq):,} bp, {len(parsed.records)} record(s))")261 if parsed.notes:262 st.info(" \u00b7 ".join(parsed.notes))263 264 m1, m2, m3, m4 = st.columns(4)265 m1.metric("GC%", f"{gc.gc_percent(main_rec.seq):.2f}")266 m2.metric("CDS analysed", f"{len(cds_seqs):,}")267 m3.metric("ENC (genome)", f"{gc.enc_from_counts(counts, ctx1):.2f}" if counts else "n/a")268 m4.metric("N50 (bp)", f"{assembly.get('N50_bp', 0):,}")269 270 fig_tabs = st.tabs(["Structure", "GC/skew", "Codon usage", "ENC/neutrality/PR2", "Tables"])271 with fig_tabs[0]:272 if "11_feature_map" in figs:273 st.pyplot(figs["11_feature_map"])274 st.pyplot(figs["04_assembly_curve"])275 with fig_tabs[1]:276 st.pyplot(figs["02_skew_and_origin"])277 st.pyplot(figs["01_composition_overview"])278 st.pyplot(figs["03_dinucleotide_oe"])279 st.caption(280 f"Predicted origin confidence: {ori_ter.get('ori_ter_confidence', 'n/a')} "281 f"(signal-to-noise {ori_ter.get('ori_ter_signal_to_noise', float('nan')):.2f})"282 )283 with fig_tabs[2]:284 st.pyplot(figs["05_rscu_heatmap"])285 st.pyplot(figs["06_codon_bar"])286 st.pyplot(figs["07_amino_acid_usage"])287 st.caption(f"CAI weight source: {weight_source}. {weight_note}")288 with fig_tabs[3]:289 cc1, cc2, cc3 = st.columns(3)290 cc1.pyplot(figs["08_enc_plot"])291 cc2.pyplot(figs["09_neutrality_plot"])292 cc3.pyplot(figs["10_pr2_plot"])293 with fig_tabs[4]:294 st.markdown("**Codon usage table**")295 st.dataframe(usage, use_container_width=True, height=280)296 st.markdown("**Per-gene codon metrics**")297 st.dataframe(per_gene, use_container_width=True, height=280)298 st.markdown("**Feature table**")299 st.dataframe(feature_df, use_container_width=True, height=200)300 301 png_figs = {f"{name}.png": fig_to_png_bytes(fig) for name, fig in figs.items()}302 pdf_figs = {f"{name}.pdf": fig_to_vector_bytes(fig) for name, fig in figs.items()}303 all_figs = {**png_figs, **pdf_figs}304 tables_out = {305 "codon_usage.csv": usage, "per_gene_metrics.csv": per_gene,306 "amino_acid_usage.csv": aa_df, "feature_table.csv": feature_df,307 "window_composition.csv": window_df, "dinucleotide_oe.csv": dinuc_df,308 "per_record_composition.csv": composition_records,309 }310 bundle = build_zip(all_figs, tables_out, METHODS_NOTE)311 st.download_button("Download complete results bundle (figures + tables + methods note)",312 data=bundle, file_name=f"{main_rec.id}_genomics_suite_results.zip",313 mime="application/zip", key="t1_download")314 for fig in figs.values():315 plt.close(fig)316 except Exception as exc:317 st.error(f"Analysis failed: {exc}")318 319 320# ============================================================================321# TOOL 2 - Gulumba Prokka annotation workbench322# (replaces app.py + genome_analysis.py + annotation_core.py)323# ============================================================================324 325with tab2:326 st.markdown("#### Prokka annotation workbench")327 st.caption(328 "Runs Prokka on an uploaded bacterial genome FASTA, then computes the same "329 "composition and codon-usage metrics as Tool 1 from the Prokka output, so "330 "results are directly comparable across tools."331 )332 333 if not ge.prokka_available():334 st.warning(335 "Prokka was not found on PATH in this environment. This tab will not run "336 "until Prokka is installed (e.g. `conda install -c bioconda -c conda-forge prokka`). "337 "Every other tool in this suite works without it."338 )339 340 genome_upload = st.file_uploader("Bacterial genome FASTA", type=["fa", "fasta", "fna", "fas"], key="t2_upload")341 cc1, cc2, cc3, cc4 = st.columns(4)342 with cc1:343 prefix = st.text_input("Sample prefix", value="sample_genome", key="t2_prefix")344 with cc2:345 genus = st.text_input("Genus", value="Escherichia", key="t2_genus")346 with cc3:347 species = st.text_input("Species", value="coli", key="t2_species")348 with cc4:349 cpus = st.number_input("CPUs", min_value=1, max_value=8, value=2, key="t2_cpus")350 351 ctx2 = genetic_code_selector("t2")352 theme2 = theme_selector("t2")353 run2 = st.button("Run Prokka annotation & analysis", type="primary", key="t2_run")354 355 if run2:356 if genome_upload is None:357 st.error("Please upload a genome FASTA file.")358 else:359 safe_prefix = "".join(c if c.isalnum() or c in "-_" else "_" for c in prefix).strip("_") or "sample_genome"360 job_dir = WORK_ROOT / f"prokka_{safe_prefix}"361 shutil.rmtree(job_dir, ignore_errors=True)362 job_dir.mkdir(parents=True, exist_ok=True)363 fasta_path = job_dir / f"{safe_prefix}.fasta"364 fasta_path.write_bytes(genome_upload.getvalue())365 366 with st.spinner("Running Prokka (this can take several minutes)..."):367 result = ge.run_prokka(fasta_path, job_dir / "prokka_out", safe_prefix, genus, species, int(cpus))368 369 with st.expander("Prokka log"):370 st.text(result["log"])371 372 if not result["success"]:373 st.error("Prokka did not complete successfully. See the log above for details.")374 else:375 if result["mode"] == "noanno":376 st.warning(377 "Full product-search annotation did not complete, so Prokka ran in "378 "--noanno mode. Gene calls, coordinates, and every metric below are "379 "still valid; CDS product names are simply unavailable."380 )381 try:382 with st.spinner("Loading Prokka output and computing metrics..."):383 parsed = gp.read_prokka_outputs(result["outdir"], safe_prefix)384 main_rec = max(parsed.records, key=lambda r: len(r.seq))385 feature_df = parsed.tables.get(main_rec.id, gp.empty_feature_table())386 lengths = [len(r.seq) for r in parsed.records]387 assembly = gc.assembly_statistics(lengths)388 window_df = gc.window_table(main_rec.seq, window=1000)389 ori_ter = gc.predict_ori_ter(window_df)390 cds_seqs = collect_cds_sequences(parsed.records, parsed.tables, 90)391 counts = gc.aggregate_codon_counts(cds_seqs, ctx2)392 stops = gc.stop_codon_counts(cds_seqs, ctx2)393 weights, weight_source, weight_note = gc.build_cai_weights(ctx2, all_seqs=cds_seqs)394 usage = gc.codon_usage_table(counts, ctx2, stops)395 per_gene = gc.per_gene_codon_table(cds_seqs, ctx2, weights, min_codons=30)396 aa_df = gc.amino_acid_usage(cds_seqs, ctx2)397 398 st.success(f"Loaded {len(parsed.records)} contig(s), {len(feature_df)} annotated features on the largest contig.")399 m1, m2, m3, m4 = st.columns(4)400 m1.metric("Total length (bp)", f"{assembly.get('total_length_bp', 0):,}")401 m2.metric("GC%", f"{gc.gc_percent(main_rec.seq):.2f}")402 m3.metric("N50 (bp)", f"{assembly.get('N50_bp', 0):,}")403 m4.metric("CDS analysed", f"{len(cds_seqs):,}")404 405 figs2 = {406 "01_feature_map": (gf.fig_feature_map(feature_df, len(main_rec.seq), main_rec.id)407 if not feature_df.empty else None),408 "02_assembly_curve": gf.fig_assembly_curve(lengths, assembly, theme2),409 "03_skew_and_origin": gf.fig_skew_and_ori(window_df, ori_ter, main_rec.id, theme2),410 "04_rscu_heatmap": gf.fig_rscu_heatmap(usage),411 "05_amino_acid_usage": gf.fig_amino_acid_usage(aa_df, theme2),412 }413 figs2 = {k: v for k, v in figs2.items() if v is not None}414 for name, fig in figs2.items():415 st.pyplot(fig)416 417 st.markdown("**Feature table**")418 st.dataframe(feature_df, use_container_width=True, height=240)419 st.markdown("**Codon usage table**")420 st.dataframe(usage, use_container_width=True, height=240)421 422 png_figs = {f"{n}.png": fig_to_png_bytes(f) for n, f in figs2.items()}423 bundle2 = build_zip(png_figs, {424 "feature_table.csv": feature_df, "codon_usage.csv": usage,425 "per_gene_metrics.csv": per_gene, "amino_acid_usage.csv": aa_df,426 }, METHODS_NOTE)427 st.download_button("Download annotation workbench results bundle", data=bundle2,428 file_name=f"{safe_prefix}_prokka_workbench_results.zip",429 mime="application/zip", key="t2_download")430 for fig in figs2.values():431 plt.close(fig)432 except Exception as exc:433 st.error(f"Post-annotation analysis failed: {exc}")434 435 436# ============================================================================437# TOOL 3 - Gulumba Pan-genome analysis438# (replaces pangenome_core.py + the two Pan-Genome Gradio apps)439# ============================================================================440 441with tab3:442 st.markdown("#### Pan-genome analysis")443 st.caption(444 "Upload two or more genome FASTA files (or a ZIP of them, or one multi-FASTA "445 "file where each record is a separate genome). Genes are called with "446 "pyrodigal and clustered with CD-HIT into core, accessory, and unique sets, "447 "then the pan-genome's openness is estimated with Heaps' law and Tettelin-"448 "style accumulation curves."449 )450 451 if not ge.cd_hit_available():452 st.warning("cd-hit was not found on PATH. This tab needs CD-HIT installed to run.")453 454 pg_upload = st.file_uploader(455 "Multi-FASTA, or ZIP of genome FASTA files", type=["fa", "fasta", "fna", "fas", "zip"], key="t3_upload"456 )457 pc1, pc2, pc3 = st.columns(3)458 with pc1:459 identity = st.slider("Clustering identity threshold", 0.70, 1.00, 0.90, 0.01, key="t3_identity")460 with pc2:461 coverage = st.slider("Clustering coverage threshold", 0.50, 1.00, 0.80, 0.01, key="t3_coverage")462 with pc3:463 core_threshold = st.slider("Core gene threshold (fraction of genomes)", 0.80, 1.00, 1.00, 0.01, key="t3_core")464 465 theme3 = theme_selector("t3")466 accumulation_perms = st.number_input("Permutations for accumulation-curve resampling", min_value=10, max_value=1000, value=100, step=10, key="t3_perms")467 run3 = st.button("Run pan-genome analysis", type="primary", key="t3_run")468 469 if run3:470 if pg_upload is None:471 st.error("Please upload a multi-FASTA file or a ZIP of genome FASTA files.")472 else:473 job_dir = WORK_ROOT / "pangenome_job"474 shutil.rmtree(job_dir, ignore_errors=True)475 job_dir.mkdir(parents=True, exist_ok=True)476 genome_dir = job_dir / "genomes"477 genome_dir.mkdir(parents=True, exist_ok=True)478 479 try:480 # Split into one FASTA file per genome, from either a ZIP or a multi-FASTA481 genome_files = []482 if pg_upload.name.lower().endswith(".zip"):483 with zipfile.ZipFile(io.BytesIO(pg_upload.getvalue())) as zf:484 for info in zf.infolist():485 if info.is_dir():486 continue487 name = os.path.basename(info.filename)488 if os.path.splitext(name)[1].lower() not in {".fa", ".fasta", ".fna", ".fas"}:489 continue490 out_path = genome_dir / name491 out_path.write_bytes(zf.read(info.filename))492 genome_files.append(out_path)493 else:494 records = gp.parse_fasta(pg_upload.getvalue())495 for rec in records:496 safe = "".join(c if c.isalnum() or c in "._-" else "_" for c in rec.id)497 out_path = genome_dir / f"{safe}.fasta"498 out_path.write_text(f">{rec.id}\n{rec.seq}\n")499 genome_files.append(out_path)500 501 if len(genome_files) < 2:502 st.error("At least two genomes are required for pan-genome analysis.")503 else:504 genome_names = [f.stem for f in genome_files]505 proteins_dir = job_dir / "predicted_proteins"506 proteins_dir.mkdir(exist_ok=True)507 508 with st.spinner(f"Calling genes in {len(genome_files)} genomes with pyrodigal..."):509 protein_fastas = []510 for genome_file, gname in zip(genome_files, genome_names):511 out_faa = proteins_dir / f"{gname}.faa"512 ge.predict_proteins(str(genome_file), str(out_faa), gname)513 protein_fastas.append(out_faa)514 515 pooled = job_dir / "all_proteins.faa"516 with open(pooled, "w") as out:517 for pf in protein_fastas:518 out.write(pf.read_text())519 520 with st.spinner("Clustering proteins with CD-HIT..."):521 clstr_prefix = job_dir / "clusters"522 clstr_file = ge.run_cd_hit(pooled, clstr_prefix, identity, coverage)523 clusters = ge.parse_cd_hit_clusters(clstr_file)524 525 presence, cluster_ids = ge.presence_absence_matrix(clusters, genome_names)526 class_df = ge.classify_clusters(presence, cluster_ids, len(genome_names), core_threshold)527 counts_cat = class_df["category"].value_counts().to_dict()528 529 with st.spinner("Estimating pan-genome openness (Heaps' law + accumulation curves)..."):530 accumulation = gc.pangenome_accumulation(presence, permutations=int(accumulation_perms))531 heaps = gc.heaps_law_fit(accumulation)532 533 st.success(534 f"{len(genome_names)} genomes, {len(cluster_ids)} gene clusters: "535 f"{counts_cat.get('core', 0)} core, {counts_cat.get('accessory', 0)} accessory, "536 f"{counts_cat.get('unique', 0)} unique."537 )538 st.info(f"Pan-genome openness: {heaps.get('interpretation', 'n/a')} (Heaps' gamma = {heaps.get('gamma', float('nan')):.3f})")539 540 fig_pie = gf.fig_pangenome_pie(counts_cat, theme3)541 fig_curves = gf.fig_pangenome_curves(accumulation, heaps, theme3)542 fig_pa = gf.fig_presence_absence(presence, genome_names)543 544 pt1, pt2, pt3 = st.columns(3)545 pt1.pyplot(fig_pie)546 pt2.pyplot(fig_curves)547 st.pyplot(fig_pa)548 549 st.markdown("**Cluster classification**")550 st.dataframe(class_df, use_container_width=True, height=280)551 552 presence_df = pd.DataFrame(presence, columns=genome_names)553 presence_df.insert(0, "cluster_id", cluster_ids)554 555 figs3 = {"pangenome_pie.png": fig_to_png_bytes(fig_pie),556 "pangenome_accumulation_curves.png": fig_to_png_bytes(fig_curves),557 "presence_absence_matrix.png": fig_to_png_bytes(fig_pa)}558 tables3 = {"cluster_classification.csv": class_df,559 "presence_absence_matrix.csv": presence_df,560 "accumulation_curve.csv": accumulation,561 "heaps_law_fit.csv": pd.DataFrame([heaps])}562 bundle3 = build_zip(figs3, tables3, METHODS_NOTE)563 st.download_button("Download pan-genome results bundle", data=bundle3,564 file_name="pangenome_results.zip", mime="application/zip", key="t3_download")565 for fig in (fig_pie, fig_curves, fig_pa):566 plt.close(fig)567 except Exception as exc:568 st.error(f"Pan-genome analysis failed: {exc}")569 570 571# ============================================================================572# TOOL 4 - Gulumba Comparative genomics573# (replaces comparative_core.py + its Streamlit app; adds an alignment-free574# mode for anything too large for pairwise Needleman-Wunsch alignment)575# ============================================================================576 577with tab4:578 st.markdown("#### Comparative genomics")579 st.caption(580 "Two modes. 'Alignment-free' (recommended default) uses MinHash/Mash "581 "sketching and scales to whole bacterial genomes and beyond. 'Detailed "582 "pairwise alignment' reuses the original Needleman-Wunsch approach, which "583 "gives base-level identity and conserved-region detail but only scales to "584 "genes, plasmids, amplicons, and other short sequences (a handful of Mb at "585 "most before memory and runtime become impractical)."586 )587 588 comp_files = st.file_uploader(589 "Two or more genome FASTA files", type=["fa", "fasta", "fna", "fas"],590 accept_multiple_files=True, key="t4_upload",591 )592 mode = st.radio(593 "Comparison mode",594 ["Alignment-free (Mash distance + ANI + tree)", "Detailed pairwise alignment (small sequences only)"],595 index=0, key="t4_mode",596 )597 598 if mode.startswith("Alignment-free"):599 k = st.slider("k-mer size", 11, 31, 21, 2, key="t4_k")600 sketch_size = st.slider("Sketch size (larger = more accurate, slower)", 200, 5000, 2000, 100, key="t4_sketch")601 else:602 min_conserved = st.number_input("Minimum conserved region length (bp)", min_value=5, max_value=10000, value=20, key="t4_minconserved")603 st.warning(604 "Pairwise alignment cost grows with the product of every pair of sequence "605 "lengths. Two 5 Mb bacterial chromosomes means roughly 25 x 10^12 alignment "606 "cells, which will exhaust memory. Use this mode for genes, plasmids, "607 "amplicons, or other short sequences; use alignment-free mode for whole "608 "genomes."609 )610 611 run4 = st.button("Run comparative analysis", type="primary", key="t4_run")612 613 if run4:614 if not comp_files or len(comp_files) < 2:615 st.error("Please upload at least 2 genome FASTA files.")616 else:617 try:618 genomes = []619 for f in comp_files:620 recs = gp.parse_fasta(f.getvalue())621 if not recs:622 st.warning(f"No sequences found in {f.name}; skipped.")623 continue624 combined = "".join(r.seq for r in recs)625 genomes.append((Path(f.name).stem, combined))626 627 if len(genomes) < 2:628 st.error("At least two valid genomes are required.")629 else:630 if mode.startswith("Alignment-free"):631 with st.spinner(f"Sketching {len(genomes)} genomes (k={k})..."):632 sketches = {name: gc.kmer_sketch(seq, k=int(k), sketch_size=int(sketch_size))[0]633 for name, seq in genomes}634 with st.spinner("Computing pairwise Mash distance and ANI..."):635 dist_df, ani_df = gc.distance_matrix_from_sketches(sketches, k=int(k))636 with st.spinner("Building neighbour-joining tree..."):637 newick = gc.neighbor_joining_newick(dist_df)638 639 st.success(f"Compared {len(genomes)} genomes using alignment-free Mash sketching.")640 fig_ani = gf.fig_distance_heatmap(ani_df, "Estimated ANI (%)", "ANI (%)", cmap="viridis")641 st.pyplot(fig_ani)642 st.markdown("**Estimated ANI matrix**")643 st.dataframe(ani_df.round(3), use_container_width=True)644 st.markdown("**Mash distance matrix**")645 st.dataframe(dist_df.round(5), use_container_width=True)646 st.markdown("**Neighbour-joining tree (Newick)**")647 st.code(newick, language="text")648 st.caption(649 "ANI here is estimated from Mash distance (Ondov et al. 2016), not from "650 "a base-level alignment. Treat values below ~90% ANI as a coarse ranking "651 "rather than a precise identity - the approximation degrades as sequences "652 "diverge further."653 )654 655 figs4 = {"ani_heatmap.png": fig_to_png_bytes(fig_ani)}656 tables4 = {"ani_matrix.csv": ani_df, "mash_distance_matrix.csv": dist_df}657 readme4 = METHODS_NOTE + f"\n\nNeighbour-joining tree (Newick format):\n{newick}\n"658 bundle4 = build_zip(figs4, tables4, readme4)659 st.download_button("Download comparative genomics results bundle", data=bundle4,660 file_name="comparative_genomics_alignment_free_results.zip",661 mime="application/zip", key="t4_download_af")662 plt.close(fig_ani)663 else:664 # Detailed pairwise mode: only sound for short sequences. Guard against665 # accidentally running it on whole genomes, which would hang or exhaust666 # memory rather than failing with a clear message.667 total_cells = 0668 n = len(genomes)669 for i in range(n):670 for j in range(i + 1, n):671 total_cells += len(genomes[i][1]) * len(genomes[j][1])672 if total_cells > 200_000_000:673 st.error(674 f"This comparison would require on the order of {total_cells:,} "675 "alignment cells, which is too large for pairwise alignment in this "676 "environment. Switch to alignment-free mode for genomes this size."677 )678 else:679 from Bio import pairwise2 # only needed for this specific mode680 681 def compute_identity(a, b):682 comparable = matches = 0683 for x, y in zip(a, b):684 if x == "-" and y == "-":685 continue686 comparable += 1687 if x == y:688 matches += 1689 return round(100.0 * matches / comparable, 2) if comparable else 0.0690 691 with st.spinner("Running pairwise alignments..."):692 names = [g[0] for g in genomes]693 mat = np.zeros((n, n))694 for i in range(n):695 mat[i, i] = 100.0696 for i in range(n):697 for j in range(i + 1, n):698 aln = pairwise2.align.globalms(699 genomes[i][1], genomes[j][1], 2, -1, -5, -0.5, one_alignment_only=True700 )[0]701 ident = compute_identity(aln.seqA, aln.seqB)702 mat[i, j] = mat[j, i] = ident703 identity_df = pd.DataFrame(mat, index=names, columns=names)704 705 st.success(f"Pairwise alignment complete for {n} sequences.")706 fig_id = gf.fig_distance_heatmap(identity_df, "Pairwise % identity", "% identity", cmap="viridis")707 st.pyplot(fig_id)708 st.dataframe(identity_df.round(2), use_container_width=True)709 710 figs4b = {"pairwise_identity_heatmap.png": fig_to_png_bytes(fig_id)}711 tables4b = {"pairwise_identity_matrix.csv": identity_df}712 bundle4b = build_zip(figs4b, tables4b, METHODS_NOTE)713 st.download_button("Download comparative genomics results bundle", data=bundle4b,714 file_name="comparative_genomics_pairwise_results.zip",715 mime="application/zip", key="t4_download_pw")716 plt.close(fig_id)717 except Exception as exc:718 st.error(f"Comparative analysis failed: {exc}")719 720 721st.markdown("---")722st.caption(723 "Gulumba Genomics Suite v1.0. All statistics are computed once in genomics_core.py "724 "and shared across every tool above, so a GC% or ENC value means the same thing "725 "everywhere it appears in this app."726)727