Ciencias de la vida: biología y medicina

from Bio.Seq import Seq
from Bio.SeqUtils import gc_fraction
from Bio.Data import CodonTable
import random

random.seed(7)
protein = "MALWMRLLPLLALLALWGPDPAAA"      # un tramo tipo péptido señal, 24 aminoácidos
table = CodonTable.unambiguous_dna_by_id[1]   # el código genético estándar
codons_for = {}
for codon, aa in table.forward_table.items():
    codons_for.setdefault(aa, []).append(codon)

gene = Seq("".join(random.choice(codons_for[aa]) for aa in protein) + "TAA")   # retrotraducir + codón de parada
print("ADN      :", gene)
print("ARN      :", gene.transcribe())
print("proteína :", gene.translate(to_stop=True), "| coincide con el diseño:", str(gene.translate(to_stop=True)) == protein)
print("contenido GC:", round(gc_fraction(gene), 3), "| complementaria inversa:", gene.reverse_complement()[:24], "...")

# la misma proteína puede escribirse con muchas secuencias de ADN: el código genético es redundante
print({aa: len(c) for aa, c in sorted(codons_for.items())})

# buscar marcos de lectura abiertos (ATG ... parada) en un tramo de ADN que esconde nuestro gen
noise = lambda n: "".join(random.choice("ACGT") for _ in range(n))
genome = Seq(noise(120) + str(gene) + noise(150))
orfs = []
for frame in range(3):
    prot = str(genome[frame:].translate())
    pos = 0
    for chunk in prot.split("*"):
        start = chunk.find("M")
        if start != -1 and len(chunk) - start >= 20:
            orfs.append((frame, frame + 3 * (pos + start), chunk[start:]))
        pos += len(chunk) + 1
print("ORF de al menos 20 aa:", [(f, s, p) for f, s, p in orfs])

from Bio import Align

def mutate(seq: str, n_subs: int, n_indels: int, rng: random.Random) -> str:
    s = list(seq)
    for _ in range(n_subs):
        i = rng.randrange(len(s)); s[i] = rng.choice([b for b in "ACGT" if b != s[i]])
    for _ in range(n_indels):
        i = rng.randrange(len(s))
        if rng.random() < 0.5: del s[i]
        else: s.insert(i, rng.choice("ACGT"))
    return "".join(s)

rng = random.Random(1)
ancestor = noise(300)
# un árbol que elegimos: ((A,B),(C,(D,E))) con ramas de distinta longitud
ab = mutate(ancestor, 18, 1, rng);  cde = mutate(ancestor, 24, 1, rng)
de = mutate(cde, 12, 1, rng)
species = {"A": mutate(ab, 6, 0, rng), "B": mutate(ab, 8, 1, rng), "C": mutate(cde, 10, 1, rng),
           "D": mutate(de, 5, 0, rng), "E": mutate(de, 7, 1, rng)}

aligner = Align.PairwiseAligner(mode="global", match_score=2, mismatch_score=-1,
                                open_gap_score=-4, extend_gap_score=-1)
best = aligner.align(species["A"], species["B"])[0]
print(f"puntuación {best.score:.0f}")
print(str(best)[:300])

from Bio.Phylo.TreeConstruction import DistanceMatrix, DistanceTreeConstructor
from Bio import Phylo
import matplotlib.pyplot as plt

names = list(species)
def p_distance(a, b):
    aln = aligner.align(a, b)[0]
    s1, s2 = aln[0], aln[1]
    pairs = [(x, y) for x, y in zip(s1, s2) if x != "-" and y != "-"]
    return sum(x != y for x, y in pairs) / len(pairs)

matrix = [[p_distance(species[names[i]], species[names[j]]) for j in range(i)] + [0.0] for i in range(len(names))]
dm = DistanceMatrix(names, matrix)
print(dm)
tree = DistanceTreeConstructor().nj(dm)
tree.root_at_midpoint()
fig, ax = plt.subplots(figsize=(6, 3.2))
Phylo.draw(tree, axes=ax, do_show=False); ax.set_title("Árbol por unión de vecinos"); plt.tight_layout(); plt.show()

from lifelines.datasets import load_gbsg2
from lifelines import KaplanMeierFitter, CoxPHFitter
from lifelines.statistics import logrank_test

gb = load_gbsg2()
print(gb.shape); print(gb["horTh"].value_counts().to_dict(), "| eventos (recaídas):", int(gb["cens"].sum()))
gb.head()

fig, ax = plt.subplots(figsize=(7, 4))
for group, d in gb.groupby("horTh"):
    KaplanMeierFitter(label=f"terapia hormonal: {group}").fit(d["time"] / 365.25, d["cens"]).plot_survival_function(ax=ax)
ax.set_xlabel("años tras la cirugía"); ax.set_ylabel("probabilidad sin recaída"); ax.set_title("Estimación de Kaplan-Meier")
plt.tight_layout(); plt.show()

yes, no = gb[gb.horTh == "yes"], gb[gb.horTh == "no"]
test = logrank_test(yes["time"], no["time"], event_observed_A=yes["cens"], event_observed_B=no["cens"])
print(f"valor p de la prueba log-rank: {test.p_value:.4f}")

X = gb.assign(horTh=(gb.horTh == "yes").astype(int), menostat=(gb.menostat == "Post").astype(int),
              tgrade=gb.tgrade.map({"I": 1, "II": 2, "III": 3}))
cph = CoxPHFitter().fit(X[["time", "cens", "horTh", "age", "menostat", "tsize", "tgrade", "pnodes", "progrec", "estrec"]],
                        duration_col="time", event_col="cens")
cph.print_summary(columns=["exp(coef)", "exp(coef) lower 95%", "exp(coef) upper 95%", "p"], decimals=3)

fig, ax = plt.subplots(figsize=(6, 3.6))
cph.plot(hazard_ratios=True, ax=ax); ax.set_title("Cocientes de riesgos (IC 95%)"); plt.tight_layout(); plt.show()
hr = cph.hazard_ratios_["horTh"]
print(f"Ajustando por los demás factores, la terapia hormonal multiplica el riesgo de recaída por {hr:.2f} "
      f"(alrededor de un {100 * (1 - hr):.0f}% menos). Más ganglios positivos y un grado tumoral mayor lo aumentan.")

from sklearn.datasets import load_breast_cancer
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_predict, StratifiedKFold
from sklearn.metrics import roc_auc_score, roc_curve, confusion_matrix

data = load_breast_cancer()
y = (data.target == 0).astype(int)                  # 1 = maligno
model = make_pipeline(StandardScaler(), LogisticRegression(max_iter=5000, C=0.5))
proba = cross_val_predict(model, data.data, y, cv=StratifiedKFold(5, shuffle=True, random_state=0), method="predict_proba")[:, 1]
print(f"AUC ROC con validación cruzada: {roc_auc_score(y, proba):.3f}")

fpr, tpr, thr = roc_curve(y, proba)
for target_sens in (0.95, 0.99):
    i = next(k for k, t in enumerate(tpr) if t >= target_sens)
    tn, fp, fn, tp = confusion_matrix(y, proba >= thr[i]).ravel()
    print(f"umbral {thr[i]:.3f}: sensibilidad {tp / (tp + fn):.3f}, especificidad {tn / (tn + fp):.3f}, "
          f"malignos no detectados {fn}, falsas alarmas {fp}")

plt.figure(figsize=(4.5, 4)); plt.plot(fpr, tpr); plt.plot([0, 1], [0, 1], "--", color="grey")
plt.xlabel("1 - especificidad"); plt.ylabel("sensibilidad"); plt.title("Curva ROC (validación cruzada de 5 particiones)"); plt.tight_layout(); plt.show()