Python PDB.NeighborSearch类代码示例

 def extract_feature(self):
     seed(self.seed)
     counter = 0
     print_info_nn(" >>> Adding D2 category based shape distribution for database {0} ... ".format(self._database.name))
     overall_time = datetime.now()
     if not os.path.exists(self._get_dir_name()):
         os.makedirs(self._get_dir_name())
     for complex_name in self._database.complexes.keys():
         protein_complex = self._database.complexes[complex_name]
         proteins = [protein_complex.unbound_formation.ligand, protein_complex.unbound_formation.receptor]
         for protein in proteins:
             shape_dist_file = self._get_dir_name() + protein.name
             if not os.path.exists(shape_dist_file + ".npy"):
                 counter += 1
                 if counter <= 15:
                     print_info_nn("{0}, ".format(protein.name))
                 else:
                     counter = 0
                     print_info("{0}".format(protein.name))
                 atoms = protein.atoms
                 neighbour_search = NeighborSearch(atoms)
                 distributions = np.zeros((len(protein.residues), self.number_of_bins))
                 for i in range(len(protein.residues)):
                     residue = protein.residues[i]
                     nearby_residues = neighbour_search.search(residue.center, self.radius, "R")
                     distributions[i, :] = self._compute_distribution(nearby_residues)
                 np.save(shape_dist_file, distributions)
             distributions = np.load(shape_dist_file + ".npy")
             for i in range(len(protein.residues)):
                 protein.residues[i].add_feature(Features.D2_CATEGORY_SHAPE_DISTRIBUTION, distributions[i, :])
     print_info("took {0} seconds.".format((datetime.now() - overall_time).seconds))

def get_residues_distance_distribution(data, r):
    """
    This function computes the distribution of number of residues in certain radius around residues.

    """
    files = glob.glob(data + "*_b.pdb")
    dist = {}
    files.sort()
    file_counter = 0
    for bound_pbd in files:
        dist[bound_pbd] = []
    l = []
    for bound_pbd in files:
        file_counter += 1
        print "Protein " + str(file_counter) + "/" + str(len(files))
        s, a, r = read_pdb_file(bound_pbd)
        ns = NeighborSearch(a)
        res_counter = 0
        for res in r:
            b = 0 * res.child_list[0].get_coord()
            for atom in res.child_list:
                b += atom.get_coord()
            center = b / len(res.child_list)
            l.append(len(ns.search(center, 100, "R")))
            res_counter += 1
            # print "Residue " + str(res_counter) + "out of " + str(len(r))
    plot(np.bincount(l))
    show()
    print files

def check_clash(str_name, v=True):
        """check_clash, fract of clashes!

        if zero contacts then error -> fix ->

        Problem, contacts, str_name: 311 505 na-prot_13536.pdb
        Sterical clashes  0.615841584158

        c is counter
        """
        print(str_name)
        structure = open(str_name)
        #model = structure[0]
        atoms_A = []
        atoms_B = []
        for line in structure.readlines():
            if line[:4] == "ATOM":
                #print line
                at_nam = line[12:16].strip()
                coor = [float(line[30:38]),float(line[38:46]), float(line[46:54])]	
                at = Atom.Atom(at_nam,coor,0.0,1.0,' ',at_nam,1,at_nam[0])
                if line[21] == "A":
                    atoms_A.append(at)
                elif line[21] == "B":
                    atoms_B.append(at)
                else: pass
        #atoms_B = Selection.unfold_entities(structure[0]['B'], 'A')
        #print len(atoms_A), len(atoms_B)
        if len(atoms_A) > len(atoms_B):
            less = atoms_B
            more = atoms_A
        else: 
            less = atoms_A
            more = atoms_B
        problem = 0
        contacts = 0 
        ns=NeighborSearch(more)
        for at in less:
             neighbors=ns.search(array(at.get_coord()),2.0,'A')
             if neighbors != []:
                 problem +=1
                 contacts +=1
             else:
                 neighbors1=ns.search(array(at.get_coord()),4.0,'A')
                 if neighbors1 != []:
                     contacts +=1
        if v:
                print('problem:', float(problem))
                print('contacts:', float(contacts))
        try:
            fract = float(problem)/float(contacts)
        except ZeroDivisionError:
            fract = problem # or skip this structure
            print('ZeroDivison -- skip:', problem, contacts, str_name)
            return fract

        #print 'Contacts, str_name:', problem, contacts, str_name, "Sterical clashes ", fract
        return fract

def get_interactions_between_chains(model, chain_id_1, chain_id_2, r_cutoff=6):
    """Calculate interactions between the residues of the two chains.

    An interaction is defines as a pair of residues where at least one pair of atom
    is closer than r_cutoff.

    .. deprecated:: 1.0
        Use python:fn:`get_interacting_residues` instead.
        It gives you both the residue index and the resnum.

    Returns
    -------
    OrderedDict
        Keys are (residue_number, residue_amino_acid) tuples
        (e.g. ('0', 'M'), ('1', 'Q'), ...).
        Values are lists of (residue_number, residue_amino_acid) tuples.
        (e.g. [('0', 'M'), ('1', 'Q'), ...]).
    """
    try:
        from Bio.PDB import NeighborSearch
    except ImportError as e:
        logger.warning('Importing Biopython NeighborSearch returned an error: {}'.format(e))
        logger.warning('Using the the slow version of the neighbour-finding algorithm...')
        return get_interactions_between_chains_slow(model, chain_id_1, chain_id_2, r_cutoff)

    # Extract the chains of interest from the model
    chain_1 = None
    chain_2 = None
    for child in model.get_list():
        if child.id == chain_id_1:
            chain_1 = child
        if child.id == chain_id_2:
            chain_2 = child
    if chain_1 is None or chain_2 is None:
        raise Exception('Chains %s and %s were not found in the model' % (chain_id_1, chain_id_2))

    ns = NeighborSearch(list(chain_2.get_atoms()))
    interactions_between_chains = OrderedDict()
    for idx, residue_1 in enumerate(chain_1):
        if residue_1.resname in AMINO_ACIDS and residue_1.id[0] == ' ':
            resnum_1 = str(residue_1.id[1]) + residue_1.id[2].strip()
            resaa_1 = convert_aa(residue_1.get_resname(), quiet=True)
            interacting_residues = set()
            for atom_1 in residue_1:
                interacting_residues.update(ns.search(atom_1.get_coord(), r_cutoff, 'R'))
            interacting_resids = []
            for residue_2 in interacting_residues:
                resnum_2 = str(residue_2.id[1]) + residue_2.id[2].strip()
                resaa_2 = convert_aa(residue_2.get_resname(), quiet=True)
                if residue_2.resname in AMINO_ACIDS and residue_2.id[0] == ' ':
                    interacting_resids.append((resnum_2, resaa_2,))
            if interacting_resids:
                interacting_resids.sort(
                    key=lambda x: int(''.join([c for c in x[0] if c.isdigit()])))
                interactions_between_chains[(resnum_1, resaa_1)] = interacting_resids
    return interactions_between_chains

 def interaction(self, pdb_id, filename, domain_1, domain_2):
     """Returns a dict with informations (atoms, residues...) if two domains
     interact with each other, and returns False if not."""
     print "Searching for interactions in "+pdb_id+"..."
     # creates a strucuture object/class to extract atoms of the two domains
     model = structure(pdb_id).get_model(pdb_id, filename)
     residues_1 = structure(pdb_id).get_residues(model, domain_1)
     residues_2 = structure(pdb_id).get_residues(model, domain_2)
     atoms_1 = Selection.unfold_entities(residues_1, 'A')
     atoms_2 = Selection.unfold_entities(residues_2, 'A')
     # gets the serial numbers of the atoms
     numbers_1 = structure(pdb_id).serial_numbers(atoms_1)
     numbers_2 = structure(pdb_id).serial_numbers(atoms_2)
     # the search starts here !
     atoms = Selection.unfold_entities(model, 'A')
     nsearch = NeighborSearch(atoms)
     interacting_atoms_1 = []
     interacting_atoms_2 = []
     for atom in atoms:
         if atom.get_serial_number() in numbers_1:
             point = atom.get_coord()
             # This is how we detect an interaction, we put 5 angstroms
             # here.
             # This is the simplest method we can use, and we're not sure
             # that it is correct.
             # Originally we have planned to go further by doing a surface
             # and accesssion analysis, but we had no time.
             # We hope we can talk about that during the talk.
             neighbors = nsearch.search(point, 5)
             for neighbor in neighbors:
                 if neighbor.get_serial_number() in numbers_2:
                     interacting_atoms_2.append(neighbor)
                     if atom not in interacting_atoms_1:
                         interacting_atoms_1.append(atom)
     # returns a dict with all residues and atoms
     if len(interacting_atoms_2) > 0:
         infos = {}
         infos['1'] = {}
         infos['2'] = {}
         # just get the parent residues for the list of atoms
         interacting_residues_1 = structure(pdb_id).atoms2residues(
                 interacting_atoms_1)
         interacting_residues_2 = structure(pdb_id).atoms2residues(
                 interacting_atoms_2)
         infos['1']['atoms'] = interacting_atoms_1
         infos['2']['atoms'] = interacting_atoms_2
         infos['1']['residues'] = interacting_residues_1
         infos['2']['residues'] = interacting_residues_2
         return infos
     else: return False

def calculate_ic(structure, d_cutoff=5.5, selection=None):
    """
    Calculates intermolecular contacts in a parsed structure object.
    """
    atom_list = list(structure.get_atoms())
    ns = NeighborSearch(atom_list)
    all_list = ns.search_all(radius=d_cutoff, level='R')

    if selection:
        _sd = selection_dict
        ic_list = [c for c in all_list if (c[0].parent.id in _sd and c[1].parent.id in _sd)
                    and (_sd[c[0].parent.id] != _sd[c[1].parent.id]) ]
    else:
        ic_list = [c for c in all_list if c[0].parent.id != c[1].parent.id]

    if not ic_list:
        raise ValueError('No contacts found for selection')

    return ic_list

 def _remove_distant_hatatms(self, new_model, hetatm_chain):
     """Detach hetatms that are more than ``self.r_cutoff`` away from the main chain(s)."""
     ns = NeighborSearch(list(new_model.get_atoms()))
     hetatm_chain.id = [
         c for c in reversed(string.ascii_uppercase) if c not in self.chain_ids][0]
     res_idx = 0
     while res_idx < len(hetatm_chain):
         res_1 = hetatm_chain.child_list[res_idx]
         in_contact = False
         for atom_1 in res_1:
             interacting_residues = ns.search(atom_1.get_coord(), self.r_cutoff, 'R')
             if interacting_residues:
                 # logger.debug(res_1.id)
                 # logger.debug(interacting_residues)
                 in_contact = True
         if in_contact:
             res_idx += 1
             continue
         # logger.debug('Detaching child: {}'.format(res_1.id))
         hetatm_chain.detach_child(res_1.id)

    def extract_feature(self):
        seed(self.seed)
        print_info_nn(" >>> Adding D1 surface shape distribution for database {0} ... ".format(self._database.name))
        overall_time = datetime.now()
        counter = 0
        if not os.path.exists(self._get_dir_name()):
            os.makedirs(self._get_dir_name())
        for complex_name in self._database.complexes.keys():
            protein_complex = self._database.complexes[complex_name]
            proteins = [protein_complex.unbound_formation.ligand, protein_complex.unbound_formation.receptor]
            for protein in proteins:
                shape_dist_file = self._get_dir_name() + protein.name
                if not os.path.exists(shape_dist_file + ".npy"):
                    counter += 1
                    if counter <= 15:
                        print_info_nn("{0}, ".format(protein.name))
                    else:
                        counter = 0
                        print_info("{0}".format(protein.name))
                    atoms = protein.atoms
                    neighbour_search = NeighborSearch(atoms)
                    distributions = np.zeros((len(protein.residues), self.number_of_bins + 1))
                    for i in range(len(protein.residues)):
                        residue = protein.residues[i]
                        nearby_residues = [protein.biopython_residues[i]]
                        temp_nearby_residues = neighbour_search.search(residue.center, self.radius, "R")
                        for nearby_residue in temp_nearby_residues:
                            if nearby_residue not in protein.biopython_residues:
                                continue
                            residues_index = protein.biopython_residues.index(nearby_residue)
                            residue = protein.residues[residues_index]

                            if residue.get_feature(Features.RELATIVE_ACCESSIBLE_SURFACE_AREA) >= self.rASA_threshold:
                                nearby_residues.append(nearby_residue)
                        distributions[i, :] = self._compute_distribution(nearby_residues, residue.center)
                    np.save(shape_dist_file, distributions)
                distributions = np.load(shape_dist_file + ".npy")
                for i in range(len(protein.residues)):
                    protein.residues[i].add_feature(Features.D1_SURFACE_SHAPE_DISTRIBUTION, distributions[i, :])
        print_info("took {0} seconds.".format((datetime.now() - overall_time).seconds))

structure = PDBParser().get_structure('X', args.pdb)

center_atoms = []
pymol_command = ""

all_atom_list = [atom for atom in structure.get_atoms() if atom.name == 'CA' ]

for k in args.chain1 :
    chain_atoms = [atom for atom in structure[0][k].get_atoms() if atom.name == 'CA' ]

    center_atoms += chain_atoms

#atom_list = [x for x in all_atom_list if x not in center_atoms]

for j in args.chain2 :
    atom_list = [atom for atom in structure[0][j].get_atoms() if atom.name == 'CA' ]
    ns = NeighborSearch(atom_list)

    nearby_residues = {res for center_atom in center_atoms
                        for res in ns.search(center_atom.coord, 8.5, 'R')}

    print "\nNeighbor residues in chain ", j, ": \n"
    print sorted(res.id[1] for res in nearby_residues)

    pymol_command = "show spheres, chain " + j + " and resi "

    for m in sorted(res.id[1] for res in nearby_residues):
        pymol_command = pymol_command + str(m) + "+"

    print pymol_command[:-1] + " and name CA \n" 

SCALAR_EXPRESSION %s_o = %s_startenergy + 4.91
SCALAR_EXPRESSION %s_s = %s
SCALAR_EXPRESSION %s_exp = exp(%s_s*(%s_currentenergy-%s_o))
SCALAR_EXPRESSION %s_k = %s
SCALAR_EXPRESSION %s_sig = (1-%s_k)+(%s_k/(1+%s_exp))

""") %(structure_id, structurefilename, correspondencefilename, resfilename, structure_id, startenergy, structure_id, structure_id, structure_id, s_value, structure_id, structure_id, structure_id, structure_id, structure_id, k_value, structure_id, structure_id, structure_id, structure_id)
		fitnessfile.write(outstring)

		fitnessstring = fitnessstring + str("*%s_sig") % (structure_id)
		

		nucleosome = structure[0]

		atom_list = Selection.unfold_entities(nucleosome, 'A') # A for atoms
		neighbor_search = NeighborSearch(atom_list)

		contacts_list = neighbor_search.search_all(radius, level = 'R')
		
		repack_residues = []
		for contact in contacts_list:
			res1 = contact[0]
			res2 = contact[1]
			res1id = int(res1.get_id()[1])
			chain1 = res1.get_parent()
			chain1id = chain1.get_id()
			res2id = int(res2.get_id()[1])
			chain2 = res2.get_parent()
			chain2id = chain2.get_id()
			
			res1_in_patch = False

def get_interacting_residues(model, r_cutoff=5, skip_hetatm_chains=True):
    """Return residue-residue interactions between all chains in `model`.

    Parameters
    ----------
    model : biopython.Model
        Model to analyse.

    Returns
    -------
    dict
        A dictionary of interactions between chains i (0..n-1) and j (i+1..n).
        Keys are (chain_idx, chain_id, residue_idx, residue_resnum, residue_amino_acid) tuples.
        (e.g. (0, 'A', 0, '0', 'M'), (0, 1, '2', 'K'), ...)
        Values are a list of tuples having the same format as the keys.

    Examples
    --------
    You can reverse the order of keys and values like this::

        complement = dict()
        for key, values in get_interacting_chains(model):
            for value in values:
                complement.setdefault(value, set()).add(key)


    You can get a list of all interacting chains using this command::

        {(key[0], value[0])
         for (key, values) in get_interacting_chains(model).items()
         for value in values}

    """
    from Bio.PDB import NeighborSearch

    interactions_between_chains = dict()

    # Chain 1
    for chain_1_idx, chain_1 in enumerate(model):
        if skip_hetatm_chains and chain_is_hetatm(chain_1):
            message = (
                "Skipping chain_1 with idx {} because it contains only hetatms."
                .format(chain_1_idx)
            )
            logger.debug(message)
            continue
        chain_1_residue_ids = get_aa_residues(chain_1)

        # Chain 2
        for j, chain_2 in enumerate(model.child_list[chain_1_idx + 1:]):
            chain_2_idx = chain_1_idx + 1 + j
            if skip_hetatm_chains and chain_is_hetatm(chain_2):
                message = (
                    "Skipping chain_2 with idx {} because it contains only hetatms."
                    .format(chain_2_idx)
                )
                logger.debug(message)
                continue
            chain_2_residue_ids = get_aa_residues(chain_2)
            ns = NeighborSearch(list(chain_2.get_atoms()))

            # Residue 1
            for residue_1 in chain_1:
                try:
                    residue_1_idx = chain_1_residue_ids.index(residue_1.id)
                except ValueError:
                    continue
                residue_1_resnum = str(residue_1.id[1]) + residue_1.id[2].strip()
                residue_1_aa = convert_aa(residue_1.resname, quiet=True)
                residue_1_key = (
                    chain_1_idx, chain_1.id, residue_1_idx, residue_1_resnum, residue_1_aa
                )
                interacting_residues = set()
                for atom_1 in residue_1:
                    interacting_residues.update(ns.search(atom_1.get_coord(), r_cutoff, 'R'))

                # Residue 2
                interacting_residue_ids = []
                for residue_2 in interacting_residues:
                    try:
                        residue_2_idx = chain_2_residue_ids.index(residue_2.id)
                    except ValueError:
                        continue
                    residue_2_resnum = str(residue_2.id[1]) + residue_2.id[2].strip()
                    residue_2_aa = convert_aa(residue_2.get_resname(), quiet=True)
                    residue_2_key = (
                        chain_2_idx, chain_2.id, residue_2_idx, residue_2_resnum, residue_2_aa
                    )
                    interacting_residue_ids.append(residue_2_key)
                if interacting_residue_ids:
                    interactions_between_chains\
                        .setdefault(residue_1_key, set())\
                        .update(interacting_residue_ids)

    return interactions_between_chains

for residue in structure.get_residues():
    if residue.resname == args.ligand:
        ligand = residue
        break

if not ligand:
    print('[!!] Ligand residue \'{0}\' not found in structure'.format(args.ligand), file=sys.stderr)
    sys.exit(1)

# Calculate center of mass of the ligand
ligand_com = map(lambda x: sum(x)/len(x), zip(*[at.coord for at in ligand]))
ligand_com = np.asarray(ligand_com, dtype=np.float32)

# Calculate neighbors considering only aminoacid/nucleotide atoms (excl. waters, other ligands, etc)
sel_atoms = [at for at in structure.get_atoms() if at.parent.id[0] == ' ']
ns = NeighborSearch(sel_atoms)
neighbors = ns.search(ligand_com, 10.0, level='R') # 10A radius, return residues

# Calculate residue closer to each ligand atom and the respective distance
ligand_atoms = ligand.child_list
min_dist_list, _seen = [], set()

for l_at in ligand_atoms:
    distances = []
    for residue in neighbors:
        for r_at in residue:
            distances.append((r_at, l_at, r_at - l_at))

    distances.sort(key=lambda x: x[-1])
    min_dist = distances[0]
    # One restraint per residue to keep the number of restraints small

    def _build_interface(self, model, id, threshold, rsa_calculation, rsa_threshold, include_waters=False, *chains):
        """
        Return the interface of a model
        """

        self.threshold=threshold

        # Recover chain list from initial unpacking
        chain_list = self.chain_list

        # Unfold atom list
        atom_list = []
        for c in model:
            if c.id in chain_list:
                atom_list.extend(Selection.unfold_entities(c,'A'))

        # Using of NeighborSearch class in order to get the list of all residues at least than
        # the threshold distance of each others
        ns=NeighborSearch(atom_list)
        pairs=ns.search_all(threshold, 'R')

        if not pairs:
            raise ValueError("No atoms found in the interface")        

        # Selection of residues pairs
        # 1. Exclude water contacts
        # 2. Filter same-chain contacts
        # 3. Filter user-defined chain pairs

        uniq_pairs=[]

        for pair in pairs:
             
            pair_resnames = (pair[0].resname, pair[1].resname)
            pair_chains = (pair[0].parent.id, pair[1].parent.id)

            if (not include_waters and 'HOH' in pair_resnames) or (pair_chains[0] == pair_chains[1]):
                continue

            if not (chains and not (pair_chains in chains)):
                uniq_pairs.append(pair)

        # Build the Interface
        # 1. Iterate over the pair list
        # 2. Add residues.

        for resA, resB in uniq_pairs:
            if resA not in self.interface:
                self._add_residue(resA)
            if resB not in self.interface:
                self._add_residue(resB)
                
        # Accessible surface area calculated for each residue
        # if naccess setup on user computer and rsa_calculation
        # argument is TRUE
        if rsa_calculation and os.system('which naccess') == 0:
            rsa_pairs=self._rsa_calculation(model, chain_list, rsa_threshold)
            
        for res in rsa_pairs:
            if res not in self.interface:
                self._add_residue(res)
        self._secondary_structure(model)
        #interface=uniq_pairs
        self.interface.uniq_pairs=uniq_pairs

        PDB = pdb_biomol_match.group(1)
        BIOMOL_ID = pdb_biomol_match.group(2)

    pdb_rcsb_asm_match = re.search(PDB_RCSB_ASM_REGEX, INPUT_FILENAME)
    if pdb_rcsb_asm_match:
        PDB = pdb_rcsb_asm_match.group(1)
        BIOMOL_ID = pdb_rcsb_asm_match.group(2)

    # LOAD STRUCTURE
    structure = PDBParser().get_structure('structure', INPUT_FILE)
    structure_atoms = list(structure.get_atoms())

    logging.info('Loaded PDB structure (BioPython).')

    # CONSTRUCT KDTREE
    neighborsearch = NeighborSearch(structure_atoms)

    logging.info('Constructured NeighborSearch.')

    # GET INTERACTIONS
    logging.info('Calculating interactions...')
    for interaction_level in 'ARC':

        if interaction_level in OUTPUTS:

            logging.info('Calculating interactions for {}s...'.format(
                LEVEL_MAP[interaction_level]))

            pairs = neighborsearch.search_all(INTERACTION_THRESHOLD,
                                              level=interaction_level)