Python Trajectory.load_from_pdb方法代码示例

# 需要导入模块: from msmbuilder import Trajectory [as 别名]
# 或者: from msmbuilder.Trajectory import load_from_pdb [as 别名]
def main(modeldir, genfile,  type, write=False):
    data=dict()
    pops=numpy.loadtxt('%s/Populations.dat' % modeldir)
    map=numpy.loadtxt('%s/Mapping.dat' % modeldir)
    frames=numpy.where(map!=-1)[0]

    unbound=numpy.loadtxt('%s/tpt-rmsd-%s/unbound_%s_states.txt' % (modeldir, type, type), dtype=int)
    bound=numpy.loadtxt('%s/tpt-rmsd-%s/bound_%s_states.txt' % (modeldir, type, type), dtype=int)

    dir=modeldir.split('Data')[0]
    name=glob.glob('%s/fkbp*xtal*pdb' % dir)
    pdb=Trajectory.load_from_pdb(name[0])
    paths=io.loadh('%s/tpt-rmsd-%s/Paths.h5' % (modeldir, type))

    committors=numpy.loadtxt('%s/commitor_states.txt' % modeldir, dtype=int)
    colors=['red', 'orange', 'green', 'cyan', 'blue', 'purple']
    colors=colors*40
    if type=='strict':
        ref=5
    elif type=='super-strict':
        ref=3
    elif type=='medium':
        ref=10
    elif type=='loose':
        ref=15
    #for p in range(0, 3):
    for p in range(0, 1):
        path=paths['Paths'][p]
        print "Bottleneck", paths['Bottlenecks'][p]
        flux=paths['fluxes'][p]/paths['fluxes'][0]
        if flux < 0.2:
            break
        print "flux %s" % flux
        frames=numpy.where(path!=-1)[0]
        path=numpy.array(path[frames], dtype=int)
        print path
        if write==True:
            size=(paths['fluxes'][p]/paths['fluxes'][0])*1000
            traj=Trajectory.load_from_xtc('%s/tpt-rmsd-%s/path%s_sample20.xtc' % (modeldir, type, p), Conf=pdb)
            data=build_metric(dir, pdb, traj)
            dir=modeldir.split('Data')[0]
            for op in sorted(data.keys()):
            #for op in residues:
                pylab.figure()
                pylab.scatter(data['rmsd'], data[op], c=colors[p], alpha=0.7) #, s=size)
                for j in paths['Bottlenecks'][p]:
                    frame=numpy.where(paths['Paths'][p]==j)[0]
                    pylab.scatter(data['rmsd'][frame*20], data[op][frame*20], marker='x', c='k', alpha=0.7, s=50)
                    location=numpy.where(committors==paths['Paths'][p][frame])[0]
                    if location.size:
                        print "path %s state %s bottleneck in committors" % (p, j)
                        print data['rmsd'][frame*20], data[op][frame*20]
                pylab.title('path %s' % p)
                pylab.xlabel('P-L RMSD')
                #pylab.xlabel('P-L COM')
                pylab.ylabel(op)
                pylab.xlim(0,max(data['rmsd'])+5)
                #pylab.ylim(0,max(data[op])+5)
                pylab.show()

# 需要导入模块: from msmbuilder import Trajectory [as 别名]
# 或者: from msmbuilder.Trajectory import load_from_pdb [as 别名]
def test_asa_3():

    traj_ref = np.loadtxt( os.path.join(reference_dir(),'g_sas_ref.dat'))
    Conf = Trajectory.load_from_pdb(os.path.join( fixtures_dir(), 'native.pdb'))

    traj = Trajectory.load_trajectory_file( os.path.join(fixtures_dir(), 'trj0.xtc') , Conf=Conf)
    traj_asa = calculate_asa(traj, probe_radius=0.14, n_sphere_points = 960)
    
    # the algorithm used by gromacs' g_sas is slightly different than the one
    # used here, so the results are not exactly the same -- see the comments
    # in src/python/geomtry/asa.py or the readme file src/ext/asa/README.txt
    # for details
    npt.assert_array_almost_equal(traj_asa, traj_ref, decimal=2)    

# 需要导入模块: from msmbuilder import Trajectory [as 别名]
# 或者: from msmbuilder.Trajectory import load_from_pdb [as 别名]
logger.setLevel(logging.INFO)
sh = logging.StreamHandler(sys.stdout)
formatter = logging.Formatter(fmt="%(asctime)s - %(message)s", datefmt="%H:%M:%S")
sh.setFormatter(formatter)
logger.addHandler(sh)
logger.propagate = False

logger.info("start")
Proj = Project.load_from(args.proj_FN)
logger.info("loaded project info")
try:
    Ass = io.loadh(args.ass_FN)["arr_0"]
except:
    Ass = io.loadh(args.ass_FN)["Data"]

pdb = Trajectory.load_from_pdb(Proj.conf_filename)

which = np.loadtxt(args.which).astype(int)

distance_cutoff = 0.32
angle_cutoff = 120


def get_hb(traj):

    # get accH - donor distance:
    dists = contact.atom_distances(traj["XYZList"], atom_contacts=which[:, 1:])

    # get angles
    angles = angle.compute_angles(traj, which)

# 需要导入模块: from msmbuilder import Trajectory [as 别名]
# 或者: from msmbuilder.Trajectory import load_from_pdb [as 别名]
def main(system, genfile, lag, volume, nproc):
    lag=int(lag)
    dir=os.path.dirname(genfile)
    if not os.path.exists('%s/target.txt' % dir.split('Data')[0]):
        print "need target.txt with experimental value"
        sys.exit()
    ref=loadtxt('%s/target.txt' % dir.split('Data')[0])
    filename=genfile.split(dir)[1].split('.lh5')[0]
    if "Coarse" in filename:
        coarse=filename.split('Coarsed')[1].split('_')[0]
        rcut=filename.split('_r')[1].split('_')[0]
        modeldir='%s/msml%s_coarse_r%s_d%s/' % ( dir.split('/')[0], lag, rcut, coarse)
    else:
        modeldir='%s/msml%s' % (dir, lag)
    print "computing PMF from %s with lag %s frames" % (filename, lag)
    conf=Trajectory.load_from_pdb('%s/%s_noh.pdb' % (dir.split('Data')[0], system))
    lig_atoms=conf['XYZList'].shape[1]
    map=loadtxt('%s/Mapping.dat' % modeldir)
    pops=loadtxt('%s/Populations.dat' % modeldir)

    ligfile='%s.vmd_ligcoords.pickle' % genfile.split('.lh5')[0]
    if os.path.exists(ligfile):
        lighandle=open(ligfile, 'rb')
    else:
        print "need to run tcl scripts and pickler for coordinates"
        sys.exit()

    ligcoors=pickle.load(lighandle)
    lighandle.close()
    mapped_ligcoors, x_range, y_range, z_range, box_volume=get_minmax(ligcoors, map)
    correction=-0.6*log(box_volume/1600.0)
    # get prot-lig distances
    com_distances=loadtxt('%s.vmd_com.dat' % genfile.split('.lh5')[0], usecols=(1,))
    ref_com=com_distances[0]
    com_distances=com_distances[1:]
    if len(com_distances)==len(map):
        print "protein-ligand COM distances per state exist"
    # PMF
    frames=where(map!=-1)[0]
    mapped_states=map[frames]
    mapped_com=range(0, len(frames))
    mapped_com_distances=com_distances[frames]
    space=PMF3D.PMF3D(x_range, y_range, z_range)
    spacetrack=space.microstate_count_grid(mapped_ligcoors)
    new_pops={ key: pops[n] for (n, key) in enumerate(mapped_ligcoors.keys())}
    GD=space.new_manual_allcoor_grid(spacetrack, mapped_ligcoors, new_pops, type='pops')
    free=array([-0.6*log(i) for i in pops])
    subtract=min(free)
    free=array([k-subtract for k in free])
    GD=space.new_manual_allcoor_grid(spacetrack, mapped_ligcoors, new_pops, type='pops')
    GDfree=-0.6*log(GD)
    GDfree=PMF3D.convert(GDfree, max(free))
    GDfree=GDfree-min(GDfree.flatten())
    space.write_dx(GDfree, modeldir)
    frees=[]
    corrs=[]
    axis=[]
    volumes=[]
    cutoffs=arange(0, 40, 1)
    if volume=='all':
        for cutoff in cutoffs:
            bound_frames=where(mapped_com_distances < cutoff)[0]
            if len(bound_frames)==0:
                print "no bound states less than reference com distance %s" % cutoff
                continue
            gen_bound_frames=[]
            for i in bound_frames:
                location=where(map==i)[0]
                gen_bound_frames.append(location)
            new_points={ key: mapped_ligcoors[key] for key in bound_frames}
            new_pops={ key: pops[key] for key in bound_frames}
            GD=space.new_manual_allcoor_grid(spacetrack, new_points, new_pops, type='pops')
            boundspace=space.pmfvolume(GD)

            new_pops={ key: 1.0 for key in bound_frames}
            GD=space.new_manual_allcoor_grid(spacetrack, new_points, new_pops, type='pops')
            boundvolume=space.pmfvolume(GD)
            volumes.append(boundvolume)
            print "count bound volume ", boundvolume

            # unbound frames are above COM cutoff
            unbound_frames=array([int(x) for x in mapped_states if x not in bound_frames])
            new_points={ key: mapped_ligcoors[key] for key in unbound_frames}
            new_pops={ key: pops[key] for key in unbound_frames}
            GD=space.new_manual_allcoor_grid(spacetrack, new_points, new_pops, type='pops')
            unboundspace=space.pmfvolume(GD)

            # for counting states
            new_pops={ key: 1.0 for key in unbound_frames}
            GD=space.new_manual_allcoor_grid(spacetrack, new_points, new_pops, type='pops')
            unboundvolume=space.pmfvolume(GD)
            #print "unbound volume ", unboundvolume

            # free energy from ratio
            depth=-0.6*log(boundspace/unboundspace)
            frees.append(depth)
            axis.append(cutoff)
            print "corrected integrated dG ratio at cutoff %s is %s" % (cutoff, depth+correction)
        k=len(pops)
        savetxt('%s/%s_k%s_l%s_pocketdistance_free.dat' % (modeldir, filename, k, lag), [x+correction for x in frees])
#.........这里部分代码省略.........

# 需要导入模块: from msmbuilder import Trajectory [as 别名]
# 或者: from msmbuilder.Trajectory import load_from_pdb [as 别名]
formatter = logging.Formatter(fmt="%(asctime)s - %(message)s", datefmt="%H:%M:%S")
sh.setFormatter(formatter)
logger.addHandler(sh)
logger.propagate = False

logger.info("start")
Proj = Project.load_from(args.proj_FN)
logger.info("loaded project info")
try:
    Ass = io.loadh(args.ass_FN)["arr_0"]
except:
    Ass = io.loadh(args.ass_FN)["Data"]

HB = metric_HB.HydrogenBond()

pdb = Trajectory.load_from_pdb(Proj.conf_filename)

donorH_ainds = np.where(((pdb["AtomNames"] == "H") | (pdb["AtomNames"] == "HN")) & (pdb["ResidueNames"] != "PRO"))[
    0
]  # _ainds correspond to atom indices in the trajectory
donor_ainds = np.where((pdb["AtomNames"] == "N") & (pdb["ResidueNames"] != "PRO"))[0]
acceptor_ainds = np.where(pdb["AtomNames"] == "O")[0]  # This excludes the C-Terminus which has naming issues...


n_res = np.unique(pdb["ResidueID"]).shape[0]
ppdb = HB.prepare_trajectory(pdb)  # Do this to set the donor indices and stuff

atom_indices = np.sort(np.concatenate((donor_ainds, donorH_ainds, acceptor_ainds)))

CMs_1d = np.zeros((Ass.max() + 1, ppdb.shape[1]))

# 需要导入模块: from msmbuilder import Trajectory [as 别名]
# 或者: from msmbuilder.Trajectory import load_from_pdb [as 别名]
residue_pairs=numpy.loadtxt('./residuepairs_0based_select_structure1.dat', dtype=int)
pairmetric2=metrics.ContinuousContact(contacts=residue_pairs, scheme='closest-heavy')
#pairmetric min 0.00255089723477 max 10.8130541983

residue_pairs=numpy.loadtxt('./residuepairs_select_0based.dat', dtype=int)
pairmetric1=metrics.ContinuousContact(contacts=residue_pairs, scheme='closest-heavy')
#pairmetric1 min 0.00447584821743 max 22.0643260692 

#pairmetric=metrics.BooleanContact(metric='matching', contacts=residue_pairs, cutoff=0.5, scheme='closest-heavy')
ligand_inds = numpy.loadtxt('./p53_Calpha_indices.dat', dtype='int')
prot_inds = numpy.loadtxt('./s100b_Calpha_indices.dat', dtype='int')

rmsdmetric2=metrics.RMSD(ligand_inds)
rmsdmetric1 = lprmsd.LPRMSD(atomindices=prot_inds, permuteindices=None, altindices=ligand_inds)
pdb=Trajectory.load_from_pdb('structure_1.pdb')
#project = Project.load_from('/home/kkappel/p53/s100b_md/msm_80ns_p53CA/ProjectInfo.yaml')
project = Project.load_from('ProjectInfo.yaml')

hybrid1 = metrics.Hybrid([rmsdmetric2, pairmetric1], weights=[1/1.16, 1/22.1])
hybrid2 = metrics.Hybrid([rmsdmetric2, pairmetric2], weights=[1/1.16, 1/10.8])


data=dict()
# check pair metric for traj0
metrics=[pairmetric1, pairmetric2, hybrid1, hybrid2 ]
names=['pair-selectall', 'pair-selectstructure1', 'hybrid-selectall', 'hybrid-selectstructure1']
for (metric, name) in zip(metrics, names):
    data[name]=dict()
    data[name]['min']=1000000
    data[name]['max']=0

# 需要导入模块: from msmbuilder import Trajectory [as 别名]
# 或者: from msmbuilder.Trajectory import load_from_pdb [as 别名]
from pyschwancr import dataIO
import os, sys, re
 
matplotlib.rcParams['font.size']=22
pp = PdfPages( args.out_plot )

HB = metric_HB.HydrogenBond()

input_dict = io.loadh( args.in_cm )
triples = input_dict['donor_h_acceptor_ainds']
CMs = input_dict['HB_maps']

num_acceptors = len( np.unique( triples[:,2] ) )
num_donors = len( np.unique( triples[:,0] ) )

nat_pdb = Trajectory.load_from_pdb( args.nat_FN )
T = mmread( args.tProb )
n_res = np.unique( nat_pdb['ResidueID'] ).shape[0]
CM_pdb = HB.prepare_trajectory( nat_pdb )

CM_pdb = CM_pdb.reshape( (1,num_donors, num_acceptors), order='C')[0]

native_locs = np.array(np.where( CM_pdb )).T

vals,vecs = MSMLib.GetEigenvectors( T, args.num_vecs+1 )
vecs = vecs.real
    
if args.lagtime != None:
   timescales = - args.lagtime / np.log( vals[1:].real ) # NOTE: This means the 1st dyn e-vec is indexed at 0
  
for i in xrange(1,args.num_vecs+1):

# 需要导入模块: from msmbuilder import Trajectory [as 别名]
# 或者: from msmbuilder.Trajectory import load_from_pdb [as 别名]
def construct_basic_metric(metric_name, args):
    if metric_name == 'rmsd':
        if args.rmsd_atom_indices != 'all':
            atom_indices = np.loadtxt(args.rmsd_atom_indices, np.int)
        else:
            atom_indices = None
        metric = RMSD(atom_indices)#, omp_parallel=args.rmsd_omp_parallel)

    elif metric_name == 'dihedral':
        metric = Dihedral(metric=args.dihedral_metric,
            p=args.dihedral_p, angles=args.dihedral_angles,
            userfilename=args.dihedral_userfilename)
    
    elif metric_name == 'contact':
        if args.contact_which != 'all':
            contact_which = np.loadtxt(args.contact_which, np.int)
        else:
            contact_which = 'all'

        if args.contact_cutoff_file != None: 
            contact_cutoff = np.loadtxt(args.contact_cutoff_file, np.float)            
        elif args.contact_cutoff != None:
            contact_cutoff = float( args.contact_cutoff )
        else:
            contact_cutoff = None
             
        if contact_cutoff != None and contact_cutoff < 0:
            metric = ContinuousContact(contacts=contact_which,
                scheme=args.contact_scheme)
        else:
            metric = BooleanContact(contacts=contact_which,
                cutoff=contact_cutoff, scheme=args.contact_scheme)
     
    elif metric_name == 'atompairs':
        if args.atompairs_which != None:
            pairs = np.loadtxt(args.atompairs_which, np.int)
        else:
            pairs = None

        metric = AtomPairs(metric=args.atompairs_metric, p=args.atompairs_p,
            atom_pairs=pairs)

    elif metric_name == 'positions':
        target = Trajectory.load_from_pdb(args.target)
        
        if args.pos_atom_indices != None:
            atom_indices = np.loadtxt(args.pos_atom_indices, np.int)
        else:
            atom_indices = None

        if args.align_indices != None:
            align_indices = np.loadtxt(args.align_indices, np.int)
        else:
            align_indices = None
        
        metric = Positions(target, atom_indices=atom_indices, align_indices=align_indices,
                           metric=args.positions_metric, p=args.positions_p)
             
    elif metric_name == 'custom':
        with open(args.picklemetric_input) as f:
            metric = pickle.load(f)
            print '#'*80
            print 'Loaded custom metric:'
            print metric
            print '#'*80
    else:
        # apply the constructor on args and take the first non-none element
        # note that using these itertools constructs, we'll only actual
        # execute the constructor until the match is achieved
        metrics = itertools.imap(lambda c: c(args), locate_metric_plugins('construct_metric'))
        try:
            metric = itertools.dropwhile(lambda c: not c, metrics).next()
        except StopIteration:
            # This means that none of the plugins acceptedthe metric
            raise RuntimeError("Bad metric. Could not be constructed by any built-in or plugin metric. Perhaps you have a poorly written plugin?")
     
    if not isinstance(metric, AbstractDistanceMetric):
        return ValueError("%s is not a AbstractDistanceMetric" % metric)

    return metric

# 需要导入模块: from msmbuilder import Trajectory [as 别名]
# 或者: from msmbuilder.Trajectory import load_from_pdb [as 别名]
    atomweights['N']=14
    atomweights['H']=1
    atomweights['S']=32
    weights=[]
    for x in names:
        try:
            int(x[0])
            weights.append(atomweights[x[1]])
        except ValueError:
            weights.append(atomweights[x[0]])
    return numpy.array(weights)


# load in trajectory, here I have just a pdb with one frame
# you can loop over frames in a trajectory
conf=Trajectory.load_from_pdb('new-aln-state847.cent.pdb')

# here i organize the groups in dictionaries
coms=dict()
groups=dict()
group_wts=dict()
residues=dict()
numbers=[1,2,3,4]
start=1
for num in numbers:
    groups[num]=dict() # group ids
    group_wts[num]=dict() # group mass weight ids
    end=start+1
    residues[num]=[start,end] # adding consecutive resids to groups in this example
    start=end+1

# 需要导入模块: from msmbuilder import Trajectory [as 别名]
# 或者: from msmbuilder.Trajectory import load_from_pdb [as 别名]
from fpop import FPOPExchanger
import numpy as np
from scipy.io import mmread
from msmbuilder import Trajectory
import IPython
from matplotlib.pyplot import *
ion()

t = mmread('tProb.mtx')

pdb = Trajectory.load_from_pdb('proteinG_wt.rename.pdb')

avg_sasas = np.loadtxt('avg_sasas.dat')

timepoints = np.arange(0, 50, 2)

print timepoints.shape

exchanger = FPOPExchanger(t, avg_sasas, 100E-6, 1E-3, 20E-3, pdb, timepoints, lagtime=50E-9, init_pops=np.loadtxt('init48.dat'), force_dense=False)

total_products, total_res_pops = exchanger.run(return_res_pops=True, num_labels=10)

for i in range(5):
    plot(timepoints * 0.05, total_products[:, i], label='%d labels' % i, lw=3)

legend(prop={'size' : 14})
xlabel(r'time ($\mu s$)')
ylabel('Population')

IPython.embed()

# 需要导入模块: from msmbuilder import Trajectory [as 别名]
# 或者: from msmbuilder.Trajectory import load_from_pdb [as 别名]
from msmbuilder.geometry import asa
from msmbuilder import Trajectory
from msmbuilder import io
import numpy as np

p=Trajectory.load_from_pdb('barstarH.pdb')

sasa = asa.calculate_asa(p, n_sphere_points=150)

io.saveh('sasa.h5', sasa=sasa)