14772793414

Committed 01 May 2025 08:56AM UTC coverage: 53.774% (-0.04%) from 53.813%

Build # 14772793414

Build Type

Pull #301

github

Committed by

web-flow

Commit Message

Merge d3b9b630f into eec26f4c5

Pull Request Pull Request #301: Issue 257

Run Details

4 of 5 new or added lines in 1 file covered. (80.0%)

1 existing line in 1 file now uncovered.

1845 of 3431 relevant lines covered (53.77%)

0.54 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

61.45

/PANDORA/PMHC/Model.py

from pyexpat import model
import os

from Bio import pairwise2
from Bio.SeqUtils import seq1
from Bio.PDB import PDBParser, PDBIO, Select

import PANDORA

class Model:

    def __init__(self, target, model_path='', output_dir=False, pdb=False, molpdf=0, dope=0):
        ''' Initiate model object
        Args:
            target: Target object
            model_path: (string) path to hypothetical model
            output_dir: (string) output directory
            pdb:  Bio.PDB object of the hypothetical model
            molpdf: (float) molpdf score
            dope:  (float) DOPE score
        '''

        self.target = target
        self.model_path = model_path
        self.molpdf = molpdf
        self.dope = dope

        # Define the output directory
        if output_dir == False:
            self.output_dir = f"{os.getcwd()}/{self.target.id}"
        else:
            self.output_dir = output_dir
        


        # Check if the user gave either the path to the model pdb or the pdb itself.
        if self.model_path == '' and not pdb:
            raise Exception('Provide the path to a model structure or a Bio.PDB object')
        # If there is a model path and no pdb, parse the pdb structure from that path.
        if not pdb:
            self.pdb = PDBParser(QUIET=True).get_structure(self.target.id, self.model_path)


    def calc_LRMSD(self, reference_pdb, atoms = ['C', 'CA', 'N', 'O'], ligand_zone="whole"):
        """Calculate the L-RMSD between the decoy and reference structure (ground truth).
            This function requires the pdb2sql module for L-RMSD calculation.

        Args:
            reference_pdb: Bio.PDB object or path to pdb file
            atoms (list, optional): The list of atoms of the ligand selected to calculte the LRMSD . Defaults to ['C', 'CA', 'N', 'O'].
            ligand_zone (str, optional): The region of the Ligand selected to calculate the LRMSD. Defaults to "whole".

        Raises:
            Exception: PDB2SQL LRMSD claulation failed
        Returns: (float) L-RMSD calculated by PDB2SQL
        """        

        #from pdb2sql import pdb2sql, superpose, StructureSimilarity
        from pdb2sql import StructureSimilarity

        # load target pdb
        if isinstance(reference_pdb, str):  # if its a string, it should be the path of the pdb, then load pdb first
            ref = PDBParser(QUIET=True).get_structure(self.target.id, reference_pdb)
        else:
            ref = reference_pdb

        # Define file names as variables
        model_name = self.model_path.split('/')[-1].split('.')[1]
        #decoy_path = '%s/%s_decoy.pdb' % (self.output_dir, model_name)
        #ref_path = '%s/%s_ref.pdb' % (self.output_dir, model_name)

        # Define zones to align
        #M_lzone = list(range(4,73))
        #N_lzone = list(range(10,80))

        ## pdb2sql needs 1 big chain and 1 ligand chain with correct numbering, for MHCII, this means merging the chains.       
        #if os.path.exists(decoy_path)==False and ligand_zone.lowe()=="whole":
        if  ligand_zone.lower()=="whole":
            decoy_path, ref_path = homogenize_pdbs(self.pdb, ref, atoms, self.output_dir, model_name)
        ## calculates core LRMSD    
        elif ligand_zone.lower() == "core":
            decoy_path, ref_path = homogenize_pdbs(self.pdb, ref, atoms, self.output_dir, model_name, anchors = self.target.anchors)
        elif ligand_zone.lower() == "flanking":
            decoy_path, ref_path = homogenize_pdbs(self.pdb, ref, atoms, self.output_dir,model_name, anchors = self.target.anchors, flanking=True)

       
        start_dir = os.getcwd()
        os.chdir(self.output_dir)
        # Produce lzone file for the l-rmsd calculation
        #lzone = get_Gdomain_lzone(ref_path, self.output_dir, self.target.MHC_class)
        #TODO: check if it's MHC I or II and adapt for chain M and N
        # Get decoy structure to superpose
        #decoy_db = pdb2sql(decoy_path)
        #decoy_lzone = np.asarray(decoy_db.get('x,y,z', resSeq=M_lzone))

        # Get ref structure to superpose
        #ref_db = pdb2sql(ref_path)
        #ref_lzone = np.asarray(ref_db.get('x,y,z', resSeq=M_lzone))

        # Align the G domains
        #superpose.superpose_selection()

        try:
            # Calculate l-rmsd between decoy and reference with pdb2sql
            sim = StructureSimilarity(decoy_path, ref_path)
            #self.lrmsd = sim.compute_lrmsd_fast(method='svd', name=atoms, lzone = lzone)
            if ligand_zone.lower()=="whole":
                self.lrmsd = sim.compute_lrmsd_pdb2sql(exportpath=None, method='svd', name = atoms)
            ## calculates LRMSD for peptide core region
            elif ligand_zone.lower() == "core":
                self.core_lrmsd = sim.compute_lrmsd_pdb2sql(exportpath=None, method='svd', name=atoms)
            ## calculates RMSD for peptide flanking region
            elif ligand_zone.lower() == "flanking":
                self.flanking_lrmsd = sim.compute_lrmsd_pdb2sql(exportpath=None, method='svd', name=atoms)
        except:
            print('An error occurred while calculating the rmsd for target %s, model %s for %s' %(ref_path, decoy_path, ligand_zone))
            raise Exception('Please check your model and ref info for model %s' %self.model_path)

        # remove intermediate files
        os.system('rm %s %s' %(decoy_path, ref_path))
        os.chdir(start_dir)

def merge_chains(pdb):
    ''' Merges two chains of MHCII to one chain. pdb2sql can only calculate L-rmsd with one chain.
    Args:
        pdb: Bio.PDB object
    Returns: Bio.PDB object with its M and N chain merged as M chain
    '''
    # Merge chains
    if 'N' in [chain.id for chain in pdb.get_chains()]:

        for j in pdb[0]['N'].get_residues():
            j.id = (j.id[0], j.id[1], 'M')
            pdb[0]['M'].add(j)

        for i in pdb.get_chains():
            for model in pdb:
                for chain in model:
                    if chain.id in ['N']:
                        model.detach_child(chain.id)
    return pdb


def renumber(pdb_ref, pdb_decoy, custom_map={"MSE":"M"}):
    ''' aligns two pdb's and renumber them accordingly.
    Args:
        pdb_ref:   Bio.PDB object
        pdb_decoy: Bio.PDB object
        
    Returns: Bio.PDB objects with renumbered residues
    '''
    ref_sequences = [[chain.id, seq1(''.join([res.resname for res in chain]), custom_map=custom_map)] for chain in pdb_ref.get_chains()]
                    #[[chain.id, ('').join([seq1(res.resname) for res in chain])]
                    #  for chain in pdb_ref.get_chains()]
    ref_sequences.sort()
    decoy_sequences = [[chain.id, seq1(''.join([res.resname for res in chain]), custom_map=custom_map)] for chain in pdb_decoy.get_chains()]
                        #[[chain.id, ('').join([seq1(res.resname) for res in chain])]
                      #for chain in pdb_decoy.get_chains()]
                        
    decoy_sequences.sort()
    
    assert(len(ref_sequences) == len(decoy_sequences))


    for ind in range(len(ref_sequences)):
        pair = pairwise2.align.globalxx(ref_sequences[ind][1], decoy_sequences[ind][1])[0]
        ref_sequences[ind][1]   = pair.seqA
        decoy_sequences[ind][1] = pair.seqB
        
    ref_sequences = [[seq[0],[i+1 for i,res in enumerate(seq[1]) if res != '-']] for seq in ref_sequences]
    decoy_sequences = [[seq[0],[i+1 for i,res in enumerate(seq[1]) if res != '-']] for seq in decoy_sequences]



    def assign(pdb, pdb_sequences):
        ''' Renumbers the pdb using aligned sequences. 
        Args:
            pdb_ref:   Bio.PDB object
            pdb_decoy: Bio.PDB object
            
        Returns: Bio.PDB objects with renumbered residues
        '''
        
        
        for chain in pdb.get_chains():
            for seq in pdb_sequences:
                if chain.id == seq[0]:
                    for ind, res in enumerate(chain):
                        res.id = ('X', seq[1][ind], res.id[2])
        for chain in pdb.get_chains():
            for res in chain:
                res.id = (' ', res.id[1], ' ')
        return pdb

    pdb_ref = assign(pdb_ref, ref_sequences)
    pdb_decoy = assign(pdb_decoy, decoy_sequences)

    return pdb_ref, pdb_decoy , ref_sequences, decoy_sequences

def homogenize_pdbs(decoy, ref, atoms, output_dir,  target_id = 'MHC' , anchors =False, flanking=False):
    ''' Make sure that the decoy and reference structure have the same structure sequences.
    Args:
        decoy: Bio.PDB object of the decoy structure
        ref: Bio.PDB object of the reference structure
        output_dir: (string) directory that is used to write intermediate files
    Returns: (tuple) Bio.PDB objects with the same structure sequence
    '''
    ref_p_len=len(ref[0]['P'])
    decoy_p_len=len(decoy[0]['P'])
    

    # If you give the anchors, the core L-RMSD will be calculated.
    # The peptide residues before and after the first and last anchor residue will be discarded.
    if anchors and not flanking:
        for x in range(len(decoy[0]['P'])):
            for i in decoy[0]['P']:
                if i.id[1] < anchors[0] or i.id[1] > anchors[-1]:
                    decoy[0]['P'].detach_child(i.id)
            for i in ref[0]['P']:
                if i.id[1] < anchors[0] or i.id[1] > anchors[-1]:
                    ref[0]['P'].detach_child(i.id)

    # If you give the anchors AND flanking = True, the flanking L-RMSD will be calculated. Only if the peptide is
    # also longer than the binding core. The peptide binding core will be discarded

    elif anchors and flanking and decoy_p_len > 9:
        for x in range(ref_p_len):
            for i in decoy[0]['P']:
                if i.id[1] >= anchors[0] and i.id[1] <= anchors[-1]:
                    decoy[0]['P'].detach_child(i.id)
            for i in ref[0]['P']:
                if i.id[1] >= anchors[0] and i.id[1] <= anchors[-1]:
                    ref[0]['P'].detach_child(i.id)


    # remove c-like domain and keep only g domain
    decoy = remove_C_like_domain(decoy)
    ref = remove_C_like_domain(ref)

    # merge chains of the decoy
    decoy = merge_chains(decoy)
    # merge chains of the reference
    ref = merge_chains(ref)

    #Renumber and trim unaligned residues
    ############### Added
    ref, decoy , ref_sequences, decoy_sequences= renumber(ref, decoy)

    ref ,decoy = trim_indels(ref, decoy, ref_sequences, decoy_sequences)

    # trim unaligned atoms
    #ref, decoy , ref_sequences, decoy_sequences= renumber(ref, decoy)
    ref, decoy = remove_mismatched_atoms_from_pdb( ref, decoy,  atoms)


    # Write pdbs
    decoy_path = '%s/%s_decoy.pdb' % (output_dir, target_id)
    io = PDBIO()
    io.set_structure(decoy)
    io.save(decoy_path, select=NotDisordered())

    ref_path = '%s/%s_ref.pdb' % (output_dir, target_id)
    io = PDBIO()
    io.set_structure(ref)
    io.save(ref_path, select=NotDisordered())

    return decoy_path, ref_path

def get_Gdomain_lzone(ref_pdb, output_dir, MHC_class):
    """ Produce a lzone file for pdb2sql.
    Args:
        ref_pdb (str): path to the pdb file to use for the lzone
        output_dir (str): output directory
        MHC_class (str): Class of the MHC
    Raises:
        Exception: In case there are unexpected chain names it raises an exception
    Returns:
        outfile (str): Path to the output file
    """

    ref_name = ref_pdb.split('/')[-1].split('.')[0]
    outfile = '%s/%s.lzone' %(output_dir, ref_name)
    if MHC_class == 'I':
        with open(outfile, 'w') as output:
            P = PDBParser(QUIET=1)
            structure = P.get_structure('r', ref_pdb)
            for chain in structure.get_chains():
                if chain.id == 'M':
                    for x in range(2,173):
                        output.write('zone %s%i-%s%i\n' %(chain.id, x, chain.id, x))
                    #output.write('zone %s2-%s172\n' %(chain.id, chain.id))
                    #output.write('zone %s2-%s172:%s2-%s172\n' %(chain.id, chain.id, chain.id, chain.id))
                elif chain.id == 'P':
                    pass
                    #output.write('fit\n')
                    #for residue in chain:
                    #    if residue.id[2] == ' ':
                    #        output.write('rzone %s%s-%s%s\n' %(chain.id, str(residue.id[1]), chain.id, str(residue.id[1])))
                else:
                    raise Exception('Unrecognized chain ID, different from M or P. Please check your file')
            #output.write('fit\n')

    elif MHC_class == 'II':
        #Chain M from 4 to 72; Chain N from 10 to 80
        with open(outfile, 'w') as output:
            P = PDBParser(QUIET=1)
            structure = P.get_structure('r', ref_pdb)
            for chain in structure.get_chains():
                if chain.id == 'M':
                    output.write('zone %s4-%s72:%s4-%s72\n' %(chain.id, chain.id, chain.id, chain.id))
                elif chain.id == 'N':
                    output.write('zone %s10-%s80:%s10-%s80\n' %(chain.id, chain.id, chain.id, chain.id))
                elif chain.id == 'P':
                    pass
                    #output.write('fit\n')
                    #for residue in chain:
                    #    if residue.id[2] == ' ':
                    #        output.write('rzone %s%s-%s%s\n' %(chain.id, str(residue.id[1]), chain.id, str(residue.id[1])))
                else:
                    raise Exception('Unrecognized chain ID, different from M, N or P. Please check your file')
            #output.write('fit\n')
    return outfile

def remove_C_like_domain(pdb, need_to_be_removed=None):
    """ Removes the C-like domain from a MHC struture and keeps only the G domain
    Args:
        pdb: (Bio.PDB): Bio.PDB object with chains names M (N for MHCII) and P
        need_to_be_removed (list, optional):list of atoms to remove from M chain. Defaults to None.
    Returns: (Bio.PDB): Bio.PDB object without the C-like domain
    """

    # If MHCII, remove the C-like domain from the M-chain (res 80 and higher) and the N-chain (res 90 and higher)
    if 'N' in [chain.id for chain in pdb.get_chains()]:

        residue_ids_to_remove_N = [res.id for res in pdb[0]['N'] if res.id[1] > 90]
        #  Remove them
        for id in residue_ids_to_remove_N:
            pdb[0]['N'].detach_child(id)

        residue_ids_to_remove_M = [res.id for res in pdb[0]['M'] if res.id[1] > 80]
        #  Remove them
        for id in residue_ids_to_remove_M:
            pdb[0]['M'].detach_child(id)

    # If MHCI, remove the C-like domain, which is from residue 180+
    if 'N' not in [chain.id for chain in pdb.get_chains()]:
        for chain in pdb.get_chains():
            if chain.id == 'M':
                if need_to_be_removed ==None:
                    need_to_be_removed = [res.id for res in chain if res.id[1] > 180]
                    _ = [chain.detach_child(x) for x in need_to_be_removed]

                else:
                    #  Remove the list of given atom names 
                    for id in need_to_be_removed:
                        #if (chain.__contains__(id)):
                            chain.detach_child((' ', id, ' '))

                
    return pdb

#ValueError: Invalid column name lzone. Possible names are
#['rowID', 'serial', 'name', 'altLoc', 'resName', 'chainID', 'resSeq',
# 'iCode', 'x', 'y', 'z', 'occ', 'temp', 'element', 'model']


def trim_indels(pdb_ref, pdb_decoy, ref_sequences, decoy_sequences):
    ''' Trim indels for both reference and decoy PDBs

    Args:
        pdb_ref (Bio.PDB): object with chains names M (N for MHCII) and P
        pdb_decoy (Bio.PDB): object with chains names M (N for MHCII) and P
        ref_sequences : List of residue numbers in reference PDB 
        decoy_sequences : List of residue numbers in decoy PDB 
        

    Returns: Bio.PDB objects (reference and decoy) with matched residues
    '''

    ref_decoy_inconsistency=list(set(ref_sequences[0][1]) - set(decoy_sequences[0][1]))
    decoy_ref_inconsistency=list(set(decoy_sequences[0][1]) - set(ref_sequences[0][1]))
   
    need_to_be_removed_ref=list()

    ref_M_length=ref_sequences[0][-1][-1]

    for i  in  ref_decoy_inconsistency:
        
        need_to_be_removed_ref.append(i)
        if i>1 and ((i-1) in ref_sequences[0][1]):
            need_to_be_removed_ref.append(i-1)
        if i<ref_M_length and  ((i+1) in ref_sequences[0][1]):
            need_to_be_removed_ref.append(i+1)
 
    need_to_be_removed_ref=list(set(need_to_be_removed_ref))

    
    decoy_M_length = decoy_sequences[0][-1][-1]
    need_to_be_removed_decoy=list()
    for i  in  decoy_ref_inconsistency:
  
        need_to_be_removed_decoy.append(i)
        if i>1 and ((i-1) in decoy_sequences[0][1]): 
            need_to_be_removed_decoy.append(i-1)
        if i<decoy_M_length and ((i+1) in decoy_sequences[0][1]):
            need_to_be_removed_decoy.append(i+1)
            
    need_to_be_removed_decoy=list(set(need_to_be_removed_decoy))

    new_remove_decoy=list(set(need_to_be_removed_ref) & set(decoy_sequences[0][1]))
    need_to_be_removed_decoy = list(set(need_to_be_removed_decoy + new_remove_decoy))

    new_remove_ref=list(set(need_to_be_removed_decoy) & set(ref_sequences[0][1]))
    need_to_be_removed_ref = list(set(need_to_be_removed_ref + new_remove_ref))
    
    try:
        pdb_ref = remove_C_like_domain(pdb_ref, need_to_be_removed_ref)
    except:
        print('These residue IDs not in ref: %s' %need_to_be_removed_ref)
        print(ref_sequences[0][1])
    try:
        pdb_decoy = remove_C_like_domain(pdb_decoy, need_to_be_removed_decoy)
    except:
        print(decoy_sequences[0][1])
        print('The residue IDs not in decoy: %s' %need_to_be_removed_decoy)
    return pdb_ref, pdb_decoy


def get_residue_atoms(decoy_res, ref_res):
    ''' Retrieves all atom names of a given residue for both reference and decoy

    Args:
        ref_res (Bio.PDB.Residue): Residue object from reference PDB
        decoy_res (Bio.PDB.Residue): Residue object from decoy PDB         

    Returns: Bio.PDB.Atom objects (reference and decoy) of the given residues
    '''
    
    decoy_atoms= list()
    for atom in decoy_res.get_atoms():
        decoy_atoms.append(atom.id)

    ref_atoms = list()
    for atom in ref_res.get_atoms():
        ref_atoms.append(atom.id)

    return decoy_atoms, ref_atoms


def remove_atoms_from_res( ref_chain,  decoy_chain, i):
    '''
    Mismatched residues in Reference is removed both from decoy and reference
    '''
    ref_chain.detach_child((' ', i+1, ' '))
    decoy_chain.detach_child((' ', i+1, ' '))
    print("Mismatch residue in ref vs decoy in resi %s ->  removed in ref & decoy" %(i+1))   

def remove_mismatched_atoms_from_res(diff_atoms, dif_res, atoms):
    '''
    Mismatched atoms in the given residue object is removed
    '''
    diff_removed = list()

    for d in diff_atoms:
        if d in atoms:
            dif_res.detach_child(d)
            diff_removed.append(d)
    return diff_removed

def remove_mismatched_atoms_from_pdb(ref, decoy, atoms):
    """ Mismatched atoms in the given PDB object is removed
    Args:
        ref (Bio.PDB): object with chains names M (N for MHCII) and P
        decoy (Bio.PDB): object with chains names M (N for MHCII) and P
        atoms (_type_): _description_

    Returns:
        _type_: _description_
    """  
    for c in decoy.get_chains():

       chain_id = c.id
       print("chain %s:" %chain_id)
       decoy_m_len=len(decoy[0][chain_id])

       for i in range(decoy_m_len):
           try:
           
                ref_res = ref[0][chain_id][(' ', i+1, ' ')]
                decoy_res = decoy[0][chain_id][(' ', i+1, ' ')]

                # if ref_res.get_resname() != decoy_res.get_resname(): # Eliminate mismatched residues
                if ref_res.get_resname() == 'MSE': # Eliminate mismatched residues
                    remove_atoms_from_res( ref[0][chain_id],  decoy[0][chain_id], i)

                else: # residues matched

                   decoy_atoms, ref_atoms = get_residue_atoms(decoy_res, ref_res)  

                   if(set(decoy_atoms) != set(ref_atoms)) :
                       diff_decoy = set(decoy_atoms) - set(ref_atoms)
                       if diff_decoy:
                           decoy_diff_occured = remove_mismatched_atoms_from_res(diff_decoy, decoy_res, atoms)
       
                       diff_ref = set(ref_atoms) - set(decoy_atoms)
                       if diff_ref:
                           ref_diff_occured = remove_mismatched_atoms_from_res(diff_ref, ref_res, atoms)

                       if decoy_diff_occured and len(decoy_diff_occured) >0:
                           print("Mismatched atoms of Decoy vs ref in residue %d is: %s" %(i+1,decoy_diff_occured))
                       if ref_diff_occured and len(ref_diff_occured) >0:
                           print("Mismatched atoms of ref vs decoy in residue %d is: %s" %(i+1,ref_diff_occured))           
           except:
               continue
    return ref, decoy

class NotDisordered(Select):  # Inherit methods from Select class
    '''
    Keep one Alternative location for the given atom
    '''
    def accept_atom(self, atom):
        keepAltID = 'A'
        if (not atom.is_disordered()) or atom.get_altloc() == keepAltID:
            atom.set_altloc(" ")  # Eliminate alt location ID before output.
            return True
        else:  # Alt location was not one to be output.
            return False

1	from pyexpat import model	1✔
2	import os	1✔
3
4	from Bio import pairwise2	1✔
5	from Bio.SeqUtils import seq1	1✔
6	from Bio.PDB import PDBParser, PDBIO, Select	1✔
7
8	import PANDORA	1✔
9
10	class Model:	1✔
11
12	def __init__(self, target, model_path='', output_dir=False, pdb=False, molpdf=0, dope=0):	1✔
13	''' Initiate model object
14	Args:
15	target: Target object
16	model_path: (string) path to hypothetical model
17	output_dir: (string) output directory
18	pdb: Bio.PDB object of the hypothetical model
19	molpdf: (float) molpdf score
20	dope: (float) DOPE score
21	'''
22
23	self.target = target	1✔
24	self.model_path = model_path	1✔
25	self.molpdf = molpdf	1✔
26	self.dope = dope	1✔
27
28	# Define the output directory
29	if output_dir == False:	1✔
30	self.output_dir = f"{os.getcwd()}/{self.target.id}"	1✔
31	else:
32	self.output_dir = output_dir	1✔
33
34
35
36	# Check if the user gave either the path to the model pdb or the pdb itself.
37	if self.model_path == '' and not pdb:	1✔
38	raise Exception('Provide the path to a model structure or a Bio.PDB object')	×
39	# If there is a model path and no pdb, parse the pdb structure from that path.
40	if not pdb:	1✔
41	self.pdb = PDBParser(QUIET=True).get_structure(self.target.id, self.model_path)	1✔
42
43
44	def calc_LRMSD(self, reference_pdb, atoms = ['C', 'CA', 'N', 'O'], ligand_zone="whole"):	1✔
45	"""Calculate the L-RMSD between the decoy and reference structure (ground truth).
46	This function requires the pdb2sql module for L-RMSD calculation.
47
48	Args:
49	reference_pdb: Bio.PDB object or path to pdb file
50	atoms (list, optional): The list of atoms of the ligand selected to calculte the LRMSD . Defaults to ['C', 'CA', 'N', 'O'].
51	ligand_zone (str, optional): The region of the Ligand selected to calculate the LRMSD. Defaults to "whole".
52
53	Raises:
54	Exception: PDB2SQL LRMSD claulation failed
55	Returns: (float) L-RMSD calculated by PDB2SQL
56	"""
57
58	#from pdb2sql import pdb2sql, superpose, StructureSimilarity
59	from pdb2sql import StructureSimilarity	1✔
60
61	# load target pdb
62	if isinstance(reference_pdb, str): # if its a string, it should be the path of the pdb, then load pdb first	1✔
63	ref = PDBParser(QUIET=True).get_structure(self.target.id, reference_pdb)	1✔
64	else:
65	ref = reference_pdb	×
66
67	# Define file names as variables
68	model_name = self.model_path.split('/')[-1].split('.')[1]	1✔
69	#decoy_path = '%s/%s_decoy.pdb' % (self.output_dir, model_name)
70	#ref_path = '%s/%s_ref.pdb' % (self.output_dir, model_name)
71
72	# Define zones to align
73	#M_lzone = list(range(4,73))
74	#N_lzone = list(range(10,80))
75
76	## pdb2sql needs 1 big chain and 1 ligand chain with correct numbering, for MHCII, this means merging the chains.
77	#if os.path.exists(decoy_path)==False and ligand_zone.lowe()=="whole":
78	if ligand_zone.lower()=="whole":	1✔
79	decoy_path, ref_path = homogenize_pdbs(self.pdb, ref, atoms, self.output_dir, model_name)	1✔
80	## calculates core LRMSD
81	elif ligand_zone.lower() == "core":	1✔
82	decoy_path, ref_path = homogenize_pdbs(self.pdb, ref, atoms, self.output_dir, model_name, anchors = self.target.anchors)	1✔
83	elif ligand_zone.lower() == "flanking":	×
84	decoy_path, ref_path = homogenize_pdbs(self.pdb, ref, atoms, self.output_dir,model_name, anchors = self.target.anchors, flanking=True)	×
85
86
87	start_dir = os.getcwd()	1✔
88	os.chdir(self.output_dir)	1✔
89	# Produce lzone file for the l-rmsd calculation
90	#lzone = get_Gdomain_lzone(ref_path, self.output_dir, self.target.MHC_class)
91	#TODO: check if it's MHC I or II and adapt for chain M and N
92	# Get decoy structure to superpose
93	#decoy_db = pdb2sql(decoy_path)
94	#decoy_lzone = np.asarray(decoy_db.get('x,y,z', resSeq=M_lzone))
95
96	# Get ref structure to superpose
97	#ref_db = pdb2sql(ref_path)
98	#ref_lzone = np.asarray(ref_db.get('x,y,z', resSeq=M_lzone))
99
100	# Align the G domains
101	#superpose.superpose_selection()
102
103	try:	1✔
104	# Calculate l-rmsd between decoy and reference with pdb2sql
105	sim = StructureSimilarity(decoy_path, ref_path)	1✔
106	#self.lrmsd = sim.compute_lrmsd_fast(method='svd', name=atoms, lzone = lzone)
107	if ligand_zone.lower()=="whole":	1✔
108	self.lrmsd = sim.compute_lrmsd_pdb2sql(exportpath=None, method='svd', name = atoms)	1✔
109	## calculates LRMSD for peptide core region
110	elif ligand_zone.lower() == "core":	1✔
111	self.core_lrmsd = sim.compute_lrmsd_pdb2sql(exportpath=None, method='svd', name=atoms)	1✔
112	## calculates RMSD for peptide flanking region
113	elif ligand_zone.lower() == "flanking":	×
114	self.flanking_lrmsd = sim.compute_lrmsd_pdb2sql(exportpath=None, method='svd', name=atoms)	×
115	except:	×
116	print('An error occurred while calculating the rmsd for target %s, model %s for %s' %(ref_path, decoy_path, ligand_zone))	×
UNCOV 117	raise Exception('Please check your model and ref info for model %s' %self.model_path)	×
118
119	# remove intermediate files
120	os.system('rm %s %s' %(decoy_path, ref_path))	1✔
121	os.chdir(start_dir)	1✔
122
123	def merge_chains(pdb):	1✔
124	''' Merges two chains of MHCII to one chain. pdb2sql can only calculate L-rmsd with one chain.
125	Args:
126	pdb: Bio.PDB object
127	Returns: Bio.PDB object with its M and N chain merged as M chain
128	'''
129	# Merge chains
130	if 'N' in [chain.id for chain in pdb.get_chains()]:	1✔
131
132	for j in pdb[0]['N'].get_residues():	×
133	j.id = (j.id[0], j.id[1], 'M')	×
134	pdb[0]['M'].add(j)	×
135
136	for i in pdb.get_chains():	×
137	for model in pdb:	×
138	for chain in model:	×
139	if chain.id in ['N']:	×
140	model.detach_child(chain.id)	×
141	return pdb	1✔
142
143
144	def renumber(pdb_ref, pdb_decoy, custom_map={"MSE":"M"}):	1✔
145	''' aligns two pdb's and renumber them accordingly.
146	Args:
147	pdb_ref: Bio.PDB object
148	pdb_decoy: Bio.PDB object
149
150	Returns: Bio.PDB objects with renumbered residues
151	'''
152	ref_sequences = [[chain.id, seq1(''.join([res.resname for res in chain]), custom_map=custom_map)] for chain in pdb_ref.get_chains()]	1✔
153	#[[chain.id, ('').join([seq1(res.resname) for res in chain])]
154	# for chain in pdb_ref.get_chains()]
155	ref_sequences.sort()	1✔
156	decoy_sequences = [[chain.id, seq1(''.join([res.resname for res in chain]), custom_map=custom_map)] for chain in pdb_decoy.get_chains()]	1✔
157	#[[chain.id, ('').join([seq1(res.resname) for res in chain])]
158	#for chain in pdb_decoy.get_chains()]
159
160	decoy_sequences.sort()	1✔
161
162	assert(len(ref_sequences) == len(decoy_sequences))	1✔
163
164
165	for ind in range(len(ref_sequences)):	1✔
166	pair = pairwise2.align.globalxx(ref_sequences[ind][1], decoy_sequences[ind][1])[0]	1✔
167	ref_sequences[ind][1] = pair.seqA	1✔
168	decoy_sequences[ind][1] = pair.seqB	1✔
169
170	ref_sequences = [[seq[0],[i+1 for i,res in enumerate(seq[1]) if res != '-']] for seq in ref_sequences]	1✔
171	decoy_sequences = [[seq[0],[i+1 for i,res in enumerate(seq[1]) if res != '-']] for seq in decoy_sequences]	1✔
172
173
174
175	def assign(pdb, pdb_sequences):	1✔
176	''' Renumbers the pdb using aligned sequences.
177	Args:
178	pdb_ref: Bio.PDB object
179	pdb_decoy: Bio.PDB object
180
181	Returns: Bio.PDB objects with renumbered residues
182	'''
183
184
185	for chain in pdb.get_chains():	1✔
186	for seq in pdb_sequences:	1✔
187	if chain.id == seq[0]:	1✔
188	for ind, res in enumerate(chain):	1✔
189	res.id = ('X', seq[1][ind], res.id[2])	1✔
190	for chain in pdb.get_chains():	1✔
191	for res in chain:	1✔
192	res.id = (' ', res.id[1], ' ')	1✔
193	return pdb	1✔
194
195	pdb_ref = assign(pdb_ref, ref_sequences)	1✔
196	pdb_decoy = assign(pdb_decoy, decoy_sequences)	1✔
197
198	return pdb_ref, pdb_decoy , ref_sequences, decoy_sequences	1✔
199
200	def homogenize_pdbs(decoy, ref, atoms, output_dir, target_id = 'MHC' , anchors =False, flanking=False):	1✔
201	''' Make sure that the decoy and reference structure have the same structure sequences.
202	Args:
203	decoy: Bio.PDB object of the decoy structure
204	ref: Bio.PDB object of the reference structure
205	output_dir: (string) directory that is used to write intermediate files
206	Returns: (tuple) Bio.PDB objects with the same structure sequence
207	'''
208	ref_p_len=len(ref[0]['P'])	1✔
209	decoy_p_len=len(decoy[0]['P'])	1✔
210
211
212	# If you give the anchors, the core L-RMSD will be calculated.
213	# The peptide residues before and after the first and last anchor residue will be discarded.
214	if anchors and not flanking:	1✔
215	for x in range(len(decoy[0]['P'])):	1✔
216	for i in decoy[0]['P']:	1✔
217	if i.id[1] < anchors[0] or i.id[1] > anchors[-1]:	1✔
218	decoy[0]['P'].detach_child(i.id)	1✔
219	for i in ref[0]['P']:	1✔
220	if i.id[1] < anchors[0] or i.id[1] > anchors[-1]:	1✔
221	ref[0]['P'].detach_child(i.id)	1✔
222
223	# If you give the anchors AND flanking = True, the flanking L-RMSD will be calculated. Only if the peptide is
224	# also longer than the binding core. The peptide binding core will be discarded
225
226	elif anchors and flanking and decoy_p_len > 9:	1✔
227	for x in range(ref_p_len):	×
228	for i in decoy[0]['P']:	×
229	if i.id[1] >= anchors[0] and i.id[1] <= anchors[-1]:	×
230	decoy[0]['P'].detach_child(i.id)	×
231	for i in ref[0]['P']:	×
232	if i.id[1] >= anchors[0] and i.id[1] <= anchors[-1]:	×
233	ref[0]['P'].detach_child(i.id)	×
234
235
236	# remove c-like domain and keep only g domain
237	decoy = remove_C_like_domain(decoy)	1✔
238	ref = remove_C_like_domain(ref)	1✔
239
240	# merge chains of the decoy
241	decoy = merge_chains(decoy)	1✔
242	# merge chains of the reference
243	ref = merge_chains(ref)	1✔
244
245	#Renumber and trim unaligned residues
246	############### Added
247	ref, decoy , ref_sequences, decoy_sequences= renumber(ref, decoy)	1✔
248
249	ref ,decoy = trim_indels(ref, decoy, ref_sequences, decoy_sequences)	1✔
250
251	# trim unaligned atoms
252	#ref, decoy , ref_sequences, decoy_sequences= renumber(ref, decoy)
253	ref, decoy = remove_mismatched_atoms_from_pdb( ref, decoy, atoms)	1✔
254
255
256	# Write pdbs
257	decoy_path = '%s/%s_decoy.pdb' % (output_dir, target_id)	1✔
258	io = PDBIO()	1✔
259	io.set_structure(decoy)	1✔
260	io.save(decoy_path, select=NotDisordered())	1✔
261
262	ref_path = '%s/%s_ref.pdb' % (output_dir, target_id)	1✔
263	io = PDBIO()	1✔
264	io.set_structure(ref)	1✔
265	io.save(ref_path, select=NotDisordered())	1✔
266
267	return decoy_path, ref_path	1✔
268
269	def get_Gdomain_lzone(ref_pdb, output_dir, MHC_class):	1✔
270	""" Produce a lzone file for pdb2sql.
271	Args:
272	ref_pdb (str): path to the pdb file to use for the lzone
273	output_dir (str): output directory
274	MHC_class (str): Class of the MHC
275	Raises:
276	Exception: In case there are unexpected chain names it raises an exception
277	Returns:
278	outfile (str): Path to the output file
279	"""
280
281	ref_name = ref_pdb.split('/')[-1].split('.')[0]	×
282	outfile = '%s/%s.lzone' %(output_dir, ref_name)	×
283	if MHC_class == 'I':	×
284	with open(outfile, 'w') as output:	×
285	P = PDBParser(QUIET=1)	×
286	structure = P.get_structure('r', ref_pdb)	×
287	for chain in structure.get_chains():	×
288	if chain.id == 'M':	×
289	for x in range(2,173):	×
290	output.write('zone %s%i-%s%i\n' %(chain.id, x, chain.id, x))	×
291	#output.write('zone %s2-%s172\n' %(chain.id, chain.id))
292	#output.write('zone %s2-%s172:%s2-%s172\n' %(chain.id, chain.id, chain.id, chain.id))
293	elif chain.id == 'P':	×
294	pass	×
295	#output.write('fit\n')
296	#for residue in chain:
297	# if residue.id[2] == ' ':
298	# output.write('rzone %s%s-%s%s\n' %(chain.id, str(residue.id[1]), chain.id, str(residue.id[1])))
299	else:
300	raise Exception('Unrecognized chain ID, different from M or P. Please check your file')	×
301	#output.write('fit\n')
302
303	elif MHC_class == 'II':	×
304	#Chain M from 4 to 72; Chain N from 10 to 80
305	with open(outfile, 'w') as output:	×
306	P = PDBParser(QUIET=1)	×
307	structure = P.get_structure('r', ref_pdb)	×
308	for chain in structure.get_chains():	×
309	if chain.id == 'M':	×
310	output.write('zone %s4-%s72:%s4-%s72\n' %(chain.id, chain.id, chain.id, chain.id))	×
311	elif chain.id == 'N':	×
312	output.write('zone %s10-%s80:%s10-%s80\n' %(chain.id, chain.id, chain.id, chain.id))	×
313	elif chain.id == 'P':	×
314	pass	×
315	#output.write('fit\n')
316	#for residue in chain:
317	# if residue.id[2] == ' ':
318	# output.write('rzone %s%s-%s%s\n' %(chain.id, str(residue.id[1]), chain.id, str(residue.id[1])))
319	else:
320	raise Exception('Unrecognized chain ID, different from M, N or P. Please check your file')	×
321	#output.write('fit\n')
322	return outfile	×
323
324	def remove_C_like_domain(pdb, need_to_be_removed=None):	1✔
325	""" Removes the C-like domain from a MHC struture and keeps only the G domain
326	Args:
327	pdb: (Bio.PDB): Bio.PDB object with chains names M (N for MHCII) and P
328	need_to_be_removed (list, optional):list of atoms to remove from M chain. Defaults to None.
329	Returns: (Bio.PDB): Bio.PDB object without the C-like domain
330	"""
331
332	# If MHCII, remove the C-like domain from the M-chain (res 80 and higher) and the N-chain (res 90 and higher)
333	if 'N' in [chain.id for chain in pdb.get_chains()]:	1✔
334
335	residue_ids_to_remove_N = [res.id for res in pdb[0]['N'] if res.id[1] > 90]	×
336	# Remove them
337	for id in residue_ids_to_remove_N:	×
338	pdb[0]['N'].detach_child(id)	×
339
340	residue_ids_to_remove_M = [res.id for res in pdb[0]['M'] if res.id[1] > 80]	×
341	# Remove them
342	for id in residue_ids_to_remove_M:	×
343	pdb[0]['M'].detach_child(id)	×
344
345	# If MHCI, remove the C-like domain, which is from residue 180+
346	if 'N' not in [chain.id for chain in pdb.get_chains()]:	1✔
347	for chain in pdb.get_chains():	1✔
348	if chain.id == 'M':	1✔
349	if need_to_be_removed ==None:	1✔
350	need_to_be_removed = [res.id for res in chain if res.id[1] > 180]	1✔
351	_ = [chain.detach_child(x) for x in need_to_be_removed]	1✔
352
353	else:
354	# Remove the list of given atom names
355	for id in need_to_be_removed:	1✔
356	#if (chain.__contains__(id)):
357	chain.detach_child((' ', id, ' '))	×
358
359
360	return pdb	1✔
361
362	#ValueError: Invalid column name lzone. Possible names are
363	#['rowID', 'serial', 'name', 'altLoc', 'resName', 'chainID', 'resSeq',
364	# 'iCode', 'x', 'y', 'z', 'occ', 'temp', 'element', 'model']
365
366
367	def trim_indels(pdb_ref, pdb_decoy, ref_sequences, decoy_sequences):	1✔
368	''' Trim indels for both reference and decoy PDBs
369
370	Args:
371	pdb_ref (Bio.PDB): object with chains names M (N for MHCII) and P
372	pdb_decoy (Bio.PDB): object with chains names M (N for MHCII) and P
373	ref_sequences : List of residue numbers in reference PDB
374	decoy_sequences : List of residue numbers in decoy PDB
375
376
377	Returns: Bio.PDB objects (reference and decoy) with matched residues
378	'''
379
380	ref_decoy_inconsistency=list(set(ref_sequences[0][1]) - set(decoy_sequences[0][1]))	1✔
381	decoy_ref_inconsistency=list(set(decoy_sequences[0][1]) - set(ref_sequences[0][1]))	1✔
382
383	need_to_be_removed_ref=list()	1✔
384
385	ref_M_length=ref_sequences[0][-1][-1]	1✔
386
387	for i in ref_decoy_inconsistency:	1✔
388
389	need_to_be_removed_ref.append(i)	×
390	if i>1 and ((i-1) in ref_sequences[0][1]):	×
391	need_to_be_removed_ref.append(i-1)	×
392	if i<ref_M_length and ((i+1) in ref_sequences[0][1]):	×
393	need_to_be_removed_ref.append(i+1)	×
394
395	need_to_be_removed_ref=list(set(need_to_be_removed_ref))	1✔
396
397
398	decoy_M_length = decoy_sequences[0][-1][-1]	1✔
399	need_to_be_removed_decoy=list()	1✔
400	for i in decoy_ref_inconsistency:	1✔
401
402	need_to_be_removed_decoy.append(i)	×
403	if i>1 and ((i-1) in decoy_sequences[0][1]):	×
404	need_to_be_removed_decoy.append(i-1)	×
405	if i<decoy_M_length and ((i+1) in decoy_sequences[0][1]):	×
406	need_to_be_removed_decoy.append(i+1)	×
407
408	need_to_be_removed_decoy=list(set(need_to_be_removed_decoy))	1✔
409
410	new_remove_decoy=list(set(need_to_be_removed_ref) & set(decoy_sequences[0][1]))	1✔
411	need_to_be_removed_decoy = list(set(need_to_be_removed_decoy + new_remove_decoy))	1✔
412
413	new_remove_ref=list(set(need_to_be_removed_decoy) & set(ref_sequences[0][1]))	1✔
414	need_to_be_removed_ref = list(set(need_to_be_removed_ref + new_remove_ref))	1✔
415
416	try:	1✔
417	pdb_ref = remove_C_like_domain(pdb_ref, need_to_be_removed_ref)	1✔
418	except:	×
419	print('These residue IDs not in ref: %s' %need_to_be_removed_ref)	×
420	print(ref_sequences[0][1])	×
421	try:	1✔
422	pdb_decoy = remove_C_like_domain(pdb_decoy, need_to_be_removed_decoy)	1✔
423	except:	×
424	print(decoy_sequences[0][1])	×
425	print('The residue IDs not in decoy: %s' %need_to_be_removed_decoy)	×
426	return pdb_ref, pdb_decoy	1✔
427
428
429	def get_residue_atoms(decoy_res, ref_res):	1✔
430	''' Retrieves all atom names of a given residue for both reference and decoy
431
432	Args:
433	ref_res (Bio.PDB.Residue): Residue object from reference PDB
434	decoy_res (Bio.PDB.Residue): Residue object from decoy PDB
435
436	Returns: Bio.PDB.Atom objects (reference and decoy) of the given residues
437	'''
438
439	decoy_atoms= list()	1✔
440	for atom in decoy_res.get_atoms():	1✔
441	decoy_atoms.append(atom.id)	1✔
442
443	ref_atoms = list()	1✔
444	for atom in ref_res.get_atoms():	1✔
445	ref_atoms.append(atom.id)	1✔
446
447	return decoy_atoms, ref_atoms	1✔
448
449
450	def remove_atoms_from_res( ref_chain, decoy_chain, i):	1✔
451	'''
452	Mismatched residues in Reference is removed both from decoy and reference
453	'''
454	ref_chain.detach_child((' ', i+1, ' '))	×
455	decoy_chain.detach_child((' ', i+1, ' '))	×
456	print("Mismatch residue in ref vs decoy in resi %s -> removed in ref & decoy" %(i+1))	×
457
458	def remove_mismatched_atoms_from_res(diff_atoms, dif_res, atoms):	1✔
459	'''
460	Mismatched atoms in the given residue object is removed
461	'''
462	diff_removed = list()	×
463
464	for d in diff_atoms:	×
465	if d in atoms:	×
466	dif_res.detach_child(d)	×
467	diff_removed.append(d)	×
468	return diff_removed	×
469
470	def remove_mismatched_atoms_from_pdb(ref, decoy, atoms):	1✔
471	""" Mismatched atoms in the given PDB object is removed
472	Args:
473	ref (Bio.PDB): object with chains names M (N for MHCII) and P
474	decoy (Bio.PDB): object with chains names M (N for MHCII) and P
475	atoms (_type_): _description_
476
477	Returns:
478	_type_: _description_
479	"""
480	for c in decoy.get_chains():	1✔
481
482	chain_id = c.id	1✔
483	print("chain %s:" %chain_id)	1✔
484	decoy_m_len=len(decoy[0][chain_id])	1✔
485
486	for i in range(decoy_m_len):	1✔
487	try:	1✔
488
489	ref_res = ref[0][chain_id][(' ', i+1, ' ')]	1✔
490	decoy_res = decoy[0][chain_id][(' ', i+1, ' ')]	1✔
491
492	# if ref_res.get_resname() != decoy_res.get_resname(): # Eliminate mismatched residues
493	if ref_res.get_resname() == 'MSE': # Eliminate mismatched residues	1✔
494	remove_atoms_from_res( ref[0][chain_id], decoy[0][chain_id], i)	×
495
496	else: # residues matched
497
498	decoy_atoms, ref_atoms = get_residue_atoms(decoy_res, ref_res)	1✔
499
500	if(set(decoy_atoms) != set(ref_atoms)) :	1✔
501	diff_decoy = set(decoy_atoms) - set(ref_atoms)	×
502	if diff_decoy:	×
503	decoy_diff_occured = remove_mismatched_atoms_from_res(diff_decoy, decoy_res, atoms)	×
504
505	diff_ref = set(ref_atoms) - set(decoy_atoms)	×
506	if diff_ref:	×
507	ref_diff_occured = remove_mismatched_atoms_from_res(diff_ref, ref_res, atoms)	×
508
509	if decoy_diff_occured and len(decoy_diff_occured) >0:	×
510	print("Mismatched atoms of Decoy vs ref in residue %d is: %s" %(i+1,decoy_diff_occured))	×
511	if ref_diff_occured and len(ref_diff_occured) >0:	×
512	print("Mismatched atoms of ref vs decoy in residue %d is: %s" %(i+1,ref_diff_occured))	×
513	except:	×
514	continue	×
515	return ref, decoy	1✔
516
517	class NotDisordered(Select): # Inherit methods from Select class	1✔
518	'''
519	Keep one Alternative location for the given atom
520	'''
521	def accept_atom(self, atom):	1✔
522	keepAltID = 'A'	1✔
523	if (not atom.is_disordered()) or atom.get_altloc() == keepAltID:	1✔
524	atom.set_altloc(" ") # Eliminate alt location ID before output.	1✔
525	return True	1✔
526	else: # Alt location was not one to be output.
NEW 527	return False	×

X-lab-3D / PANDORA / 14772793414

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous