Localizar una secuencia de aminoácidos en una estructura (PyMol).


El siguiente script de Python se ejecuta desde PyMol y busca en la cadena "C" del objeto "7VDP" la secuencia de aminoácidos "KEAFWDRCLSVINLMSSKMLQINA" y si la encuentra cree un nuevo objeto llamado "peptido-P5". Logicamente se puede cambiar la cadena, el objeto, la secuencia y el nombre del objeto nuevo con el peptido que encuentre. Se permite hasta 6 errores en la secuencia (mismatches_max=6).



from pymol import cmd

# Diccionario 3→1
aa_3to1 = {
    "ALA":"A","ARG":"R","ASN":"N","ASP":"D","CYS":"C",
    "GLN":"Q","GLU":"E","GLY":"G","HIS":"H","ILE":"I",
    "LEU":"L","LYS":"K","MET":"M","PHE":"F","PRO":"P",
    "SER":"S","THR":"T","TRP":"W","TYR":"Y","VAL":"V"
}

def hamming(seq1, seq2):
    return sum(a != b for a, b in zip(seq1, seq2))

def crear_peptido_desde_secuencia(objeto="7VDP", cadena="C",
                                 secuencia_objetivo="KEAFWDRCLSVINLMSSKMLQINA",
                                 mismatches_max=6,
                                 nombre_objeto="peptido-P5"):

    modelo = cmd.get_model(f"{objeto} and chain {cadena} and name CA")
    
    residuos = []
    secuencia_real = ""
    
    for atom in modelo.atom:
        resi = atom.resi
        resn = atom.resn
        
        if len(residuos) == 0 or residuos[-1][0] != resi:
            residuos.append((resi, resn))
            secuencia_real += aa_3to1.get(resn, "X")
    
    L = len(secuencia_objetivo)
    mejor_match = None
    mejor_score = L + 1  # peor caso
    
    # Sliding window
    for i in range(len(secuencia_real) - L + 1):
        ventana = secuencia_real[i:i+L]
        dist = hamming(ventana, secuencia_objetivo)
        
        if dist < mejor_score:
            mejor_score = dist
            mejor_match = i
    
    if mejor_match is None or mejor_score > mismatches_max:
        print("No se encontró ninguna región dentro del umbral de mismatches.")
        print(f"Mejor coincidencia tenía {mejor_score} mismatches.")
        return
    
    print(f"Mejor match en posición: {mejor_match}")
    print(f"Mismatches: {mejor_score} / {L}")
    
    inicio = mejor_match
    fin = mejor_match + L
    
    resis_seleccionados = [residuos[i][0] for i in range(inicio, fin)]
    
    seleccion = f"{objeto} and chain {cadena} and resi " + "+".join(resis_seleccionados)
    
    cmd.create(nombre_objeto, seleccion)
    
    # Información extra útil
    seq_encontrada = secuencia_real[inicio:fin]
    print("Secuencia encontrada:")
    print(seq_encontrada)
    print("Secuencia objetivo:")
    print(secuencia_objetivo)

# Ejecutar
crear_peptido_desde_secuencia()