Python pybel.readfile函数代码示例

def do_for_each(name, num_threads):
    _format = name[name.rfind(".")+1:]
    basename = name[:name.rfind(".")]
    out_name = basename + ".num"
    pdb_id = basename.split("_")[-1]
    pdb_name = "/home/xmluo/jlpeng/cMet/pdb_protein/%s_protein.pdb"%pdb_id
    global ligand,protein,outf,ap_iter
    ligand = pybel.readfile(_format, name).next()
    protein = pybel.readfile("pdb", pdb_name).next()
    outf = open(out_name, "w")
    #make_pairs()
    #print "after creating pair_queue"
    ap_iter = AtomPairs()
    threads = []
    for i in xrange(num_threads):
        threads.append(AtomPairWalker())
    print "len(threads)=%d"%len(threads)
    for t in threads:
        t.start()
    for t in threads:
        t.join()
        if t.isAlive():
            print "thread %s is alive"%t.getName()
        else:
            print "thread %s is dead"%t.getName()
    outf.close()

def mp_helper(file_one, file_two, iteration_depth = 1):
    """
        Helper function for the multiprocessing library.
        Two fragment files gets passed and we merge all against all molecules in that two files.
    """
    results = list()
    fragments = list()
    for mol_one in pybel.readfile( 'smi', file_one ):
        for i,mol_two in enumerate(pybel.readfile( 'smi', file_two )):
            #print 'merge:', str(mol_one).strip(), str(mol_two).strip()
            result, fragment = merge(mol_two,mol_one, options, iteration_depth)
            if result:
                #print '\tr', result
                results.extend( result )
            if fragment:
                #print '\tf', fragment
                fragments.extend( fragment )

    fragment_return, molecule_return = None, None
    if fragments:
        fragment_file = tempfile.NamedTemporaryFile(dir=temp_dir, delete=False)
        fragment_file.writelines( fragments )
        fragment_file.close()
        fragment_return = fragment_file.name
    if results:
        result_file = tempfile.NamedTemporaryFile(dir=temp_dir, delete=False)
        result_file.writelines( results )
        result_file.close()
        molecule_return = result_file.name
    return molecule_return, fragment_return

def dock_ligands_to_receptors(docking_dir, worker_pool=False, exhaustiveness=None, chosen_receptor=None, restrict_box=True):
  subdirs = glob.glob(os.path.join(docking_dir, '*/'))
  for subdir in subdirs:
    subdir = subdir.rstrip('/')
    receptor_name = os.path.basename(subdir)
    if chosen_receptor is not None and chosen_receptor != receptor_name:
        continue
    print("receptor name = %s" % receptor_name)
    receptor_filename = os.path.join(subdir, "%s.pdbqt" % receptor_name)
    if not os.path.exists(receptor_filename):
      continue

    print("Examining %s" % receptor_filename)

    receptor_pybel = next(pybel.readfile("pdb", 
        os.path.join(subdir, "%s.pdb" % receptor_name)))
    protein_centroid, protein_range = get_molecule_data(receptor_pybel)

    box_dims = protein_range + 5.0

    ligands = sorted(glob.glob(os.path.join(subdir, '*_prepared.pdbqt')))
    print("Num ligands = %d" % len(ligands))


    dock_ligand_to_receptor_partial = partial(dock_ligand_to_receptor, receptor_filename=receptor_filename,
                                              protein_centroid=protein_centroid, box_dims=box_dims,
                                              subdir=subdir, exhaustiveness=exhaustiveness)

    if restrict_box:
      active_ligand = ""
      for ligand in ligands:
        if "CHEM" in ligand:
          active_ligand = ligand
          break

      print("Docking to %s first to ascertain centroid and box dimensions" % active_ligand)

      out_pdb_qt = dock_ligand_to_receptor_partial(active_ligand)
      ligand_pybel = next(pybel.readfile("pdbqt", 
                                    out_pdb_qt))
      ligand_centroid, _ = get_molecule_data(ligand_pybel)
      print("Protein centroid = %s" %(str(protein_centroid)))
      print("Ligand centroid = %s" %(str(ligand_centroid)))
      box_dims = np.array([20., 20., 20.])
      dock_ligand_to_receptor_partial = partial(dock_ligand_to_receptor, receptor_filename=receptor_filename,
                                          protein_centroid=ligand_centroid, box_dims=box_dims,
                                          subdir=subdir, exhaustiveness=exhaustiveness)

      print("Finished docking to %s, docking to remainder of ligands now." % active_ligand)

    if worker_pool is False:
      for i, ligand_file in enumerate(ligands):
        a = time.time()
        dock_ligand_to_receptor_partial(ligand)
        print("took %f seconds to dock single ligand." %(time.time() - a))
    else:
      print("parallelizing docking over worker pool")

      worker_pool.map_sync(dock_ligand_to_receptor_partial, ligands)

    def helper_run(self, result):
        result_pdbqt = result.output().path
        result_lig = pybel.readfile("pdbqt", result_pdbqt).next()
        native_lig = pybel.readfile("pdbqt",
                                    LigPdbqt(self.tname).output().path).next()
        result_rmsd = rmsd_between(native_lig, result_lig)

        return result_rmsd

def create_inchi(job, chemid, xyz_file=''):
    if xyz_file == '':
        xyz_file = os.path.expanduser(job) + 'xyz/' + chemid + '.xyz'
    obmol = list(pybel.readfile('xyz', xyz_file))[0]
    try:
        obmol = list(pybel.readfile('xyz', xyz_file))[0]
    except NameError:
        logging.error('Pybel is not installed or loaded correctly.')
        sys.exit()

    return obmol.write("inchi", opt={'T': 'nostereo'}).split()[0]

def prepare_dsstox_dataset(root, name, dest=None, overwrite=False):
    """ This method bootstraps the analysis of DSSTox data.
       - Rename the compounds
       - Merge train/test
       - Generate 3D conformations
       - Save "master" and "saliviewer" tables
       - Redirects stdout/stderr to a "prepare.log" file
    """
    if not dest: dest = root

    dataset_root = op.join(dest, name)
    dest_sdf = op.join(dataset_root, name + '.sdf')
    if op.exists(dest_sdf) and not overwrite:
        print '%s is already there and not overwriting requested' % dest_sdf
        return

    print 'Reading %s' % name
    train_mols = list(pybel.readfile('sdf', op.join(root, name + '_training.sdf')))
    test_mols = list(pybel.readfile('sdf', op.join(root, name + '_prediction.sdf')))

    print '\tCreating dataset root: %s' % dataset_root
    if not op.exists(dataset_root):
        os.makedirs(dataset_root)

    print '\tRenaming the compounds to keep track of the provenance'
    rename_mols_by_index(train_mols, name + '-train-')
    rename_mols_by_index(test_mols, name + '-test-')

    print '\tGenerating conformations'
    for mol in train_mols + test_mols:
        #Some molecules from mutagenicity produce segfault on make3D
        #See bug report at https://sourceforge.net/tracker/?func=detail&aid=3374324&group_id=40728&atid=428740
        #Train 3988: OC(=O)[[email protected]]1(C)CCC[[email protected]]2(C1CC[[email protected]]13C2CC[[email protected]](C3)([[email protected]]2(C1)OC2)O)C
        #Train 4205: CC(CCC[[email protected]]([[email protected]]1CC[[email protected]@H]2[[email protected]]1(C)CC[[email protected]]1[[email protected]]2CC2([[email protected]@H]3[[email protected]]1(C)CC[[email protected]@H](C3)Br)S(=O)(=O)CCS2(=O)=O)C)C
        #This kind of fatal errors are worrying, is there any robust way of controlling them in python/java? Will need to create one
        if not any(name in mol.title for name in ('train-3988', 'train-4205')):
            try:
                print 'Conformation for %s' % mol.title
                mol.make3D()
            except Exception:
                print 'Error computing a 3D conformation for %s' % mol.title

    print '\tSaving compounds'
    save_mols(train_mols + test_mols, dest_sdf)

    master_table = op.join(dataset_root, name + '-master.csv')
    print '\tCreating \"master\" table: %s' % master_table
    create_master_table(dest_sdf, master_table)

    sali_table = op.join(dataset_root, name + '-saliviewer.csv')
    print '\tCreating \"saliviewer\" table: %s' % sali_table
    create_saliviewer_input(master_table, sali_table)

 def __init__(self, filename, ext):
     self.data = {'name':'','index':'', 'metal':'', 'topology':'', 'parent':'',
             'atomic_info':'', 'bond_table':'', 'connectivity':'',
             'connect_flag':'', 'connect_sym':''}
     name = os.path.split(filename)[-1]
     self.name = clean(name, ext) 
     self.update(name=self.name)
     # may be a source of error.. untested
     if version_info.major >= 3:
         self.mol = next(pybel.readfile(ext, filename))
     else:
         self.mol = pybel.readfile(ext, filename).next()
     self._reset_formal_charges()

    def pred_foo(tname):
        try:
            geaux_sdf = os.path.join(GEAUX_OUTPUT, tname,
                                     tname + '_pred_0.sdf')
            native_sdf = VinaPath(tname).lig_sdf

            geaux_lig = pybel.readfile("sdf", geaux_sdf).next()
            native_lig = pybel.readfile("sdf", native_sdf).next()
            result_rmsd = rmsd_between(native_lig, geaux_lig)
            return result_rmsd
        except Exception as e:
            print(e)
            return None

    def caculateRMSD(self):
        vina_task = self.requires()
        predicted_pdbqt = vina_task.output().path
        predicted_mol = pybel.readfile('pdbqt', predicted_pdbqt).next()
        crystal_pdbqt = vina_task.lig_pdbqt
        crystal_mol = pybel.readfile('pdbqt', crystal_pdbqt).next()

        def rmsd(m1, m2):
            c1 = [a.coords for a in m1 if not a.OBAtom.IsHydrogen()]
            c2 = [a.coords for a in m2 if not a.OBAtom.IsHydrogen()]
            return dockedpose.rmsd(c1, c2)

        return rmsd(predicted_mol, crystal_mol)

 def helper_geauxdock(tname, version="0.7"):
     try:
         native_sdf = VinaPath(tname).lig_sdf
         geaux_sdf = os.path.join(
             ModelPath(tname,
                       version=version).work_dir,
             "{}_{}_0.sdf".format(tname, version))
         geaux_lig = pybel.readfile("sdf", geaux_sdf).next()
         native_lig = pybel.readfile("sdf", native_sdf).next()
         result_rmsd = rmsd_between(native_lig, geaux_lig)
         return result_rmsd
     except Exception as e:
         print(e)
         return None

 def run(self):
     data = {}
     for sdf_id in self.getSdfs():
         path = Path(sdf_id)
         lig = pybel.readfile('sdf', path.astex_sdf()).next()
         lig.removeh()
         prt = pybel.readfile('pdb', path.astex_pdb()).next()
         prt.removeh()
         lig_sz = len(lig.atoms)
         prt_sz = len(prt.atoms)
         data[sdf_id] = {"lig_sz": lig_sz,
                         "prt_sz": prt_sz}
     dset = pd.DataFrame(data)
     dset.to_csv(self.output().path)

def main():
    file_list=GetFileList()
    dc_fn='dist_splited_train_0525.dat'
    dc=dist.UnformatedInput(dc_fn)
    pts=dist.PTS(dc)
    for i in range(len(file_list)):
        lig_fn=file_list[i][0]
        pro_fn=file_list[i][1]
        pdbid=file_list[i][2]
        if os.path.isfile(lig_fn):
            lig=pybel.readfile('sdf',lig_fn).next()
            pro=pybel.readfile('pdb',pro_fn).next()
            score,scorev=GetFinger(lig,pro,pts,pdbid)
            print score,scorev

def main( args ):

    for infile in args.infile:
        file_extension = args.format or os.path.splitext( infile )[-1].lstrip('.')

        if not args.format and file_extension not in ['smi', 'sdf', 'inchi', 'mol']:
            sys.exit('Could not guess the format from the file extension please specify with the --format option.')

        molecules = pybel.readfile(file_extension, infile)
        for mol in molecules:
            args.outfile.write( 't # id %s\n' % mol.title.strip() )
            for atom in openbabel.OBMolAtomIter( mol.OBMol):
                label = atom.GetAtomicNum()
                vertex_index = atom.GetIdx()
                args.outfile.write('v %s %s\n' % (vertex_index, label))

            for bond in openbabel.OBMolBondIter( mol.OBMol):
                src_index = bond.GetBeginAtomIdx()
                dest_index = bond.GetEndAtomIdx()
                assert(src_index > 0)
                assert(dest_index > 0)
                if bond.IsAromatic():
                    label = 'a'
                elif bond.IsSingle():
                    label = 's'
                elif bond.IsDouble():
                    label = 'd'
                elif bond.IsTriple():
                    label = 't'
                atom1 = bond.GetBeginAtom()
                atom2 = bond.GetEndAtom()
                args.outfile.write('e %s %s %s\n' % (src_index, dest_index, label))

def run():
	
	inputfile=pybel.readfile(sys.argv[1].split(".")[-1],sys.argv[1])
	value=()
	for mol in inputfile:
		descvalues=mol.calcdesc()
		value= value+(descvalues.get('TPSA'),)
		value= value+(descvalues.get('HBD'),)
		value= value+(descvalues.get('logP'),)
		value= value+(descvalues.get('MW'),)
		value= value+(descvalues.get('tbonds'),)
		value= value+(descvalues.get('nF'),)
		value= value+(descvalues.get('bonds'),)
		value= value+(descvalues.get('atoms'),)
		value= value+(descvalues.get('HBA1'),)
		value= value+(descvalues.get('HBA2'),)
		value= value+(descvalues.get('sbonds'),)
		value= value+(descvalues.get('dbonds'),)
		value= value+(descvalues.get('MR'),)
		value= value+(descvalues.get('abonds'),)
	
		smarts = pybel.Smarts("[+]")
		num=smarts.findall(mol)				
		value= value+(len(num),)			
	
		smarts = pybel.Smarts("[-]")
		num=smarts.findall(mol)				
		value= value+(len(num),)


	model=joblib.load('volume_model/volume.pkl')
	for result in model.predict(value):
		return round(result,2)

def main():
    fa = sys.stdin.read()
    fa = fa.rstrip()
    of = open(args['outfile'], 'w')
    of.write('cid,')
    labels = ''
    for i in properties:
        labels = labels + i + ','
    labels = re.match(r"^(.*),", labels).group(1)
    of.write(labels + '\n')
    inputTemp = tempfile.NamedTemporaryFile(suffix='.sdf', delete=False)
    inputTemp.write(fa)
    inputTempName = inputTemp.name
    inputTemp.close()
    for mol in pybel.readfile('sdf', inputTempName):
        info = ''
        myid = mol.title
        mol.addh()
        mol.make3D()
        desc = mol.calcdesc()
        for thisdesc in properties:
            info = info + str(desc[thisdesc]) + ','
        info = re.match(r"^(.*),", info).group(1)
        of.write(myid.strip() + ',' + info + '\n')
    of.close()
    os.unlink(inputTempName)