Zoom in too ligand and its interactions.
192047
def zoom_to_ligand(self):
        """"""
        cmd.center(self.ligname)
        cmd.orient(self.ligname)
        cmd.turn('x', 110)  # If the ligand is aligned with the longest axis, aromatic rings are hidden
        if 'AllBSRes' in cmd.get_names("selections"):
            cmd.zoom('%s or AllBSRes' % self.ligname, 3)
        else:
            if self.object_exists(self.ligname):
                cmd.zoom(self.ligname, 3)
        cmd.origin(self.ligname)
1019199
def get_atom_mapping(self):
        """"""
        # Atom mappings
        smiles_to_pdb_mapping = self.bindingsite.xpath('mappings/smiles_to_pdb/text()')
        if smiles_to_pdb_mapping == []:
            self.mappings = {'smiles_to_pdb': None, 'pdb_to_smiles': None}
        else:
            smiles_to_pdb_mapping = {int(y[0]): int(y[1]) for y in [x.split(':')
                                                                    for x in smiles_to_pdb_mapping[0].split(',')]}
            self.mappings = {'smiles_to_pdb': smiles_to_pdb_mapping}
            self.mappings['pdb_to_smiles'] = {v: k for k, v in self.mappings['smiles_to_pdb'].items()}
1089600
def show_metal(self):
        """"""
        metal_complexes = self.plcomplex.metal_complexes
        if not len(metal_complexes) == 0:
            self.select_by_ids('Metal-M', self.metal_ids)
            for metal_complex in metal_complexes:
                cmd.select('tmp_m', 'id %i' % metal_complex.metal_id)
                cmd.select('tmp_t', 'id %i' % metal_complex.target_id)
                if metal_complex.location == 'water':
                    cmd.select('Metal-W', 'Metal-W or id %s' % metal_complex.target_id)
                if metal_complex.location.startswith('protein'):
                    cmd.select('tmp_t', 'tmp_t & %s' % self.protname)
                    cmd.select('Metal-P', 'Metal-P or (id %s & %s)' % (metal_complex.target_id, self.protname))
                if metal_complex.location == 'ligand':
                    cmd.select('tmp_t', 'tmp_t & %s' % self.ligname)
                    cmd.select('Metal-L', 'Metal-L or (id %s & %s)' % (metal_complex.target_id, self.ligname))
                cmd.distance('MetalComplexes', 'tmp_m', 'tmp_t')
                cmd.delete('tmp_m or tmp_t')
        if self.object_exists('MetalComplexes'):
            cmd.set('dash_color', 'violetpurple', 'MetalComplexes')
            cmd.set('dash_gap', 0.5, 'MetalComplexes')
            # Show water molecules for metal complexes
            cmd.show('spheres', 'Metal-W')
            cmd.color('lightblue', 'Metal-W')
123913
def find_charged(self, mol):
        """"""
        data = namedtuple('pcharge', 'atoms atoms_orig_idx type center restype resnr reschain')
        a_set = []
        # Iterate through all residue, exclude those in chains defined as peptides
        for res in [r for r in pybel.ob.OBResidueIter(mol.OBMol) if not r.GetChain() in config.PEPTIDES]:
            if config.INTRA is not None:
                if res.GetChain() != config.INTRA:
                    continue
            a_contributing = []
            a_contributing_orig_idx = []
            if res.GetName() in ('ARG', 'HIS', 'LYS'):  # Arginine, Histidine or Lysine have charged sidechains
                for a in pybel.ob.OBResidueAtomIter(res):
                    if a.GetType().startswith('N') and res.GetAtomProperty(a, 8) \
                            and not self.Mapper.mapid(a.GetIdx(), mtype='protein') in self.altconf:
                        a_contributing.append(pybel.Atom(a))
                        a_contributing_orig_idx.append(self.Mapper.mapid(a.GetIdx(), mtype='protein'))
                if not len(a_contributing) == 0:
                    a_set.append(data(atoms=a_contributing,
                                      atoms_orig_idx=a_contributing_orig_idx,
                                      type='positive',
                                      center=centroid([ac.coords for ac in a_contributing]),
                                      restype=res.GetName(),
                                      resnr=res.GetNum(),
                                      reschain=res.GetChain()))
            if res.GetName() in ('GLU', 'ASP'):  # Aspartic or Glutamic Acid
                for a in pybel.ob.OBResidueAtomIter(res):
                    if a.GetType().startswith('O') and res.GetAtomProperty(a, 8) \
                            and not self.Mapper.mapid(a.GetIdx(), mtype='protein') in self.altconf:
                        a_contributing.append(pybel.Atom(a))
                        a_contributing_orig_idx.append(self.Mapper.mapid(a.GetIdx(), mtype='protein'))
                if not len(a_contributing) == 0:
                    a_set.append(data(atoms=a_contributing,
                                      atoms_orig_idx=a_contributing_orig_idx,
                                      type='negative',
                                      center=centroid([ac.coords for ac in a_contributing]),
                                      restype=res.GetName(),
                                      resnr=res.GetNum(),
                                      reschain=res.GetChain()))
        return a_set
13963
def is_within(self, cutoff_dist, point, ligands=True):
        """"""
        return find_atoms_within_distance(self.get_atoms(ligands=ligands), cutoff_dist, point)
472203
def _approximate_unkown_bond_lengths(self):
        """"""
        dataset = self.lengths[BOND_SINGLE]
        for n1 in periodic.iter_numbers():
            for n2 in periodic.iter_numbers():
                if n1 <= n2:
                    pair = frozenset([n1, n2])
                    atom1 = periodic[n1]
                    atom2 = periodic[n2]
                    #if (pair not in dataset) and hasattr(atom1, "covalent_radius") and hasattr(atom2, "covalent_radius"):
                    if (pair not in dataset) and (atom1.covalent_radius is not None) and (atom2.covalent_radius is not None):
                        dataset[pair] = (atom1.covalent_radius + atom2.covalent_radius)
308190
def find_the_closest_atoms(self,topology):
        

        # The measurements are made to ligand molecule without hydrogen atoms (ligand_noH) because the
        # hydrogen atoms are not plotted in the final graph
        self.universe.load_new(topology)
        self.universe.ligand_noH = self.universe.ligand.select_atoms("not name H*")
        ligand_positions = self.universe.ligand_noH.positions

        for residue in self.dict_of_plotted_res.keys():
            residue_selection = self.universe.select_atoms("resname "+residue[0]+" and resid "+residue[1]+" and segid "+ residue[2])
            residue_positions = residue_selection.positions
            dist_array = MDAnalysis.analysis.distances.distance_array(ligand_positions,residue_positions)
            min_values_per_atom={}
            i=0
            for atom in self.universe.ligand_noH:
                min_values_per_atom[atom.name]=dist_array[i].min()
                i+=1
            sorted_min_values = sorted(min_values_per_atom.items(), key=operator.itemgetter(1))
            self.closest_atoms[residue]=[(sorted_min_values[0][0],sorted_min_values[0][1])]
148047
def run(self):
        """"""

        # Load compound information
        def compound_name(id):
            if id not in self._model.compounds:
                return id
            return self._model.compounds[id].properties.get('name', id)

        epsilon = self._args.epsilon

        if self._args.unrestricted:
            # Allow all exchange reactions with no flux limits
            for reaction in self._mm.reactions:
                if self._mm.is_exchange(reaction):
                    del self._mm.limits[reaction].bounds

        loop_removal = self._get_loop_removal_option()
        enable_tfba = loop_removal == 'tfba'
        enable_fastcore = self._args.fastcore

        if enable_tfba and enable_fastcore:
            self.argument_error(
                'Using Fastcore with thermodynamic constraints'
                ' is not supported!')
        start_time = time.time()

        if enable_fastcore:
            solver = self._get_solver()
            try:
                inconsistent = set(fastcore.fastcc(
                    self._mm, epsilon, solver=solver))
            except fluxanalysis.FluxBalanceError as e:
                self.report_flux_balance_error(e)
        else:
            if enable_tfba:
                solver = self._get_solver(integer=True)
            else:
                solver = self._get_solver()

            if self._args.reduce_lp:
                logger.info('Running with reduced number of LP problems.')
                try:
                    inconsistent = set(
                        fluxanalysis.consistency_check(
                            self._mm, self._mm.reactions, epsilon,
                            tfba=enable_tfba, solver=solver))
                except fluxanalysis.FluxBalanceError as e:
                    self.report_flux_balance_error(e)
            else:
                logger.info('Using flux bounds to determine consistency.')
                try:
                    inconsistent = set(self._run_fva_fluxcheck(
                        self._mm, solver, enable_tfba, epsilon))
                except FluxCheckFVATaskError:
                    self.report_flux_balance_error()

        logger.info('Solving took {:.2f} seconds'.format(
            time.time() - start_time))

        # Count the number of reactions that are fixed at zero. While these
        # reactions are still inconsistent, they are inconsistent because they
        # have been explicitly disabled.
        disabled_exchange = 0
        disabled_internal = 0

        count_exchange = 0
        total_exchange = 0

        count_internal = 0
        total_internal = 0

        # Print result
        for reaction in sorted(self._mm.reactions):
            disabled = self._mm.limits[reaction].bounds == (0, 0)

            if self._mm.is_exchange(reaction):
                total_exchange += 1
                count_exchange += int(reaction in inconsistent)
                disabled_exchange += int(disabled)
            else:
                total_internal += 1
                count_internal += int(reaction in inconsistent)
                disabled_internal += int(disabled)

            if reaction in inconsistent:
                rx = self._mm.get_reaction(reaction)
                rxt = rx.translated_compounds(compound_name)
                print('{}\t{}'.format(reaction, rxt))

        logger.info('Model has {}/{} inconsistent internal reactions'
                    ' ({} disabled by user)'.format(
                        count_internal, total_internal, disabled_internal))
        logger.info('Model has {}/{} inconsistent exchange reactions'
                    ' ({} disabled by user)'.format(
                        count_exchange, total_exchange, disabled_exchange))
971180
def run(self):
        """"""

        # Load compound information
        def compound_name(id):
            if id not in self._model.compounds:
                return id
            return self._model.compounds[id].properties.get('name', id)

        # Calculate penalty if penalty file exists
        penalties = {}
        if self._args.penalty is not None:
            for line in self._args.penalty:
                line, _, comment = line.partition('#')
                line = line.strip()
                if line == '':
                    continue
                rxnid, penalty = line.split(None, 1)
                penalties[rxnid] = float(penalty)

        core = set(self._mm.reactions)

        solver = self._get_solver(integer=True)
        default_comp = self._model.default_compartment
        epsilon = self._args.epsilon
        v_max = float(self._model.default_flux_limit)

        blocked = set()
        for compound in self._args.compound:
            if compound.compartment is None:
                compound = compound.in_compartment(default_comp)
            blocked.add(compound)

        if len(blocked) > 0:
            logger.info('Unblocking compounds: {}...'.format(
                ', '.join(text_type(c) for c in sorted(blocked))))
        else:
            logger.info(
                'Unblocking all compounds in model. Use --compound option to'
                ' unblock specific compounds.')
            blocked = set(self._mm.compounds)

        exclude = set()
        if self._model.biomass_reaction is not None:
            exclude.add(self._model.biomass_reaction)

        # Add exchange and transport reactions to database
        model_complete, weights = create_extended_model(
            self._model,
            db_penalty=self._args.db_penalty,
            ex_penalty=self._args.ex_penalty,
            tp_penalty=self._args.tp_penalty,
            penalties=penalties)

        implicit_sinks = not self._args.no_implicit_sinks

        logger.info('Searching for reactions to fill gaps')
        try:
            added_reactions, no_bounds_reactions = gapfill(
                model_complete, core, blocked, exclude, solver=solver,
                epsilon=epsilon, v_max=v_max, weights=weights,
                implicit_sinks=implicit_sinks,
                allow_bounds_expansion=self._args.allow_bounds_expansion)
        except GapFillError as e:
            self._log_epsilon_and_fail(epsilon, e)

        for reaction_id in sorted(self._mm.reactions):
            rx = self._mm.get_reaction(reaction_id)
            rxt = rx.translated_compounds(compound_name)
            print('{}\t{}\t{}\t{}'.format(reaction_id, 'Model', 0, rxt))

        for rxnid in sorted(added_reactions):
            rx = model_complete.get_reaction(rxnid)
            rxt = rx.translated_compounds(compound_name)
            print('{}\t{}\t{}\t{}'.format(
                rxnid, 'Add', weights.get(rxnid, 1), rxt))

        for rxnid in sorted(no_bounds_reactions):
            rx = model_complete.get_reaction(rxnid)
            rxt = rx.translated_compounds(compound_name)
            print('{}\t{}\t{}\t{}'.format(
                rxnid, 'Remove bounds', weights.get(rxnid, 1), rxt))
708006
def find_metal_binding(self, mol):
        """"""
        data = namedtuple('metal_binding', 'atom atom_orig_idx type restype resnr reschain location')
        a_set = []
        for res in pybel.ob.OBResidueIter(mol.OBMol):
            restype, reschain, resnr = res.GetName().upper(), res.GetChain(), res.GetNum()
            if restype in ['ASP', 'GLU', 'SER', 'THR', 'TYR']:  # Look for oxygens here
                for a in pybel.ob.OBResidueAtomIter(res):
                    if a.GetType().startswith('O') and res.GetAtomProperty(a, 8) \
                            and not self.Mapper.mapid(a.GetIdx(), mtype='protein') in self.altconf:
                        atom_orig_idx = self.Mapper.mapid(a.GetIdx(), mtype=self.mtype, bsid=self.bsid)
                        a_set.append(data(atom=pybel.Atom(a), atom_orig_idx=atom_orig_idx, type='O', restype=restype,
                                          resnr=resnr, reschain=reschain,
                                          location='protein.sidechain'))
            if restype == 'HIS':  # Look for nitrogen here
                for a in pybel.ob.OBResidueAtomIter(res):
                    if a.GetType().startswith('N') and res.GetAtomProperty(a, 8) \
                            and not self.Mapper.mapid(a.GetIdx(), mtype='protein') in self.altconf:
                        atom_orig_idx = self.Mapper.mapid(a.GetIdx(), mtype=self.mtype, bsid=self.bsid)
                        a_set.append(data(atom=pybel.Atom(a), atom_orig_idx=atom_orig_idx, type='N', restype=restype,
                                          resnr=resnr, reschain=reschain,
                                          location='protein.sidechain'))
            if restype == 'CYS':  # Look for sulfur here
                for a in pybel.ob.OBResidueAtomIter(res):
                    if a.GetType().startswith('S') and res.GetAtomProperty(a, 8) \
                            and not self.Mapper.mapid(a.GetIdx(), mtype='protein') in self.altconf:
                        atom_orig_idx = self.Mapper.mapid(a.GetIdx(), mtype=self.mtype, bsid=self.bsid)
                        a_set.append(data(atom=pybel.Atom(a), atom_orig_idx=atom_orig_idx, type='S', restype=restype,
                                          resnr=resnr, reschain=reschain,
                                          location='protein.sidechain'))
            for a in pybel.ob.OBResidueAtomIter(res):  # All main chain oxygens
                if a.GetType().startswith('O') and res.GetAtomProperty(a, 2) \
                        and not self.Mapper.mapid(a.GetIdx(), mtype='protein') in self.altconf and restype != 'HOH':
                    atom_orig_idx = self.Mapper.mapid(a.GetIdx(), mtype=self.mtype, bsid=self.bsid)
                    a_set.append(data(atom=pybel.Atom(a), atom_orig_idx=atom_orig_idx, type='O', restype=res.GetName(),
                                      resnr=res.GetNum(), reschain=res.GetChain(),
                                      location='protein.mainchain'))
        return a_set
64581
def hydrophobic_atoms(self, all_atoms):
        """"""
        atom_set = []
        data = namedtuple('hydrophobic', 'atom orig_atom orig_idx')
        atm = [a for a in all_atoms if a.atomicnum == 6 and set([natom.GetAtomicNum() for natom
                                                                 in pybel.ob.OBAtomAtomIter(a.OBAtom)]).issubset(
            {1, 6})]
        for atom in atm:
            orig_idx = self.Mapper.mapid(atom.idx, mtype=self.mtype, bsid=self.bsid)
            orig_atom = self.Mapper.id_to_atom(orig_idx)
            if atom.idx not in self.altconf:
                atom_set.append(data(atom=atom, orig_atom=orig_atom, orig_idx=orig_idx))
        return atom_set
1252182
def from_tuples_dict(pair_dict):
        '''
        '''
        lm = LigandMap()
        for k, v in pair_dict.iteritems():
            lm.add(k[0], k[1], v[0], v[1])
        return lm
987272
private void reflect(int v, IBond bond) {
        visited[v] = true;
        IAtom atom = container.getAtom(v);
        atom.setPoint2d(reflect(atom.getPoint2d(), bond));
        for (int w : graph[v]) {
            if (!visited[w]) reflect(w, bond);
        }
    }
69116
public void calculateOverlapsAndReduce(IAtomContainer molecule1, IAtomContainer molecule2, boolean shouldMatchBonds)
            throws CDKException {

        setSource(molecule1);
        setTarget(molecule2);

        setMappings(new ArrayList<Map<Integer, Integer>>());

        if ((getSource().getAtomCount() == 1) || (getTarget().getAtomCount() == 1)) {
            List<CDKRMap> overlaps = CDKMCS.checkSingleAtomCases(getSource(), getTarget());
            int nAtomsMatched = overlaps.size();
            nAtomsMatched = (nAtomsMatched > 0) ? 1 : 0;
            if (nAtomsMatched > 0) {
                /* UnComment this to get one Unique Mapping */
                //List reducedList = removeRedundantMappingsForSingleAtomCase(overlaps);
                //int counter = 0;
                identifySingleAtomsMatchedParts(overlaps, getSource(), getTarget());

            }

        } else {
            List<List<CDKRMap>> overlaps = CDKMCS.search(getSource(), getTarget(), new BitSet(), new BitSet(), true,
                    true, shouldMatchBonds);

            List<List<CDKRMap>> reducedList = removeSubGraph(overlaps);
            Stack<List<CDKRMap>> allMaxOverlaps = getAllMaximum(reducedList);
            while (!allMaxOverlaps.empty()) {
                //                System.out.println("source: " + source.getAtomCount() + ", target: " + target.getAtomCount() + ", overl: " + allMaxOverlaps.peek().size());
                List<List<CDKRMap>> maxOverlapsAtoms = makeAtomsMapOfBondsMap(allMaxOverlaps.peek(), getSource(),
                        getTarget());
                //                System.out.println("size of maxOverlaps: " + maxOverlapsAtoms.size());
                identifyMatchedParts(maxOverlapsAtoms, getSource(), getTarget());
                //                identifyMatchedParts(allMaxOverlaps.peek(), source, target);
                allMaxOverlaps.pop();
            }
        }

        FinalMappings.getInstance().set(getMappings());

    }
1341521
def analyse_ligand_sasa(self):
		""""""
		i=0
		start = timer()
		if self.trajectory == []:
			self.trajectory = [self.topology_data.universe.filename]
		try:
			for traj in self.trajectory:
				new_traj = mdtraj.load(traj,top=self.topology_data.universe.filename)
				#Analyse only non-H ligand
				ligand_slice = new_traj.atom_slice(atom_indices=self.topology_data.universe.ligand_noH.ids)

				self.sasa = mdtraj.shrake_rupley(ligand_slice)
				self.atom_sasa[i]=self.assign_per_atom_sasa()
				i+=1
			self.total_sasa = self.get_total_per_atom_sasa()
		except KeyError as e:
			print "WARNING: SASA analysis cannot be performed due to incorrect atom names in"
			print "the topology ", e

		print "SASA: "+str(timer()-start)
849539
def characterize_complex(self, ligand):
        """"""

        single_sites = []
        for member in ligand.members:
            single_sites.append(':'.join([str(x) for x in member]))
        site = ' + '.join(single_sites)
        site = site if not len(site) > 20 else site[:20] + '...'
        longname = ligand.longname if not len(ligand.longname) > 20 else ligand.longname[:20] + '...'
        ligtype = 'Unspecified type' if ligand.type == 'UNSPECIFIED' else ligand.type
        ligtext = "\n%s [%s] -- %s" % (longname, ligtype, site)
        if ligtype == 'PEPTIDE':
            ligtext = '\n Chain %s [PEPTIDE / INTER-CHAIN]' % ligand.chain
        if ligtype == 'INTRA':
            ligtext = "\n Chain %s [INTRA-CHAIN]" % ligand.chain
        any_in_biolip = len(set([x[0] for x in ligand.members]).intersection(config.biolip_list)) != 0
        write_message(ligtext)
        write_message('\n' + '-' * len(ligtext) + '\n')

        if ligtype not in ['POLYMER', 'DNA', 'ION', 'DNA+ION', 'RNA+ION', 'SMALLMOLECULE+ION'] and any_in_biolip:
            write_message('may be biologically irrelevant\n', mtype='info', indent=True)

        lig_obj = Ligand(self, ligand)
        cutoff = lig_obj.max_dist_to_center + config.BS_DIST
        bs_res = self.extract_bs(cutoff, lig_obj.centroid, self.resis)
        # Get a list of all atoms belonging to the binding site, search by idx
        bs_atoms = [self.atoms[idx] for idx in [i for i in self.atoms.keys()
                                                if self.atoms[i].OBAtom.GetResidue().GetIdx() in bs_res]
                    if idx in self.Mapper.proteinmap and self.Mapper.mapid(idx, mtype='protein') not in self.altconf]
        if ligand.type == 'PEPTIDE':
            # If peptide, don't consider the peptide chain as part of the protein binding site
            bs_atoms = [a for a in bs_atoms if a.OBAtom.GetResidue().GetChain() != lig_obj.chain]
        if ligand.type == 'INTRA':
            # Interactions within the chain
            bs_atoms = [a for a in bs_atoms if a.OBAtom.GetResidue().GetChain() == lig_obj.chain]
        bs_atoms_refined = []

        # Create hash with BSRES -> (MINDIST_TO_LIG, AA_TYPE)
        # and refine binding site atom selection with exact threshold
        min_dist = {}
        for r in bs_atoms:
            bs_res_id = ''.join([str(whichresnumber(r)), whichchain(r)])
            for l in ligand.mol.atoms:
                distance = euclidean3d(r.coords, l.coords)
                if bs_res_id not in min_dist:
                    min_dist[bs_res_id] = (distance, whichrestype(r))
                elif min_dist[bs_res_id][0] > distance:
                    min_dist[bs_res_id] = (distance, whichrestype(r))
                if distance <= config.BS_DIST and r not in bs_atoms_refined:
                    bs_atoms_refined.append(r)
        num_bs_atoms = len(bs_atoms_refined)
        write_message('Binding site atoms in vicinity (%.1f A max. dist: %i).\n' % (config.BS_DIST, num_bs_atoms),
                      indent=True)

        bs_obj = BindingSite(bs_atoms_refined, self.protcomplex, self, self.altconf, min_dist, self.Mapper)
        pli_obj = PLInteraction(lig_obj, bs_obj, self)
        self.interaction_sets[ligand.mol.title] = pli_obj
1282083
def detect_aromatic_rings_in_ligand(self):
        """"""
        self.ligrings = {}
        try:
            ring_info = self.topology_data.mol.GetRingInfo()
            self.ligand_ring_num = ring_info.NumRings()
        except Exception as e:
            m = Chem.MolFromPDBFile("lig.pdb")
            ring_info = m.GetRingInfo()
            self.ligand_ring_num = ring_info.NumRings()
        i=0
        for ring in range(self.ligand_ring_num):
            if 4 < len(ring_info.AtomRings()[ring]) <= 6 and  False not in [self.topology_data.mol.GetAtomWithIdx(x).GetIsAromatic() for x in ring_info.AtomRings()[ring]]: #narrow ring definition
                atom_ids_in_ring = []
                for atom in ring_info.AtomRings()[ring]:
                    atom_ids_in_ring.append(self.topology_data.universe.ligand.atoms[atom].name)
                self.ligrings[i]=atom_ids_in_ring
                i+=1