Sequence alignment using PyMOL

align.py

#!/usr/bin/env python

# Sequence alignment using PyMOL
# The purpose of this script is to generate a sequence alignment between
# the original crystal structure of the apo and holo models, and the sequence
# of the finalised, ungapped Rosetta models. This allows us to get a 1 to 1
# corresponcence between the residue numberings in both structures.

# USAGE: Run once from the project root.
# "pockets.csv" contains the information about apo holo pairs.

import __main__
import glob
import os
import pandas as pd
import pymol
import time


def align(mobile, reference, filename="alignment.aln"):
    """Do a sequence alignment between two pdbs and write the output.
    Parameters
    ----------
    mobile : A PDB file. Sequence will be first in the alignment.
    reference : A PDB file. Sequence will be second in alignment.
    filename : Output name of alignment file.
    """
    __main__.pymol_argv = ['pymol', '-qc']
    pymol.finish_launching()
    pymol.cmd.load(mobile)
    pymol.cmd.load(reference)
    obj = list(pymol.cmd.get_object_list('all'))
    pymol.cmd.align(obj[0], obj[1], object='aln', transform=0)
    pymol.cmd.save(filename, 'aln')
    pymol.cmd.reinitialize()
    time.sleep(1)


if __name__ == "__main__":
    # Read data frame with holo, apo pairs
    cols = ["APO", "HOLO", "ligand", "size", "APO_accession",
            "HOLO_accession", "same"]
    pairs = pd.read_csv("pockets.csv", index_col=None, names=cols)

    for index, row in pairs.iterrows():
        apo = row['APO']
        holo = row['HOLO']
        if os.path.exists(apo + "/selected_models") and os.path.exists(
                holo + "/selected_models"):
            print("Aligning:", apo, holo)
            apo_model = glob.glob(apo + "/selected_models/S_????*.pdb")[0]
            holo_model = glob.glob(holo + "/selected_models/S_????*.pdb")[0]
            apo_crystal = "{}/{}.pdb".format(apo, apo[:-2])
            holo_crystal = "{}/{}.pdb".format(holo, holo[:-2])

            align(apo_crystal, apo_model, filename=apo +
                  "/apo_crystal.aln")
            align(holo_crystal, apo_model, filename=apo +
                  "/holo_crystal.aln")
            align(apo_crystal, holo_model, filename=holo +
                  "/apo_crystal.aln")
            align(holo_crystal, holo_model, filename=holo +
                  "/holo_crystal.aln")