<jupyter_start><jupyter_text>Searching PDB for binders to desired proteinIdea: - complete partial complexes- start with single chain/ sequence- find homologs which are in a complex- check if complex could also be formed with query instead of the homolog- use as template for homology modellingDan Farrell & <NAME><jupyter_code>from pyspark import SparkContext
from pyspark.sql import SparkSession
from mmtfPyspark.io import mmtfReader
from mmtfPyspark.webfilters import SequenceSimilarity
from mmtfPyspark.mappers import StructureToPolymerChains<jupyter_output><empty_output><jupyter_text>__Configure Spark__<jupyter_code>spark = SparkSession.builder.master("local[*]").appName("Search").getOrCreate()
sc = spark.sparkContext<jupyter_output><empty_output><jupyter_text>Prepare an example queryDownload a structure and extract its sequence for the query.<jupyter_code>query_id = '4LGR.A'
query_tf = mmtfReader.download_reduced_mmtf_files([query_id[:4]], sc)
primary = query_tf.flatMap(StructureToPolymerChains())
query_seq = primary.lookup(query_id)[0].entity_list[0]['sequence']<jupyter_output><empty_output><jupyter_text>Run BLAST against PDBLoad the PDB and use the query sequence on it with BLAST.<jupyter_code>path = "../mmtf-workshop-2018/resources/mmtf_reduced_sample"
pdb = mmtfReader.read_sequence_file(path, sc).cache()
pdb = pdb.flatMap(StructureToPolymerChains())\
      .filter(SequenceSimilarity(sequence=query_seq, searchTool=SequenceSimilarity.BLAST,
                                 eValueCutoff=0.1, sequenceIdentityCutoff=10,
                                 maskLowComplexity=True)).collect()<jupyter_output><empty_output><jupyter_text>Print out hits:<jupyter_code>for pdbId, structure in pdb:
    print(f"{pdbId} :     {structure.entity_list[0]['sequence']}")<jupyter_output>1RL0.A :     ATAYTLNLANPSASQYSSFLDQIRNNVRDTSLIYGGTDVAVIGAPSTTDKFLRLNFQGPRGTVSLGLRRENLYVVAYLAMDNANVNRAYYFKNQITSAELTALFPEVVVANQKQLEYGEDYQAIEKNAKITTGDQSRKELGLGINLLITMIDGVNKKVRVVKDEARFLLIAIQMTAEAARFRYIQNLVTKNFPNKFDSENKVIQFQVSWSKISTAIFGDCKNGVFNKDYDFGFGKVRQAKDLQMGLLKYLGRPKS
3BWH.A :     NVRFDLSSATSSSYKTFIKNLREALPKDGKVYDIPVLLSTVMDSRRFILIDLVNYDGQSITAAIDVLNVYIVAYSTGTVSYFFQQVPAQAPKLLFKGTQQRTLPYTGNYENLQTAAKKLRENIELGLPALDSAITTLFHYNAEAAASALLVLIQTTSEAARFRYIELQIANNVGTKFKPSQTIISLENNWSALSKQIQIAKNKNGQFETPVILIDPQGNRVQITNVTSNVVTQNIQLLLNIGATA
3U6T.A :     DVSFRLSGADPSSYGMFIKDLRNALPHTEKVYNIPLLLPSVSGAGRYLLMHLFNYDGNTITVAVDVTNVYIMGYLALTTSYFFNEPAADLASQYVFRSARRKITLPYSGNYERLQIAAGKPREKIPIGLPALDTAISTLLHYDSTAAAGALLVLIQTTAEAARFKYIEQQIQERAYRDEVPSSATISLENSWSGLSKQIQLAQGNNGVFRTPTVLVDSKGNRVQITNVTSNVVTSNIQLLLNTKNI
4M1U.A :     REFTIDFSTQQSYVSSLNSIRTEISTPLEHISQGTTSVSVINHTPPGSYFAVDIRGLDVYQARFDHLRLIIEQNNLYVAGFVNTATNTFYRFSDFTHISVPGVTTVSMTTDSSYTTLQRVAALERSGMQISRHSLVSSYLALMEFSGNTMTRDASRAVLRFVTVTAEALRFRQIQREFRQALSETAPVYTMTPGDVDLT[...]<jupyter_text>TODO: Find homologs which are oligomers TODO: Extract cahins bound to homologs Part 2: General idea:1. Align your possible binder to the query protein's ideal binder.    * Do sequence alignment between the two    * Do structural alignment based on the sequence alignment2. Combine your protein of interest with a possible binder, then combine them into one structure3. Generate a matrix around the interface for **both** (anchor on the constant protein) and fill the matrix with vdw spheres that represent the CAlphas.    * **both**= input structure + ideal binder && input structure + possible binder4. Then compare the two matricies and give it a score based on the cross correlation5. Sort by score and see if you have found any unexpected binders!<jupyter_code>import os
import Bio.PDB
from Bio.PDB import PDBIO
from Bio.Align.Applications import ClustalOmegaCommandline
from Bio.PDB.mmtf import MMTFParser
from Bio.PDB.Polypeptide import is_aa, three_to_one
from Bio import AlignIO
import py3Dmol
def get_sequence_from_biopython_and_chain_name(bp_structure, chain_name):
    """Get the AA sequence of a chain from biopython

    Args:
        bp_structure: A biopython structure
        chain_name: A chain in that biopython structure to get the AA sequence of

    Returns:
        sequence

    Notes:
        Some horrible pdbs only have one atom for whole residues (IE just N, instead of NCC)
        so you MUST check for CA since that's what we're aligning on
    """
    sequence = ""
    for chain in bp_structure:
        if chain.id != chain_name:
            continue
        for res in chain:
            if is_aa(res):
                for atom in res:
                    if atom.name == "CA":
                        sequence += three_to_one(res.get_resname())
    return sequence

def get_aln_atoms(bp_structure, residues, chain_name):
    """Get CA atoms from aligned residues

    Alignments are residue based, but we want to do our alignment based on CA atoms only
    essentially convert residue to atom['CA']

    Args:
        bp_structure: a biopython structure
        residues: the index of aligned residues that you're interested in
        chain_name: The chain name of interest, (tells us what chains to iterate over)
    
    Retruns:
        biopython CA atoms that we want to use in the alignment
    """
    ret_atoms = []
    residue_count = 0
    for chain in bp_structure:
        if chain.id != chain_name:
            continue
        for res in chain:
            if is_aa(res):
                for atom in res:
                    if atom.name == "CA" and residue_count in residues:
                        ret_atoms.append(atom)
                residue_count += 1
    return ret_atoms

def make_combo_pose(input_chainID, target_chain, ref_model, target_model):
    """Combine the best parts of our target/reference

    Args:
        input_chainID: The original quereied seuqence's chainID
        target_chain: The chainID of the target
        ref_model: The biopython model of the reference(input)
        target_model: the biopython model of the target pdb

    Returns:
        a biopython structure containing the reference_model(input_chainID)
    """
    base_structure = Bio.PDB.Structure.Structure("combo_pose")
    base_model = Bio.PDB.Model.Model(0)
    base_structure.add(base_model)
    ref_chains_to_mv = []
    for i, chain in enumerate(ref_model):
        if chain.id != input_chainID:
            continue
        ref_chains_to_mv.append(i)
    ref_chains = [x for x in ref_model]
    for x in reversed(ref_chains_to_mv):
        ref_chains[x].detach_parent()
        ref_chains[x].id = 'A'
        base_structure[0].add(ref_chains[x])

    target_chains_to_mv = []
    for i, chain in enumerate(target_model):
        if chain.id != target_chain:
            continue
        target_chains_to_mv.append(i)
    target_chains = [x for x in target_model]
    for x in reversed(target_chains_to_mv):
        target_chains[x].detach_parent()
        target_chains[x].id = 'B'
        base_structure[0].add(target_chains[x])

    io = Bio.PDB.PDBIO()
    io.set_structure(base_structure)
    io.save("test_align.pdb")
    return base_structure

def align_two_structures(input_pdbID_chainID, reference, target):
    mmtf_parser = MMTFParser()
    # Only get first model
    ref_structure = mmtf_parser.get_structure_from_url(reference.split('_')[0])
    ref_model = ref_structure[0]
    ref_sequence = get_sequence_from_biopython_and_chain_name(ref_model, reference.split('_')[1])
    target_structure = mmtf_parser.get_structure_from_url(target.split('_')[0])
    target_model = target_structure[0]
    target_sequence = get_sequence_from_biopython_and_chain_name(target_model, target.split('_')[1])
    if os.path.isfile("out.clust.aln"):
        os.remove("out.clust.aln")
    with open("ref_targ.fasta", 'w') as f:
        f.write(">reference\n{}\n".format(ref_sequence))
        f.write(">target\n{}\n".format(target_sequence))
    clustalomega_cline = ClustalOmegaCommandline(infile="ref_targ.fasta", outfile="out.clust.aln",
                                                 verbose=True, auto=True)
    clustalomega_cline()
    align = AlignIO.read("out.clust.aln", "fasta")


    # We look for aligned residues in the clustalO alignment, and mark them to be used in super imposition
    ref_atoms = []
    ref_count = 0
    target_atoms = []
    target_count = 0
    for i, _ in enumerate(align[0].seq):
        use = True
        # Reference on top
        if align[0].seq[i] == '-':
            use = False
        else:
            ref_count += 1
        if align[1].seq[i] == '-':
            use = False
        else:
            target_count += 1
        if use:
            ref_atoms.append(ref_count-1)
            target_atoms.append(target_count-1)
    ref_atoms_bp = get_aln_atoms(ref_model, ref_atoms, reference.split('_')[1])
    target_atoms_bp = get_aln_atoms(target_model, target_atoms, target.split('_')[1])

    super_imposer = Bio.PDB.Superimposer()
    super_imposer.set_atoms(ref_atoms_bp, target_atoms_bp)
    super_imposer.apply(target_model.get_atoms())

    combo_pose = make_combo_pose(input_pdbID_chainID.split('_')[-1], target.split('_')[-1], ref_model, target_model)[0]
    # A = reference protein, B = target protein

def part_2_main(inputs):
    """Run alignment and scoring
    
    Args:
        inputs: a list of tuples, each tuple has the format
        (input pdbID_chainID of base protein of interest,
            found pdbID_chainID of protein that interacts with tuple[0],
            protein that is homologous, or structurally similar to tuple[1])

    Notes:
        General Idea:
            combine your protein of interest with a possible binder, then combine them into one,
            generate a matrix around the interface for both (anchor on the constant protein) and fill the matrix
            with vdw spheres that represent the CAlphas.  then compare that matrix to the matrix of the ideal
            binder, and give it a score based on the cross correlation

    TODO:
        SCORING based on vdw spheres cross correlation (overlap)
        of result and input 'answer'
    """
    for pair in inputs:
        align_two_structures(pair[0], pair[1], pair[2])

inputs = [("2BJN_A", "3KXC_A", "3KXC_C")]
# part_2_main(inputs)
viewer = py3Dmol.view(query='pdb:2BJN')

# setting styles will be covered in the next tutorial
viewer.setStyle({'cartoon': {'color': 'spectrum'}})
viewer.setStyle({'hetflag': True}, {'stick':{'radius': 0.3, 'singleBond': False}})
viewer.zoomTo()
viewer.show()

viewer = py3Dmol.view(query='pdb:3KXC')

# setting styles will be covered in the next tutorial
viewer.setStyle({'cartoon': {'color': 'spectrum'}})
viewer.setStyle({'hetflag': True}, {'stick':{'radius': 0.3, 'singleBond': False}})
viewer.zoomTo()
viewer.show()

structure = open('test_align.pdb','r').read()
viewer = py3Dmol.view()
viewer.addModel(structure,'pdb')

viewer.setStyle({'cartoon': {'color': 'spectrum'}})
viewer.setStyle({'hetflag': True}, {'stick':{'radius': 0.3, 'singleBond': False}})
viewer.zoomTo()
viewer.show()
spark.stop()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08779910836138659, 'fraction_numerical': 0.006641797834591939, 'mean_word_length': 3.7874564459930316, 'pattern_counts': {'":': 1, '<': 25, '<?xml version=': 0, '>': 27, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '6', 'dir': 'jupyter-structured-clean-dedup', 'id': '397618', 'license': 'Apache-2.0', 'n_tokens_mistral': 3830, 'n_tokens_neox': 3621, 'n_words': 873, 'path': 'finding_protein_complexes_from_single_sequence/Searching_PDB_for_binders_to_desired_proteins.ipynb', 'repo_name': 'danpf/mmtf-hackathon-2018', 'size': '10999'}