"
)
for species in species_names:
s = scores[species]
html_parts.append(
f'
{species}
{s["common_name"]}
'
f'
{s["length"]}bp
{s["gc_content"]:.1%}
'
f'
{s["protein_length"]}aa
'
f'
{s["log_likelihood"]:.2f}
'
)
html_parts.append("
")
html_parts.append('
')
html_parts.append(_make_bar_chart(
species_names,
{"Log-Likelihood": [scores[s]["log_likelihood"] for s in species_names]},
title="Carbon Log-Likelihood per Species",
value_format=".1f",
))
html_parts.append("
")
await flyte.report.replace.aio(_wrap_report("\n".join(html_parts)), do_flush=True)
result = {
"gene_name": data["gene_name"],
"description": data["description"],
"species": species_names,
"scores": scores,
}
return json.dumps(result)
# ------------------------------------------------------------------
# Task 3: Align sequences and compute similarity
# ------------------------------------------------------------------
@cpu_env.task(report=True)
async def align_and_compare(
scores_json: str,
genes_dir: flyte.io.Dir,
) -> str:
"""Align sequences with Needleman-Wunsch and compute pairwise identity.
Translates DNA to protein, builds DNA and protein identity matrices,
and generates phylogenetic trees from sequence divergence.
"""
scores_data = json.loads(scores_json)
species_names = scores_data["species"]
scores = scores_data["scores"]
gene_name = scores_data["gene_name"]
n = len(species_names)
genes_path = await genes_dir.download()
with open(os.path.join(genes_path, "genes.json")) as f:
data = json.load(f)
await flyte.report.replace.aio(_wrap_report(
f"
{gene_name} - Sequence Alignment
"
f"
Aligning {n} species with Needleman-Wunsch...
"
), do_flush=True)
# Pairwise DNA identity matrix
identity_matrix = []
for sp1 in species_names:
row = []
for sp2 in species_names:
dna1 = data["sequences"][sp1]["dna"]
dna2 = data["sequences"][sp2]["dna"]
identity = _sequence_identity(dna1, dna2)
row.append(round(identity, 4))
identity_matrix.append(row)
# Pairwise protein identity matrix
protein_matrix = []
for sp1 in species_names:
row = []
for sp2 in species_names:
identity = _sequence_identity(scores[sp1]["protein"], scores[sp2]["protein"])
row.append(round(identity, 4))
protein_matrix.append(row)
# Average pairwise identities (exclude diagonal)
dna_pairs = [identity_matrix[i][j] for i in range(n) for j in range(i + 1, n)]
prot_pairs = [protein_matrix[i][j] for i in range(n) for j in range(i + 1, n)]
avg_dna = sum(dna_pairs) / len(dna_pairs) if dna_pairs else 0
avg_prot = sum(prot_pairs) / len(prot_pairs) if prot_pairs else 0
# Most/least similar pair
best_pair = max(range(len(dna_pairs)), key=lambda k: dna_pairs[k])
worst_pair = min(range(len(dna_pairs)), key=lambda k: dna_pairs[k])
pair_indices = [(i, j) for i in range(n) for j in range(i + 1, n)]
best_sp = f"{species_names[pair_indices[best_pair][0]]}-{species_names[pair_indices[best_pair][1]]}"
worst_sp = f"{species_names[pair_indices[worst_pair][0]]}-{species_names[pair_indices[worst_pair][1]]}"
# Report
html_parts = [
f"
{gene_name} - Sequence Alignment
",
f'
Pairwise alignment using Needleman-Wunsch (match=2, mismatch=-1, gap=-2). '
f"Identity is computed as matches / aligned length from the optimal global alignment.
")
# DNA vs Protein conservation comparison
html_parts.append('
')
html_parts.append(_make_bar_chart(
species_names,
{
"DNA vs Human": [identity_matrix[0][j] for j in range(n)],
"Protein vs Human": [protein_matrix[0][j] for j in range(n)],
},
title=f"Conservation vs {species_names[0]} (DNA and Protein)",
value_format=".0%",
))
html_parts.append("
")
await flyte.report.replace.aio(_wrap_report("\n".join(html_parts)), do_flush=True)
result = {
"gene_name": gene_name,
"description": scores_data["description"],
"species": species_names,
"scores": scores,
"dna_identity_matrix": identity_matrix,
"protein_identity_matrix": protein_matrix,
}
return json.dumps(result)
# ------------------------------------------------------------------
# Task 4: Fold proteins with ESMFold
# ------------------------------------------------------------------
@gpu_env.task(report=True)
async def fold_proteins(
comparison_json: str,
max_length: int = 400,
) -> str:
"""Fold each species' translated protein with ESMFold for 3D comparison.
Returns PDB strings and pLDDT confidence scores for each species.
"""
import torch
import numpy as np
from transformers import AutoTokenizer, EsmForProteinFolding
comparison = json.loads(comparison_json)
species_names = comparison["species"]
scores = comparison["scores"]
log.info("Loading ESMFold model...")
device = "cuda" if torch.cuda.is_available() else "cpu"
tokenizer = AutoTokenizer.from_pretrained("facebook/esmfold_v1")
model = EsmForProteinFolding.from_pretrained("facebook/esmfold_v1", low_cpu_mem_usage=True)
model = model.to(device)
model.eval()
structure_data = {}
n = len(species_names)
for idx, species in enumerate(species_names):
protein = scores[species]["protein"]
if len(protein) > max_length:
log.info(f"Skipping {species} ({len(protein)} aa > {max_length} max)")
continue
log.info(f"ESMFold [{idx + 1}/{n}]: {species} ({len(protein)} aa)")
await flyte.report.replace.aio(_wrap_report(
f"
"
), do_flush=True)
inputs = tokenizer(protein, return_tensors="pt", add_special_tokens=False).to(device)
with torch.no_grad():
outputs = model(**inputs)
pdb_str = _outputs_to_pdb(outputs, protein)
plddt_raw = outputs.plddt[0].cpu().numpy()
if plddt_raw.ndim == 2:
plddt_raw = plddt_raw[-1]
plddt = plddt_raw.flatten()[:len(protein)]
if plddt.max() <= 1.0:
plddt = plddt * 100
plddt_mean = float(np.mean(plddt))
structure_data[species] = {
"pdb_str": pdb_str,
"plddt_mean": round(plddt_mean, 1),
"plddt_per_residue": [round(float(v), 1) for v in plddt[:len(protein)]],
"protein_length": len(protein),
}
log.info(f" → mean pLDDT: {plddt_mean:.1f}")
# Report with 3D viewers
n_folded = len(structure_data)
avg_plddt = sum(d["plddt_mean"] for d in structure_data.values()) / n_folded if n_folded else 0
threeDmol_script = ''
stats_html = f"""
ESMFold - Cross-Species Structure Comparison
ESMFold predicts 3D structure directly from amino acid sequence.
Comparing structures across species reveals which parts of the protein are
structurally conserved (functional core) vs divergent (surface loops, species-specific adaptations).
',
]
# Key metrics
# Average pairwise identity (exclude diagonal)
n = len(species)
dna_pairs = [dna_matrix[i][j] for i in range(n) for j in range(i + 1, n)]
protein_pairs = [protein_matrix[i][j] for i in range(n) for j in range(i + 1, n)]
avg_dna_id = sum(dna_pairs) / len(dna_pairs) if dna_pairs else 0
avg_protein_id = sum(protein_pairs) / len(protein_pairs) if protein_pairs else 0
avg_plddt = sum(d["plddt_mean"] for d in structures.values()) / len(structures) if structures else 0
html_parts.append('
')
html_parts.append(f'
{n}
Species
')
html_parts.append(f'
{avg_dna_id:.0%}
Avg DNA Identity
')
html_parts.append(f'
{avg_protein_id:.0%}
Avg Protein Identity
')
html_parts.append(f'
{avg_plddt:.1f}
Avg pLDDT
')
html_parts.append(f'
{len(structures)}
Structures Folded
')
html_parts.append("
")
# Full comparison table
html_parts.append("
Per-Species Detail
")
html_parts.append(
"
Species
Scientific Name
DNA (bp)
"
"
Protein (aa)
GC%
Carbon LL
pLDDT
"
)
for sp in species:
s = scores[sp]
plddt = structures.get(sp, {}).get("plddt_mean", "N/A")
plddt_str = f"{plddt:.1f}" if isinstance(plddt, float) else plddt
html_parts.append(
f'
{sp}
{s["common_name"]}
'
f'
{s["length"]}
{s["protein_length"]}
'
f'
{s["gc_content"]:.1%}
{s["log_likelihood"]:.2f}
'
f'
{plddt_str}
'
)
html_parts.append("
")
# GC content comparison
html_parts.append('
')
html_parts.append(_make_bar_chart(
species,
{"GC Content": [scores[s]["gc_content"] for s in species]},
title="GC Content Across Species",
value_format=".2f",
))
html_parts.append("
")
# DNA phylogenetic tree
html_parts.append("
Phylogenetic Relationships
")
html_parts.append(
'
'
"Trees built from pairwise sequence identity using UPGMA clustering. "
"Species that diverged more recently cluster together. DNA and protein trees "
"may differ when synonymous mutations dominate."
"
All 4 pipeline stages completed. View individual task reports for DNA/protein
identity heatmaps, phylogenetic trees, interactive 3D protein structures with
pLDDT confidence, Carbon log-likelihood scores, and evolutionary analysis.
"""
await flyte.report.replace.aio(_wrap_report(final_html), do_flush=True)
log.info("Pipeline complete.")
return comparison_json, structures_json
# {{/docs-fragment pipeline}}
if __name__ == "__main__":
flyte.init_from_config()
run = flyte.run(pipeline)
print(run.url)
run.wait()
```
*Source: https://github.com/unionai/unionai-examples/blob/main/v2/tutorials/genomic_gene_comparison/genomic_gene_comparison.py*
Dependencies are declared at the top of the file using the `uv` script style:
```
# /// script
# requires-python = ">=3.12"
# dependencies = [
# "flyte>=2.4.0",
# "torch>=2.9.0",
# "transformers>=4.49.0",
# "accelerate>=0.34.0",
# "numpy",
# ]
# ///
```
## Orchestrate the pipeline
The top-level `pipeline` task chains four stages: load genes, Carbon scoring, sequence alignment, ESMFold folding, and a cross-species summary report.
```
# /// script
# requires-python = ">=3.12"
# dependencies = [
# "flyte>=2.4.0",
# "torch>=2.9.0",
# "transformers>=4.49.0",
# "accelerate>=0.34.0",
# "numpy",
# ]
# main = "pipeline"
# params = ""
# ///
import json
import logging
import math
import os
import tempfile
import flyte
import flyte.io
import flyte.report
# {{docs-fragment env}}
main_img = flyte.Image.from_uv_script(__file__, name="genomic-gene-comparison", pre=True)
gpu_env = flyte.TaskEnvironment(
name="genomic-gene-comparison-gpu",
image=main_img,
resources=flyte.Resources(cpu=4, memory="32Gi", gpu=1),
)
cpu_env = flyte.TaskEnvironment(
name="genomic-gene-comparison-cpu",
image=main_img,
resources=flyte.Resources(cpu=2, memory="8Gi"),
depends_on=[gpu_env],
)
# {{/docs-fragment env}}
logging.basicConfig(level=logging.WARNING, format="%(message)s", force=True)
log = logging.getLogger(__name__)
log.setLevel(logging.INFO)
# ------------------------------------------------------------------
# Homologous gene sets - same gene across species
# ------------------------------------------------------------------
# Full-length coding sequences from NCBI RefSeq (stop codon excluded).
GENE_SETS = {
"insulin": {
"gene_name": "Insulin",
"description": "Insulin regulates blood sugar in all vertebrates. Highly conserved across 500M+ years of evolution - even fish insulin can lower blood sugar in humans. Comparing across species reveals which regions are functionally essential (conserved) vs free to drift.",
"sequences": {
"Human": {
"dna": "ATGGCCCTGTGGATGCGCCTCCTGCCCCTGCTGGCGCTGCTGGCCCTCTGGGGACCTGACCCAGCCGCAGCCTTTGTGAACCAACACCTGTGCGGCTCACACCTGGTGGAAGCTCTCTACCTAGTGTGCGGGGAACGAGGCTTCTTCTACACACCCAAGACCCGCCGGGAGGCAGAGGACCTGCAGGTGGGGCAGGTGGAGCTGGGCGGGGGCCCTGGTGCAGGCAGCCTGCAGCCCTTGGCCCTGGAGGGGTCCCTGCAGAAGCGTGGCATTGTGGAACAATGCTGTACCAGCATCTGCTCCCTCTACCAGCTGGAGAACTACTGCAAC",
"common_name": "Homo sapiens",
},
"Mouse": {
"dna": "ATGGCCCTGTGGATGCGCTTCCTGCCCCTGCTGGCCCTGCTCTTCCTCTGGGAGTCCCACCCCACCCAGGCTTTTGTCAAGCAGCACCTTTGTGGTTCCCACCTGGTGGAGGCTCTCTACCTGGTGTGTGGGGAGCGTGGCTTCTTCTACACACCCATGTCCCGCCGTGAAGTGGAGGACCCACAAGTGGCACAACTGGAGCTGGGTGGAGGCCCGGGAGCAGGTGACCTTCAGACCTTGGCACTGGAGGTGGCCCAGCAGAAGCGTGGCATTGTAGATCAGTGCTGCACCAGCATCTGCTCCCTCTACCAGCTGGAGAACTACTGCAAC",
"common_name": "Mus musculus",
},
"Chicken": {
"dna": "ATGGCTCTCTGGATCCGATCACTGCCTCTTCTGGCTCTCCTTGTCTTTTCTGGCCCTGGAACCAGCTATGCAGCTGCCAACCAGCACCTCTGTGGCTCCCACTTGGTGGAGGCTCTCTACCTGGTGTGTGGAGAGCGTGGCTTCTTCTACTCCCCCAAAGCCCGACGGGATGTCGAGCAGCCCCTAGTGAGCAGTCCCTTGCGTGGCGAGGCAGGAGTGCTGCCTTTCCAGCAGGAGGAATACGAGAAAGTCAAGCGAGGGATTGTTGAGCAATGCTGCCATAACACGTGTTCCCTCTACCAACTGGAGAACTACTGCAAC",
"common_name": "Gallus gallus",
},
"Zebrafish": {
"dna": "ATGGCAGTGTGGCTTCAGGCTGGTGCTCTGTTGGTCCTGTTGGTCGTGTCCAGTGTAAGCACTAACCCAGGCACACCGCAGCACCTGTGTGGATCTCATCTGGTCGATGCCCTTTATCTGGTCTGTGGCCCAACAGGCTTCTTCTACAACCCCAAGAGAGACGTTGAGCCCCTTCTGGGTTTCCTTCCTCCTAAATCTGCCCAGGAAACTGAGGTGGCTGACTTTGCATTTAAAGATCATGCCGAGCTGATAAGGAAGAGAGGCATTGTAGAGCAGTGCTGCCACAAACCCTGCAGCATCTTTGAGCTGCAGAACTACTGTAAC",
"common_name": "Danio rerio",
},
"Frog": {
"dna": "ATGGCTCTATGGATGCAGTGTCTGCCCCTGGTTCTTGTCCTCTTTTTCTCTACACCCAACACCGAAGCTCTAGTTAACCAGCACTTGTGTGGGTCTCACCTGGTAGAAGCCCTGTACTTAGTATGTGGGGATCGAGGCTTCTTCTACTACCCTAAGGTCAAACGGGACATGGAACAAGCACTTGTCAGTGGACCCCAGGATAATGAGTTGGATGGAATGCAGCTCCAGCCTCAGGAGTATCAGAAAATGAAGAGGGGGATTGTGGAGCAATGTTGCCACAGCACATGTTCTCTCTTCCAGCTGGAGAGTTACTGCAAC",
"common_name": "Xenopus laevis",
},
"Cow": {
"dna": "ATGGCCCTGTGGACACGCCTGGCGCCCCTGCTGGCCCTGCTGGCGCTCTGGGCCCCCGCCCCGGCCCGCGCCTTCGTCAACCAGCATCTGTGTGGCTCCCACCTGGTGGAGGCGCTGTACCTGGTGTGCGGAGAGCGCGGCTTCTTCTACACGCCCAAGGCCCGCCGGGAGGTGGAGGGCCCCCAGGTGGGGGCGCTGGAGCTGGCCGGAGGCCCGGGCGCGGGCGGCCTGGAGGGGCCCCCGCAGAAGCGTGGCATCGTGGAGCAGTGCTGTGCCAGCGTCTGCTCGCTCTACCAGCTGGAGAACTACTGTAAC",
"common_name": "Bos taurus",
},
},
},
"hemoglobin": {
"gene_name": "Hemoglobin Beta",
"description": "Beta-globin carries oxygen from lungs to tissues. The most studied gene in molecular evolution - sequence differences power the 'molecular clock' hypothesis. Sickle cell mutation (E6V) in humans shows how a single base change creates devastating disease.",
"sequences": {
"Human": {
"dna": "ATGGTGCATCTGACTCCTGAGGAGAAGTCTGCCGTTACTGCCCTGTGGGGCAAGGTGAACGTGGATGAAGTTGGTGGTGAGGCCCTGGGCAGGCTGCTGGTGGTCTACCCTTGGACCCAGAGGTTCTTTGAGTCCTTTGGGGATCTGTCCACTCCTGATGCTGTTATGGGCAACCCTAAGGTGAAGGCTCATGGCAAGAAAGTGCTCGGTGCCTTTAGTGATGGCCTGGCTCACCTGGACAACCTCAAGGGCACCTTTGCCACACTGAGTGAGCTGCACTGTGACAAGCTGCACGTGGATCCTGAGAACTTCAGGCTCCTGGGCAACGTGCTGGTCTGTGTGCTGGCCCATCACTTTGGCAAAGAATTCACCCCACCAGTGCAGGCTGCCTATCAGAAAGTGGTGGCTGGTGTGGCTAATGCCCTGGCCCACAAGTATCAC",
"common_name": "Homo sapiens",
},
"Mouse": {
"dna": "ATGGTGCACCTGACTGATGCTGAGAAGGCTGCTGTCTCTGGCCTGTGGGGAAAGGTGAACGCCGATGAAGTTGGTGGTGAGGCCCTGGGCAGGCTGCTGGTTGTCTACCCTTGGACCCAGCGGTACTTTGATAGCTTTGGAGACCTATCCTCTGCCTCTGCTATCATGGGTAATGCCAAAGTGAAGGCCCATGGCAAGAAAGTGATAACTGCCTTTAACGATGGCCTGAATCACTTGGACAGCCTCAAGGGCACCTTTGCCAGCCTCAGTGAGCTCCACTGTGACAAGCTGCATGTGGATCCTGAGAACTTCAGGCTCCTGGGCAATATGATCGTGATTGTGCTGGGCCACCACCTGGGCAAGGATTTCACCCCCGCTGCACAGGCTGCCTTCCAGAAGGTGGTGGCTGGAGTGGCTGCTGCCCTGGCTCACAAGTACCAC",
"common_name": "Mus musculus",
},
"Chicken": {
"dna": "ATGGTGCACTGGACTGCTGAGGAGAAGCAGCTCATCACCGGCCTCTGGGGCAAGGTCAATGTGGCCGAATGTGGGGCTGAAGCCCTGGCCAGGCTGCTGATCGTCTACCCCTGGACCCAGAGGTTCTTTGCGTCCTTTGGGAACCTCTCCAGCCCCACTGCCATCCTTGGCAACCCCATGGTCCGCGCCCATGGCAAGAAAGTGCTCACCTCCTTTGGGGATGCTGTGAAGAACCTGGACAACATCAAGAACACCTTCTCCCAACTGTCCGAACTGCATTGTGACAAGCTGCATGTGGACCCCGAGAACTTCAGGCTCCTGGGTGACATCCTCATCATTGTCCTGGCCGCCCACTTCAGCAAGGACTTCACTCCTGAATGCCAGGCTGCCTGGCAGAAGCTGGTCCGCGTGGTGGCCCATGCCCTGGCTCGCAAGTACCAC",
"common_name": "Gallus gallus",
},
"Zebrafish": {
"dna": "ATGGTTGAGTGGACAGATGCCGAGCGCACAGCCATCCTTGGCCTGTGGGGAAAGCTCAATATCGATGAAATCGGACCTCAGGCCCTATCCAGATGTCTGATCGTGTATCCCTGGACTCAGAGATATTTCGCCACATTCGGCAACCTGTCAAGCCCCGCTGCGATCATGGGTAACCCCAAAGTGGCAGCTCATGGGAGGACTGTGATGGGAGGTCTTGAGAGAGCCATCAAGAACATGGACAACGTCAAGAACACCTATGCCGCCCTCAGTGTGATGCACTCTGAGAAACTGCATGTGGATCCCGACAACTTCAGGCTTCTCGCTGATTGCATCACCGTTTGCGCTGCCATGAAGTTCGGCCAAGCTGGTTTCAATGCTGATGTCCAGGAGGCCTGGCAGAAGTTTCTGGCTGTGGTCGTTTCTGCTCTGTGCAGACAGTACCAC",
"common_name": "Danio rerio",
},
"Frog": {
"dna": "ATGGTTCATTGGACAGCTGAAGAGAAGGCCGCCATCACCTCTGTGTGGCAGGAGGTCAACCAGGAGCAAGATGGCCATGATGCACTCACAAGGCTGCTGGTTGTGTACCCCTGGACCCAGAGATACTTCAGCAGTTTTGGAAATCTCGGTAATGCCACAGCTATTGCTGGAAATGTCAAGGTGCGTGCCCATGGCAAGAAGGTTCTTTCAGCTGTTGGTGATGCCATCGCCCATCTTGACAACGTGAAGGGAACTCTCCATGACCTCAGTGTGGTCCACGCCTTCAAGCTCTATGTGGATCCTGAGAACTTCAAGCGTCTTGGTGAAGTTCTGGTCATTGTCTTGGCTTCCAAACTGGGATCAGCCTTTACTCCTCAAGTCCAGGGAGCCTGGGAGAAATTTGTTGCTGTTCTGGTTGATGCCCTCAGCCAAGGATACAAC",
"common_name": "Xenopus laevis",
},
"Cow": {
"dna": "ATGCTGACTGCTGAGGAGAAGGCTGCCGTCACCGCCTTTTGGGGCAAGGTGAAAGTGGATGAAGTTGGTGGTGAGGCCCTGGGCAGGCTGCTGGTTGTCTACCCCTGGACTCAGAGGTTCTTTGAGTCCTTTGGGGACTTGTCCACTGCTGATGCTGTTATGAACAACCCTAAGGTGAAGGCCCATGGCAAGAAGGTGCTAGATTCCTTTAGTAATGGCATGAAGCATCTCGATGACCTCAAGGGCACCTTTGCTGCGCTGAGTGAGCTGCACTGTGATAAGCTGCATGTGGATCCTGAGAACTTCAAGCTCCTGGGCAACGTGCTAGTGGTTGTGCTGGCTCGCAATTTTGGCAAGGAATTCACCCCGGTGCTGCAGGCTGACTTTCAGAAGGTGGTGGCTGGTGTGGCCAATGCCCTGGCCCACAGATATCAT",
"common_name": "Bos taurus",
},
},
},
"p53": {
"gene_name": "p53 (TP53)",
"description": "The 'guardian of the genome' - p53 detects DNA damage and triggers repair or cell death. Mutated in >50% of human cancers. Elephants have 20 copies of p53 (humans have 1), which may explain their extremely low cancer rates despite their size (Peto's paradox).",
"sequences": {
"Human": {
"dna": "ATGGAGGAGCCGCAGTCAGATCCTAGCGTCGAGCCCCCTCTGAGTCAGGAAACATTTTCAGACCTATGGAAACTACTTCCTGAAAACAACGTTCTGTCCCCCTTGCCGTCCCAAGCAATGGATGATTTGATGCTGTCCCCGGACGATATTGAACAATGGTTCACTGAAGACCCAGGTCCAGATGAAGCTCCCAGAATGCCAGAGGCTGCTCCCCCCGTGGCCCCTGCACCAGCAGCTCCTACACCGGCGGCCCCTGCACCAGCCCCCTCCTGGCCCCTGTCATCTTCTGTCCCTTCCCAGAAAACCTACCAGGGCAGCTACGGTTTCCGTCTGGGCTTCTTGCATTCTGGGACAGCCAAGTCTGTGACTTGCACGTACTCCCCTGCCCTCAACAAGATGTTTTGCCAACTGGCCAAGACCTGCCCTGTGCAGCTGTGGGTTGATTCCACACCCCCGCCCGGCACCCGCGTCCGCGCCATGGCCATCTACAAGCAGTCACAGCACATGACGGAGGTTGTGAGGCGCTGCCCCCACCATGAGCGCTGCTCAGATAGCGATGGTCTGGCCCCTCCTCAGCATCTTATCCGAGTGGAAGGAAATTTGCGTGTGGAGTATTTGGATGACAGAAACACTTTTCGACATAGTGTGGTGGTGCCCTATGAGCCGCCTGAGGTTGGCTCTGACTGTACCACCATCCACTACAACTACATGTGTAACAGTTCCTGCATGGGCGGCATGAACCGGAGGCCCATCCTCACCATCATCACACTGGAAGACTCCAGTGGTAATCTACTGGGACGGAACAGCTTTGAGGTGCGTGTTTGTGCCTGTCCTGGGAGAGACCGGCGCACAGAGGAAGAGAATCTCCGCAAGAAAGGGGAGCCTCACCACGAGCTGCCCCCAGGGAGCACTAAGCGAGCACTGCCCAACAACACCAGCTCCTCTCCCCAGCCAAAGAAGAAACCACTGGATGGAGAATATTTCACCCTTCAGATCCGTGGGCGTGAGCGCTTCGAGATGTTCCGAGAGCTGAATGAGGCCTTGGAACTCAAGGATGCCCAGGCTGGGAAGGAGCCAGGGGGGAGCAGGGCTCACTCCAGCCACCTGAAGTCCAAAAAGGGTCAGTCTACCTCCCGCCATAAAAAACTCATGTTCAAGACAGAAGGGCCTGACTCAGAC",
"common_name": "Homo sapiens",
},
"Mouse": {
"dna": "ATGACTGCCATGGAGGAGTCACAGTCGGATATCAGCCTCGAGCTCCCTCTGAGCCAGGAGACATTTTCAGGCTTATGGAAACTACTTCCTCCAGAAGATATCCTGCCATCACCTCACTGCATGGACGATCTGTTGCTGCCCCAGGATGTTGAGGAGTTTTTTGAAGGCCCAAGTGAAGCCCTCCGAGTGTCAGGAGCTCCTGCAGCACAGGACCCTGTCACCGAGACCCCTGGGCCAGTGGCCCCTGCCCCAGCCACTCCATGGCCCCTGTCATCTTTTGTCCCTTCTCAAAAAACTTACCAGGGCAACTATGGCTTCCACCTGGGCTTCCTGCAGTCTGGGACAGCCAAGTCTGTTATGTGCACGTACTCTCCTCCCCTCAATAAGCTATTCTGCCAGCTGGCGAAGACGTGCCCTGTGCAGTTGTGGGTCAGCGCCACACCTCCAGCTGGGAGCCGTGTCCGCGCCATGGCCATCTACAAGAAGTCACAGCACATGACGGAGGTCGTGAGACGCTGCCCCCACCATGAGCGCTGCTCCGATGGTGATGGCCTGGCTCCTCCCCAGCATCTTATCCGGGTGGAAGGAAATTTGTATCCCGAGTATCTGGAAGACAGGCAGACTTTTCGCCACAGCGTGGTGGTACCTTATGAGCCACCCGAGGCCGGCTCTGAGTATACCACCATCCACTACAAGTACATGTGTAATAGCTCCTGCATGGGGGGCATGAACCGCCGACCTATCCTTACCATCATCACACTGGAAGACTCCAGTGGGAACCTTCTGGGACGGGACAGCTTTGAGGTTCGTGTTTGTGCCTGCCCTGGGAGAGACCGCCGTACAGAAGAAGAAAATTTCCGCAAAAAGGAAGTCCTTTGCCCTGAACTGCCCCCAGGGAGCGCAAAGAGAGCGCTGCCCACCTGCACAAGCGCCTCTCCCCCGCAAAAGAAAAAACCACTTGATGGAGAGTATTTCACCCTCAAGATCCGCGGGCGTAAACGCTTCGAGATGTTCCGGGAGCTGAATGAGGCCTTAGAGTTAAAGGATGCCCATGCTACAGAGGAGTCTGGAGACAGCAGGGCTCACTCCAGCTACCTGAAGACCAAGAAGGGCCAGTCTACTTCCCGCCATAAAAAAACAATGGTCAAGAAAGTGGGGCCTGACTCAGAC",
"common_name": "Mus musculus",
},
"Chicken": {
"dna": "ATGGCGGAGGAGATGGAACCATTGCTGGAACCCACTGAGGTCTTCATGGACCTCTGGAGCATGCTCCCCTATAGCATGCAACAGCTGCCCCTCCCTGAGGATCACAGCAACTGGCAGGAGCTGAGCCCCCTGGAACCCAGCGACCCCCCCCCACCACCGCCACCACCACCTCTGCCATTGGCCGCCGCCGCCCCCCCCCCATTAAACCCCCCCACCCCCCCCCGCGCTGCCCCCTCCCCGGTGGTCCCATCCACGGAGGATTATGGGGGGGACTTCGACTTCCGGGTGGGGTTCGTGGAGGCGGGCACAGCCAAATCGGTCACCTGCACTTACTCCCCGGTGCTGAATAAGGTCTATTGCCGCCTGGCCAAGCCGTGCCCGGTGCAGGTGAGGGTGGGGGTGGCGCCCCCCCCCGGTTCCTCCCTCCGCGCCGTGGCCGTCTATAAGAAATCAGAGCACGTGGCCGAAGTGGTGCGGCGCTGCCCCCACCACGAGCGCTGCGGGGGGGGCACCGACGGCCTGGCCCCCGCACAGCACCTCATCCGGGTGGAGGGGAACCCCCAGGCGCGTTACCACGACGACGAGACCACCAAACGGCACAGCGTCGTCGTCCCCTATGAGCCCCCCGAGGTGGGCTCTGACTGTACCACGGTGCTGTACAACTTCATGTGCAACAGTTCCTGCATGGGGGGGATGAACCGCCGCCCCATCCTCACCATCCTTACACTGGAGGGGCCGGGGGGGCAGCTGTTGGGGCGGCGCTGCTTCGAGGTGCGCGTGTGCGCATGTCCGGGGAGGGACCGCAAGATCGAGGAGGAGAACTTCCGCAAGAGGGGCGGGGCCGGGGGCGTGGCTAAGCGAGCCATGTCGCCCCCAACCGAAGCCCCCGAGCCCCCCAAGAAGCGCGTGCTGAACCCCGACAATGAGATATTCTACCTGCAGGTGCGCGGGCGCCGCCGCTATGAGATGCTGAAGGAGATCAATGAGGCGCTGCAGCTCGCCGAGGGGGGGTCCGCACCGCGGCCTTCCAAAGGCCGCCGTGTGAAGGTGGAGGGACCCCAACCCAGCTGCGGGAAGAAACTGCTGCAAAAAGGCTCGGAC",
"common_name": "Gallus gallus",
},
"Zebrafish": {
"dna": "ATGGCGCAAAACGACAGCCAAGAGTTCGCGGAGCTCTGGGAGAAGAATTTGATTATTCAGCCCCCAGGTGGTGGCTCTTGCTGGGACATCATTAATGATGAGGAGTACTTGCCGGGATCGTTTGACCCCAATTTTTTTGAAAATGTGCTTGAAGAACAGCCTCAGCCATCCACTCTCCCACCAACATCCACTGTTCCGGAGACAAGCGACTATCCCGGCGATCATGGATTTAGGCTCAGGTTCCCGCAGTCTGGCACAGCAAAATCTGTAACTTGCACTTATTCACCGGACCTGAATAAACTCTTCTGTCAGCTGGCAAAAACTTGCCCCGTTCAAATGGTGGTGGACGTTGCCCCTCCACAGGGCTCCGTGGTTCGAGCCACTGCCATCTATAAGAAGTCCGAGCATGTGGCTGAAGTGGTCCGCAGATGCCCCCATCATGAGCGAACCCCGGATGGAGATAACTTGGCGCCTGCTGGTCATTTGATAAGAGTGGAGGGCAATCAGCGAGCAAATTACAGGGAAGATAACATCACTTTAAGGCATAGTGTTTTTGTCCCATATGAAGCACCACAGCTTGGTGCTGAATGGACAACTGTGCTACTAAACTACATGTGCAATAGCAGCTGCATGGGGGGGATGAACCGCAGGCCCATCCTCACAATCATCACTCTGGAGACTCAGGAAGGTCAGTTGCTGGGCCGGAGGTCTTTTGAGGTGCGTGTGTGTGCATGTCCAGGCAGAGACAGGAAAACTGAGGAGAGCAACTTCAAGAAAGACCAAGAGACCAAAACCATGGCCAAAACCACCACTGGGACCAAACGTAGTTTGGTGAAAGAATCTTCTTCAGCTACATTACGACCTGAGGGGAGCAAAAAGGCCAAGGGCTCCAGCAGCGATGAGGAGATCTTTACCCTGCAGGTGAGGGGCAGGGAGCGTTATGAAATTTTAAAGAAATTGAACGACAGTCTGGAGTTAAGTGATGTGGTGCCTGCCTCAGATGCTGAAAAGTATCGTCAGAAATTCATGACAAAAAACAAAAAAGAGAATCGTGAATCATCTGAGCCCAAACAGGGAAAGAAGCTGATGGTGAAGGACGAAGGAAGAAGCGACTCTGAT",
"common_name": "Danio rerio",
},
"Elephant": {
"dna": "ATGGAGGAGCCCCAGTCAGATCTCAGCACCGAGCTCCCTCTGAGTCAAGAGACGTTTTCATACTTATGGGAACTCCTTCCTGAGAATCCGGTTCTGTCCCCCACACTACCCCCGGCAGTGGAGGTCATGGACGATCTGCTACTCTCAGAAGACACTGCAAACTGGCTAGAAAGCCAAGTTGAGGCTCAGGGAATGTCCACAACCCCTGCACCAGCCACCCCTACACCGGTGGCCCCCGCACCAGCCACCTCCTGGACCCTGTCATCTTCCGTCCCTTCCCAAAAGACCTACCCTGGCACCTATGGTTTCCGTCTGGGCTTCCTACATTCTGGGACAGCCAAGTCCGTCACCTGCACGTACTCCCCTGACCTTAACAAGCTGTTTTGCCAGCTGGCAAAAACCTGCCCAGTGCAGCTGTGGGTCGCCTCACCACCCCCGCCCGGCACCCGTGTTCGCACCATGGCCATCTACAAGAAGTCAGAGCATATGACGGAGGTCGTCAAGCGCTGCCCCCACCATGAGCGCTGCTCTGACTCTAGCGATGGCCTGGCCCCTCCTCAGCACCTCATCCGGGTGGAAGGAAACCTGCGTGCTGAGTATCTGGAGGACAGCATCACTCTCCGACACAGTGTGGTGGTGCCCTACGAGCCGCCCGAGGTTGGGTCTGACTGTACCACCATCCACTTCAACTTCATGTGTAACAGCTCCTGCATGGGGGGCATGAACCGGCGGCCCATCCTCACCATCATCACACTGGAAGACTCCAGTGGTAATCTGCTGGGACGTAACAGCTTTGAGGTGCGCATTTGTGCCTGTCCTGGAAGAGACAGACGTACAGAAGAAGAAAATTTCCACAAGAAGGGAGAGCCTTGCCCAGAGCCGCCACCCCCTGGGAGGAGCACTAAGCGAGCACTGCCCACCAACACCAGCTCCTCTACCCAGCCAAAGAAGAAGCCACTGGATGAAGAATATTTCACCCTTCAGATCCGTGGGCGTGAACGCTTCAAGATGTTCCTAGAGCTAAATGAGGCCTTGGAGCTGAAGGATGCCCAGGCTGGGAAGGAGCCAGAGGGGAGCCGGGCTCACTCCAGCCCTTCGAAGTCTAAGAAGGGACAGTCTACCTCCCGCCATAAAAAACCAATGTTCAAGAGAGAGGGACCTGACTCAGAC",
"common_name": "Loxodonta africana",
},
"Dog": {
"dna": "ATGGAGGAGTCGCAGTCAGAGCTCAATATCGACCCCCCTCTGAGCCAGGAGACATTTTCAGAATTGTGGAACCTGCTTCCTGAAAACAATGTTCTGTCTTCGGAGCTGTGCCCAGCAGTGGATGAGCTGCTGCTCCCAGAGAGCGTCGTGAACTGGCTAGACGAAGACTCAGATGATGCTCCCAGGATGCCAGCCACTTCTGCCCCCACAGCCCCTGGACCGGCCCCCTCGTGGCCCCTATCATCCTCTGTCCCTTCCCCGAAGACCTACCCTGGCACCTATGGGTTCCGTTTGGGGTTCCTGCATTCCGGGACAGCCAAGTCTGTTACTTGGACGTACTCCCCTCTCCTCAACAAGTTGTTTTGCCAGCTGGCGAAGACCTGCCCCGTGCAGCTGTGGGTCAGCTCCCCACCCCCACCCAATACCTGCGTCCGCGCTATGGCCATCTATAAGAAGTCGGAGTTCGTGACCGAGGTTGTGCGGCGCTGCCCCCACCATGAACGCTGCTCTGACAGTAGTGACGGTCTTGCCCCTCCTCAGCATCTCATCCGAGTGGAAGGAAATTTGCGGGCCAAGTACCTGGACGACAGAAACACTTTTCGACACAGTGTGGTGGTGCCTTATGAGCCACCCGAGGTTGGCTCTGACTATACCACCATCCACTACAACTACATGTGTAACAGTTCCTGCATGGGAGGCATGAACCGGCGGCCCATCCTCACTATCATCACCCTGGAAGACTCCAGTGGAAACGTGCTGGGACGCAACAGCTTTGAGGTACGCGTTTGTGCCTGTCCCGGGAGAGACCGCCGGACTGAGGAGGAGAATTTCCACAAGAAGGGGGAGCCTTGTCCTGAGCCACCCCCCGGGAGTACCAAGCGAGCACTGCCTCCCAGCACCAGCTCCTCTCCCCCGCAAAAGAAGAAGCCACTAGATGGAGAATATTTCACCCTTCAGATCCGTGGGCGTGAACGCTATGAGATGTTCAGGAATCTGAATGAAGCCTTGGAGCTGAAGGATGCCCAGAGTGGAAAGGAGCCAGGGGGAAGCAGGGCTCACTCCAGCCACCTGAAGGCAAAGAAGGGGCAATCTACCTCTCGCCATAAAAAACTGATGTTCAAGAGAGAAGGGCTTGACTCAGAC",
"common_name": "Canis lupus familiaris",
},
},
},
}
# Standard genetic code
CODON_TABLE = {
"TTT": "F", "TTC": "F", "TTA": "L", "TTG": "L", "CTT": "L", "CTC": "L",
"CTA": "L", "CTG": "L", "ATT": "I", "ATC": "I", "ATA": "I", "ATG": "M",
"GTT": "V", "GTC": "V", "GTA": "V", "GTG": "V", "TCT": "S", "TCC": "S",
"TCA": "S", "TCG": "S", "CCT": "P", "CCC": "P", "CCA": "P", "CCG": "P",
"ACT": "T", "ACC": "T", "ACA": "T", "ACG": "T", "GCT": "A", "GCC": "A",
"GCA": "A", "GCG": "A", "TAT": "Y", "TAC": "Y", "TAA": "*", "TAG": "*",
"CAT": "H", "CAC": "H", "CAA": "Q", "CAG": "Q", "AAT": "N", "AAC": "N",
"AAA": "K", "AAG": "K", "GAT": "D", "GAC": "D", "GAA": "E", "GAG": "E",
"TGT": "C", "TGC": "C", "TGA": "*", "TGG": "W", "CGT": "R", "CGC": "R",
"CGA": "R", "CGG": "R", "AGT": "S", "AGC": "S", "AGA": "R", "AGG": "R",
"GGT": "G", "GGC": "G", "GGA": "G", "GGG": "G",
}
BASE_COLORS = {"A": "#2ecc71", "T": "#e74c3c", "G": "#f39c12", "C": "#3498db"}
def _translate(dna: str) -> str:
"""Translate DNA to protein in reading frame 0."""
dna = dna.upper()
protein = []
for i in range(0, len(dna) - 2, 3):
codon = dna[i:i + 3]
aa = CODON_TABLE.get(codon, "X")
if aa == "*":
break
protein.append(aa)
return "".join(protein)
def _gc_content(seq: str) -> float:
if not seq:
return 0.0
return sum(1 for b in seq.upper() if b in "GC") / len(seq)
def _sequence_identity(seq1: str, seq2: str, match: int = 2, mismatch: int = -1, gap: int = -2) -> float:
"""Percent identity via Needleman-Wunsch global alignment."""
if not seq1 or not seq2:
return 0.0
n, m = len(seq1), len(seq2)
# Build score matrix
dp = [[0] * (m + 1) for _ in range(n + 1)]
for i in range(1, n + 1):
dp[i][0] = dp[i - 1][0] + gap
for j in range(1, m + 1):
dp[0][j] = dp[0][j - 1] + gap
for i in range(1, n + 1):
for j in range(1, m + 1):
s = match if seq1[i - 1] == seq2[j - 1] else mismatch
dp[i][j] = max(dp[i - 1][j - 1] + s, dp[i - 1][j] + gap, dp[i][j - 1] + gap)
# Traceback to count matches and alignment length
i, j = n, m
matches = 0
aligned = 0
while i > 0 or j > 0:
if i > 0 and j > 0:
s = match if seq1[i - 1] == seq2[j - 1] else mismatch
if dp[i][j] == dp[i - 1][j - 1] + s:
if seq1[i - 1] == seq2[j - 1]:
matches += 1
aligned += 1
i -= 1
j -= 1
continue
if i > 0 and dp[i][j] == dp[i - 1][j] + gap:
aligned += 1
i -= 1
else:
aligned += 1
j -= 1
return matches / aligned if aligned else 0.0
# ------------------------------------------------------------------
# Report styling
# ------------------------------------------------------------------
REPORT_CSS = """
"""
def _wrap_report(html: str) -> str:
return f'{REPORT_CSS}
{html}
'
# ------------------------------------------------------------------
# SVG chart helpers
# ------------------------------------------------------------------
def _make_heatmap(
matrix: list[list[float]],
row_labels: list[str],
col_labels: list[str],
title: str = "",
width: int = 600,
height: int = 500,
value_format: str = ".1f",
color_scale: str = "blue",
) -> str:
"""Generate an SVG heatmap."""
n_rows = len(matrix)
n_cols = len(matrix[0]) if matrix else 0
if not n_rows or not n_cols:
return ""
show_values = n_rows <= 10 and n_cols <= 10
flat = [v for row in matrix for v in row]
v_min = min(flat)
v_max = max(flat)
v_range = v_max - v_min or 1
if color_scale == "blue":
def get_color(v):
t = (v - v_min) / v_range
r = int(255 - t * (255 - 30))
g = int(255 - t * (255 - 58))
b = int(255 - t * (255 - 95))
return f"rgb({r},{g},{b})"
else: # green
def get_color(v):
t = (v - v_min) / v_range
r = int(255 - t * (255 - 6))
g = int(255 - t * (255 - 95))
b = int(255 - t * (255 - 70))
return f"rgb({r},{g},{b})"
ml = max(80, max(len(l) for l in row_labels) * 7 + 10) if row_labels else 80
mr = 20
mt = 80
mb = 20
cw = width - ml - mr
ch = height - mt - mb
cell_w = cw / n_cols
cell_h = ch / n_rows
svg = [
f'")
return "\n".join(svg)
def _make_dendrogram(
names: list[str],
matrix: list[list[float]],
title: str = "",
width: int = 700,
height: int = 350,
color: str = "#2563eb",
) -> str:
"""Generate an SVG dendrogram from a similarity matrix using UPGMA."""
n = len(names)
if n < 2:
return ""
dist = [[1.0 - matrix[i][j] for j in range(n)] for i in range(n)]
clusters = [{"members": [i], "height": 0.0, "left": None, "right": None} for i in range(n)]
active = list(range(n))
while len(active) > 1:
best_d = float("inf")
bi, bj = 0, 1
for ii in range(len(active)):
for jj in range(ii + 1, len(active)):
ci, cj = active[ii], active[jj]
d = 0
count = 0
for mi in clusters[ci]["members"]:
for mj in clusters[cj]["members"]:
d += dist[mi][mj]
count += 1
avg_d = d / count if count else 0
if avg_d < best_d:
best_d = avg_d
bi, bj = ii, jj
ci, cj = active[bi], active[bj]
new_cluster = {
"members": clusters[ci]["members"] + clusters[cj]["members"],
"height": best_d,
"left": clusters[ci],
"right": clusters[cj],
}
clusters.append(new_cluster)
new_idx = len(clusters) - 1
active.pop(bj)
active.pop(bi)
active.append(new_idx)
root = clusters[active[0]]
max_label_len = max((len(n) for n in names), default=0)
ml, mr, mt, mb = max(50, max_label_len * 5 + 10), 30, 40, 80
cw = width - ml - mr
ch = height - mt - mb
max_h = root["height"] or 1
leaf_positions = {}
leaf_counter = [0]
def assign_leaves(node):
if node["left"] is None and node["right"] is None:
leaf_positions[node["members"][0]] = leaf_counter[0]
leaf_counter[0] += 1
else:
if node["left"]:
assign_leaves(node["left"])
if node["right"]:
assign_leaves(node["right"])
assign_leaves(root)
n_leaves = len(leaf_positions)
leaf_spacing = cw / max(n_leaves - 1, 1)
svg = [
f'")
return "\n".join(svg)
def _make_bar_chart(
labels: list[str],
series: dict[str, list[float]],
title: str = "",
colors: list[str] | None = None,
width: int = 700,
height: int = 300,
value_format: str = ".1f",
) -> str:
"""Generate an SVG grouped bar chart."""
if not labels:
return ""
default_colors = ["#2563eb", "#059669", "#f59e0b", "#dc2626", "#7c3aed"]
colors = colors or default_colors
ml, mr, mt, mb = 60, 20, 40, 80
cw = width - ml - mr
ch = height - mt - mb
all_vals = [v for vals in series.values() for v in vals]
y_min = min(all_vals) if all_vals else 0
y_max = max(all_vals) if all_vals else 1
if y_min >= 0:
y_min_plot = 0
y_max_plot = y_max * 1.15 or 1
else:
y_range = y_max - y_min or 1
y_min_plot = y_min - y_range * 0.05
y_max_plot = y_max + y_range * 0.15
n_groups = len(labels)
n_series = len(series)
group_width = cw / n_groups
bar_width = group_width * 0.7 / max(n_series, 1)
gap = group_width * 0.15
def sy(v):
return mt + ch - ((v - y_min_plot) / (y_max_plot - y_min_plot)) * ch
svg = [
f'")
return "\n".join(svg)
def _make_plddt_sparkline(values: list[float], width: int = 400, height: int = 50) -> str:
"""pLDDT sparkline with AlphaFold-style coloring."""
if not values or len(values) < 2:
return ""
pad = 4
cw = width - 2 * pad
ch = height - 2 * pad
seg_w = cw / len(values)
svg = [
f'")
return "\n".join(svg)
def _outputs_to_pdb(outputs, sequence: str) -> str:
"""Convert ESMFold outputs to PDB format string."""
import numpy as np
pos = outputs.positions[0]
if pos.dim() == 4:
pos = pos[-1]
positions = pos.cpu().numpy()
atom_names = ["N", "CA", "C", "O"]
aa_3letter = {
"A": "ALA", "R": "ARG", "N": "ASN", "D": "ASP", "C": "CYS",
"Q": "GLN", "E": "GLU", "G": "GLY", "H": "HIS", "I": "ILE",
"L": "LEU", "K": "LYS", "M": "MET", "F": "PHE", "P": "PRO",
"S": "SER", "T": "THR", "W": "TRP", "Y": "TYR", "V": "VAL",
}
pdb_lines = []
atom_idx = 1
for res_idx, aa in enumerate(sequence):
res_name = aa_3letter.get(aa, "UNK")
for atom_i, atom_name in enumerate(atom_names):
if atom_i >= positions.shape[1]:
break
x, y, z = positions[res_idx, atom_i]
if any(math.isnan(c) for c in (x, y, z)):
continue
pdb_lines.append(
f"ATOM {atom_idx:5d} {atom_name:<3s} {res_name} A{res_idx + 1:4d} "
f"{x:8.3f}{y:8.3f}{z:8.3f} 1.00 0.00 {atom_name[0]:>2s}"
)
atom_idx += 1
pdb_lines.append("END")
return "\n".join(pdb_lines)
# ------------------------------------------------------------------
# Task 1: Load gene set
# ------------------------------------------------------------------
@cpu_env.task()
async def load_genes(
gene_set: str = "insulin",
custom_json: str = "",
) -> flyte.io.Dir:
"""Load a set of homologous genes from different species."""
if custom_json:
data = json.loads(custom_json)
elif gene_set in GENE_SETS:
data = GENE_SETS[gene_set]
else:
available = ", ".join(GENE_SETS.keys())
raise ValueError(f"Unknown gene set '{gene_set}'. Available: {available}")
log.info(f"Loaded gene set: {data['gene_name']} - {len(data['sequences'])} species")
out_dir = tempfile.mkdtemp(prefix="gene_compare_")
with open(os.path.join(out_dir, "genes.json"), "w") as f:
json.dump(data, f)
return await flyte.io.Dir.from_local(out_dir)
# ------------------------------------------------------------------
# Task 2: Score sequences with Carbon
# ------------------------------------------------------------------
@gpu_env.task(report=True)
async def score_sequences(
genes_dir: flyte.io.Dir,
model_name: str = "HuggingFaceBio/Carbon-3B",
) -> str:
"""Score each species' gene with Carbon-3B genomic language model.
Returns per-species log-likelihood scores and sequence metadata.
"""
import torch
from transformers import AutoModelForCausalLM, AutoTokenizer
log.info(f"Loading Carbon model: {model_name}")
device = "cuda" if torch.cuda.is_available() else "cpu"
tokenizer = AutoTokenizer.from_pretrained(model_name, trust_remote_code=True)
model = AutoModelForCausalLM.from_pretrained(
model_name, trust_remote_code=True,
dtype=torch.bfloat16 if device == "cuda" else torch.float32,
).to(device)
model.eval()
genes_path = await genes_dir.download()
with open(os.path.join(genes_path, "genes.json")) as f:
data = json.load(f)
species_names = list(data["sequences"].keys())
n = len(species_names)
scores = {}
for i, species in enumerate(species_names):
await flyte.report.replace.aio(_wrap_report(
f"
"
)
for species in species_names:
s = scores[species]
html_parts.append(
f'
{species}
{s["common_name"]}
'
f'
{s["length"]}bp
{s["gc_content"]:.1%}
'
f'
{s["protein_length"]}aa
'
f'
{s["log_likelihood"]:.2f}
'
)
html_parts.append("
")
html_parts.append('
')
html_parts.append(_make_bar_chart(
species_names,
{"Log-Likelihood": [scores[s]["log_likelihood"] for s in species_names]},
title="Carbon Log-Likelihood per Species",
value_format=".1f",
))
html_parts.append("
")
await flyte.report.replace.aio(_wrap_report("\n".join(html_parts)), do_flush=True)
result = {
"gene_name": data["gene_name"],
"description": data["description"],
"species": species_names,
"scores": scores,
}
return json.dumps(result)
# ------------------------------------------------------------------
# Task 3: Align sequences and compute similarity
# ------------------------------------------------------------------
@cpu_env.task(report=True)
async def align_and_compare(
scores_json: str,
genes_dir: flyte.io.Dir,
) -> str:
"""Align sequences with Needleman-Wunsch and compute pairwise identity.
Translates DNA to protein, builds DNA and protein identity matrices,
and generates phylogenetic trees from sequence divergence.
"""
scores_data = json.loads(scores_json)
species_names = scores_data["species"]
scores = scores_data["scores"]
gene_name = scores_data["gene_name"]
n = len(species_names)
genes_path = await genes_dir.download()
with open(os.path.join(genes_path, "genes.json")) as f:
data = json.load(f)
await flyte.report.replace.aio(_wrap_report(
f"
{gene_name} - Sequence Alignment
"
f"
Aligning {n} species with Needleman-Wunsch...
"
), do_flush=True)
# Pairwise DNA identity matrix
identity_matrix = []
for sp1 in species_names:
row = []
for sp2 in species_names:
dna1 = data["sequences"][sp1]["dna"]
dna2 = data["sequences"][sp2]["dna"]
identity = _sequence_identity(dna1, dna2)
row.append(round(identity, 4))
identity_matrix.append(row)
# Pairwise protein identity matrix
protein_matrix = []
for sp1 in species_names:
row = []
for sp2 in species_names:
identity = _sequence_identity(scores[sp1]["protein"], scores[sp2]["protein"])
row.append(round(identity, 4))
protein_matrix.append(row)
# Average pairwise identities (exclude diagonal)
dna_pairs = [identity_matrix[i][j] for i in range(n) for j in range(i + 1, n)]
prot_pairs = [protein_matrix[i][j] for i in range(n) for j in range(i + 1, n)]
avg_dna = sum(dna_pairs) / len(dna_pairs) if dna_pairs else 0
avg_prot = sum(prot_pairs) / len(prot_pairs) if prot_pairs else 0
# Most/least similar pair
best_pair = max(range(len(dna_pairs)), key=lambda k: dna_pairs[k])
worst_pair = min(range(len(dna_pairs)), key=lambda k: dna_pairs[k])
pair_indices = [(i, j) for i in range(n) for j in range(i + 1, n)]
best_sp = f"{species_names[pair_indices[best_pair][0]]}-{species_names[pair_indices[best_pair][1]]}"
worst_sp = f"{species_names[pair_indices[worst_pair][0]]}-{species_names[pair_indices[worst_pair][1]]}"
# Report
html_parts = [
f"
{gene_name} - Sequence Alignment
",
f'
Pairwise alignment using Needleman-Wunsch (match=2, mismatch=-1, gap=-2). '
f"Identity is computed as matches / aligned length from the optimal global alignment.
")
# DNA vs Protein conservation comparison
html_parts.append('
')
html_parts.append(_make_bar_chart(
species_names,
{
"DNA vs Human": [identity_matrix[0][j] for j in range(n)],
"Protein vs Human": [protein_matrix[0][j] for j in range(n)],
},
title=f"Conservation vs {species_names[0]} (DNA and Protein)",
value_format=".0%",
))
html_parts.append("
")
await flyte.report.replace.aio(_wrap_report("\n".join(html_parts)), do_flush=True)
result = {
"gene_name": gene_name,
"description": scores_data["description"],
"species": species_names,
"scores": scores,
"dna_identity_matrix": identity_matrix,
"protein_identity_matrix": protein_matrix,
}
return json.dumps(result)
# ------------------------------------------------------------------
# Task 4: Fold proteins with ESMFold
# ------------------------------------------------------------------
@gpu_env.task(report=True)
async def fold_proteins(
comparison_json: str,
max_length: int = 400,
) -> str:
"""Fold each species' translated protein with ESMFold for 3D comparison.
Returns PDB strings and pLDDT confidence scores for each species.
"""
import torch
import numpy as np
from transformers import AutoTokenizer, EsmForProteinFolding
comparison = json.loads(comparison_json)
species_names = comparison["species"]
scores = comparison["scores"]
log.info("Loading ESMFold model...")
device = "cuda" if torch.cuda.is_available() else "cpu"
tokenizer = AutoTokenizer.from_pretrained("facebook/esmfold_v1")
model = EsmForProteinFolding.from_pretrained("facebook/esmfold_v1", low_cpu_mem_usage=True)
model = model.to(device)
model.eval()
structure_data = {}
n = len(species_names)
for idx, species in enumerate(species_names):
protein = scores[species]["protein"]
if len(protein) > max_length:
log.info(f"Skipping {species} ({len(protein)} aa > {max_length} max)")
continue
log.info(f"ESMFold [{idx + 1}/{n}]: {species} ({len(protein)} aa)")
await flyte.report.replace.aio(_wrap_report(
f"
"
), do_flush=True)
inputs = tokenizer(protein, return_tensors="pt", add_special_tokens=False).to(device)
with torch.no_grad():
outputs = model(**inputs)
pdb_str = _outputs_to_pdb(outputs, protein)
plddt_raw = outputs.plddt[0].cpu().numpy()
if plddt_raw.ndim == 2:
plddt_raw = plddt_raw[-1]
plddt = plddt_raw.flatten()[:len(protein)]
if plddt.max() <= 1.0:
plddt = plddt * 100
plddt_mean = float(np.mean(plddt))
structure_data[species] = {
"pdb_str": pdb_str,
"plddt_mean": round(plddt_mean, 1),
"plddt_per_residue": [round(float(v), 1) for v in plddt[:len(protein)]],
"protein_length": len(protein),
}
log.info(f" → mean pLDDT: {plddt_mean:.1f}")
# Report with 3D viewers
n_folded = len(structure_data)
avg_plddt = sum(d["plddt_mean"] for d in structure_data.values()) / n_folded if n_folded else 0
threeDmol_script = ''
stats_html = f"""
ESMFold - Cross-Species Structure Comparison
ESMFold predicts 3D structure directly from amino acid sequence.
Comparing structures across species reveals which parts of the protein are
structurally conserved (functional core) vs divergent (surface loops, species-specific adaptations).
',
]
# Key metrics
# Average pairwise identity (exclude diagonal)
n = len(species)
dna_pairs = [dna_matrix[i][j] for i in range(n) for j in range(i + 1, n)]
protein_pairs = [protein_matrix[i][j] for i in range(n) for j in range(i + 1, n)]
avg_dna_id = sum(dna_pairs) / len(dna_pairs) if dna_pairs else 0
avg_protein_id = sum(protein_pairs) / len(protein_pairs) if protein_pairs else 0
avg_plddt = sum(d["plddt_mean"] for d in structures.values()) / len(structures) if structures else 0
html_parts.append('
')
html_parts.append(f'
{n}
Species
')
html_parts.append(f'
{avg_dna_id:.0%}
Avg DNA Identity
')
html_parts.append(f'
{avg_protein_id:.0%}
Avg Protein Identity
')
html_parts.append(f'
{avg_plddt:.1f}
Avg pLDDT
')
html_parts.append(f'
{len(structures)}
Structures Folded
')
html_parts.append("
")
# Full comparison table
html_parts.append("
Per-Species Detail
")
html_parts.append(
"
Species
Scientific Name
DNA (bp)
"
"
Protein (aa)
GC%
Carbon LL
pLDDT
"
)
for sp in species:
s = scores[sp]
plddt = structures.get(sp, {}).get("plddt_mean", "N/A")
plddt_str = f"{plddt:.1f}" if isinstance(plddt, float) else plddt
html_parts.append(
f'
{sp}
{s["common_name"]}
'
f'
{s["length"]}
{s["protein_length"]}
'
f'
{s["gc_content"]:.1%}
{s["log_likelihood"]:.2f}
'
f'
{plddt_str}
'
)
html_parts.append("
")
# GC content comparison
html_parts.append('
')
html_parts.append(_make_bar_chart(
species,
{"GC Content": [scores[s]["gc_content"] for s in species]},
title="GC Content Across Species",
value_format=".2f",
))
html_parts.append("
")
# DNA phylogenetic tree
html_parts.append("
Phylogenetic Relationships
")
html_parts.append(
'
'
"Trees built from pairwise sequence identity using UPGMA clustering. "
"Species that diverged more recently cluster together. DNA and protein trees "
"may differ when synonymous mutations dominate."
"
All 4 pipeline stages completed. View individual task reports for DNA/protein
identity heatmaps, phylogenetic trees, interactive 3D protein structures with
pLDDT confidence, Carbon log-likelihood scores, and evolutionary analysis.
"""
await flyte.report.replace.aio(_wrap_report(final_html), do_flush=True)
log.info("Pipeline complete.")
return comparison_json, structures_json
# {{/docs-fragment pipeline}}
if __name__ == "__main__":
flyte.init_from_config()
run = flyte.run(pipeline)
print(run.url)
run.wait()
```
*Source: https://github.com/unionai/unionai-examples/blob/main/v2/tutorials/genomic_gene_comparison/genomic_gene_comparison.py*
## Run the workflow
From the [example directory](https://github.com/unionai/unionai-examples/tree/main/v2/tutorials/genomic_gene_comparison):
```
cd v2/tutorials/genomic_gene_comparison
uv run --script genomic_gene_comparison.py
```
Or submit a specific gene set with the Flyte CLI:
```
flyte run genomic_gene_comparison.py pipeline --gene_set hemoglobin
```
This example needs a GPU for Carbon and ESMFold. Open the run URL and check each task's report tab for heatmaps, dendrograms, and interactive 3D viewers.
---
**Source**: https://github.com/unionai/unionai-docs/blob/main/content/tutorials/biotech-healthcare/genomic-gene-comparison/_index.md
**HTML**: https://www.union.ai/docs/v2/union/tutorials/biotech-healthcare/genomic-gene-comparison/