Python General.getBase方法代码示例

# 需要导入模块: import General [as 别名]
# 或者: from General import getBase [as 别名]
def alternativeAbundance(matchf, pds, topn = 50):
    MASTER = '/home/anthill/fzheng/home/scripts/termanal_updating'

    # generate the original structures of topn hits
    # pds = General.changeExt(pdb, 'pds')
    # cmd = [MASTER + '/createPDS', '--type', 'query', '--pdb', pdb, '--pds', pds]
    # cmd = ' '.join(cmd)
    # os.system(cmd)
    cmd = [MASTER + '/master', '--query', pds, '--matchIn', matchf, '--structOut', General.getBase(pds) + 'tmp', '--outType', 'match', '--bbRMSD', '--topN', str(topn)]
    cmd = ' '.join(cmd)
    os.system(cmd)

    # for these N structures, calculate RMSD between any two. should be O(N^2)
    odir = os.getcwd()
    ndir = General.getBase(pds) + 'tmp'
    os.chdir(ndir)
    mpdbs = glob.glob('*.pdb')
    mpdbs.sort()
    RMSDs = []
    print 'calculating pairwise RMSD'
    for i in range(len(mpdbs)-1):
        for j in range(i+1, len(mpdbs)):
            mol1, mol2 = parsePDB(mpdbs[i]), parsePDB(mpdbs[j])
            bbAtoms1, bbAtoms2 = mol1.select('backbone').copy(), mol2.select('backbone').copy()
            trans = calcTransformation(bbAtoms2, bbAtoms1)
            bbAtoms2_t = applyTransformation(trans, bbAtoms2)
            rmsd = calcRMSD(bbAtoms1, bbAtoms2_t)
            RMSDs.append(round(rmsd, 3))
    print 'finish calculating RMSD'
    os.chdir(odir)
    # now calculate the average Z-score of all the rmsds of the query
    RMSDs = np.array(RMSDs)
    qRMSD = Analyze.readColumn(matchf, 0, top = topn)
    qRMSD = np.array([float(x) for x in qRMSD])
    meanRMSD, stdRMSD = np.mean(RMSDs), np.std(RMSDs)
    Z_qRMSD = (qRMSD - meanRMSD) / stdRMSD
    return round(np.median(Z_qRMSD), 3)

# 需要导入模块: import General [as 别名]
# 或者: from General import getBase [as 别名]
par.add_argument("--opts", default="", help="other options specifies to run.py")
args = par.parse_args()

selfbin = os.path.dirname(os.path.realpath(sys.argv[0]))
src = selfbin + "/run.py"
pdbs = args.p
pdbs.sort()
jobs = []
odir = os.getcwd()

# homo = '/home/anthill/fzheng/home/termanal_upgrade/caspHomoStrict/'
homo = "/home/anthill/fzheng/ironfs/prepareCASP/CASP10_benchmark/homo/"
for pdb in pdbs:
    caspid = pdb.split("/")[0]
    # caspid = pdb.split('.')[0]
    label = args.d + "/" + caspid + "/" + General.getBase(pdb) + "." + args.label
    if os.path.isfile(label):
        continue
    cmd = ["python", src, "--p", pdb, "--o", args.d + "/" + caspid, "--dontuse", homo + caspid + ".homo", "--hrs", "3"]
    cmd.append(args.opts)
    cmd = " ".join(cmd)
    job = General.jobOnCluster([cmd], os.getcwd(), os.path.realpath(label))
    job.submit(3)
    jobs.append(job)
    sleep(0.5)

while len(jobs) > 0:
    sleep(120)
    for j in jobs:
        j.checkjob()
        if j.running == 0:

# 需要导入模块: import General [as 别名]
# 或者: from General import getBase [as 别名]
    info = l.strip().split('/')
    subdir, name = info[-2], info[-1]
    # use the environment file for the whole protein but cmap file for single chain, don't know if this is good
    envfile = General.changeExt(subdir + '/' + name.split('_')[0], args.e)
    cmapfile = General.changeExt(subdir + '/' + name, args.c)
    
    if not (os.path.isfile(envfile) and os.path.isfile(cmapfile)):
        continue
    
    env = {}
    with open(envfile) as ef:
        f_csv = csv.DictReader(ef, delimiter = '\t')
        for row in f_csv:
            env[row['residue']] = row['environment_score']
    
    cf = open(cmapfile)
    for cfl in cf:
        if not cfl.startswith('contact'):
            continue
        cfa = cfl.strip().split()
        res1, res2, cond, aa1, aa2 = cfa[1:]
        if (not res1 in env) or (not res2 in env):
            continue
        if (float(cond) >= args.range[0]) and (float(cond) <= args.range[1]):
            outstr = '\t'.join(map(str, [General.getBase(name), cond, PDB.t2s(aa1), PDB.t2s(aa2), env[res1], env[res2]])) + '\n'
            out.write(outstr)

         
        
        
    

# 需要导入模块: import General [as 别名]
# 或者: from General import getBase [as 别名]
lines = open(args.tab).readlines()
odir = os.getcwd()
jobs = []
seen = {}
    
for line in lines:
    os.chdir(odir)
    if line.startswith('#'):
        continue
    mut = mustpress.getMutation(line)
    os.chdir(mut.dir)
    pdbs = glob.glob('*.pdb')
    pdbs.sort()

    for pdb in pdbs:
        base = General.getBase(pdb)
        matchf = args.head + '_' + base + '.match'
        if not os.path.isfile(matchf):
            continue
        outname = General.getBase(matchf) + '.' + args.o
        if os.path.isfile(outname):
            continue
        if outname in seen:
            continue
        seen[outname] = 1

        # for the gap test only
        if args.wgap != None:
            pdb = args.wgap + '/' + mut.dir + '/' + pdb
        ##
        pos = PDB.findPositionInPDB(pdb, str(mut.n), mut.c)

# 需要导入模块: import General [as 别名]
# 或者: from General import getBase [as 别名]
import General, PDB

par = argparse.ArgumentParser()
par.add_argument('--l', required = True, help = 'list file')
par.add_argument('--o', required = True, help = 'output file')
par.add_argument('--ext', required = True, help = 'extension of files')
args = par.parse_args()

out = open(args.o, 'w')

for line in open(args.l):
    line = line.strip()
    path = os.path.basename(os.path.dirname(line))
    file = os.path.basename(line)

    name = General.getBase(file)
    cid = name.split('_')[1]
    dbf = path + '/' + name + '.' + args.ext + '.db'
    if not os.path.isfile(db):
        continue
    db = shelve.open(dbf)
    keys = list(db.keys())
    keys_sort = sorted(keys, key = lambda x: int(x.split(',')[1]))
    for k in keys_sort:
        dbr = db[k]
        if (dbr['phi'] > 180.000) or (dbr['phi'] < -180.000) or (dbr['psi'] > 180.000) or (dbr[psi] < -180.000):
            continue
        outstr = [PDB.t2s(dbr['aa']), dbr['sumcond'], dbr['crwdnes'], dbr['freedom'], dbr['phi'], dbr['psi']]
        outstr = '\t'.join(map(str, outstr))
        out.write(outstr +'\n')
    db.close()

# 需要导入模块: import General [as 别名]
# 或者: from General import getBase [as 别名]
# def worker(cmd):
#     sub.call(cmd, shell = True)
    
src = '/home/anthill/fzheng/home/scripts/termanal_updating/updating_scripts/run.py'

pdbs.sort()
jobs = []
odir = os.getcwd()
minhits = 250

# homo = '/home/anthill/fzheng/home/termanal_upgrade/caspHomoStrict/'
for pdb in pdbs:

    ## standard ##
    label = outdir + '/' + General.getBase(pdb) +'.ssc50'
    if os.path.isfile(label):
        continue
    cmd = ['python', src, '--p', pdb, '--o', outdir, '--env', '--minhits', str(minhits)]
    ##

    cmd = ' '.join(cmd)
    job = General.jobOnCluster([cmd], os.getcwd(), os.path.realpath(label))
    job.submit(3)    
    jobs.append(job)
    sleep(0.5)

while (len(jobs) > 0):
    sleep(120)
    for j in jobs:
        j.checkjob()

# 需要导入模块: import General [as 别名]
# 或者: from General import getBase [as 别名]
        cmd_matchInFile = [MASTER + '/master', '--query', pds_file, '--matchIn', matchf1, '--structOut', General.changeExt(seqf1, 'sout'), '--outType', 'match', '--bbRMSD', '--topN', str(topN)]
        cmd_matchInFile = sub.call(cmd_matchInFile)
        match_lines = []
        nmatch = 0
        for match_l in open(matchf1):
            nmatch += 1
            match_lines.append(match_l)
            if nmatch == topN:
                break
        # determine the index of original residue in its own TERM
        rn = findPositionInPDB(frag_pdb, resnum, cid)
        envf = open(General.changeExt(seqf1, 'envi'), 'w')
        for i in range(1, len(match_lines)+1):
            midx = str(i).rjust(len(str(len(match_lines))), '0')
            targetid = General.getBase( General.removePath(match_lines[i-1].split()[1]) )
            env_dict = env_data_path + '/' + targetid[1:3] + '/' + targetid.split('_')[0] + '.env.db'
            db = shelve.open(env_dict, 'r')
            mfrag_pdb = General.changeExt(seqf1, 'sout') + '/match' + midx + '.pdb' 
            mfrag_conres = ConRes(mfrag_pdb)
            mres = mfrag_conres[rn-1]
            resid = mres.getChid() + ',' + str(mres.getResnum()) + res.getIcode()
            if not resid in db:
                envf.write('NA\tNA\n')
            else:
                cond, crwd = db[resid]
                envf.write(cond+'\t'+crwd+'\n')
        shutil.rmtree(General.changeExt(seqf1, 'sout'))

    ematrix = '/home/anthill/fzheng/home/searchDB/statistics/environment_10x10.npy'
    envtable = np.load(ematrix)

# 需要导入模块: import General [as 别名]
# 或者: from General import getBase [as 别名]
if len(sys.argv) -1 != 2:
	print '<usage> [designscore or abundance][a head of files]'
	exit(0)
score, head = sys.argv[1:]

dirs = [x for x in os.listdir('.') if os.path.isdir(x)]
dirs.sort()
for d in dirs:
	outf = d + '/' + score + '.' + head + '.info'
	outfh = open(outf, 'w')
	outfh.write('\t'.join(['residue_number', 'number_of_hits', 'number_of_contacts', 'number_of_residues'])+'\n')
	lists  = glob.glob(d + '/fragments/*.list')
	lists.sort()
	info = []
	for lst in lists:
		base = General.removePath(General.getBase(lst))
		resnum = base.split('_')[-1][1:]
		ncon = len(open(lst).readlines())
		pdb = d + '/fragments/' + base + '.pdb'
		nres = len(PDB.ConRes(pdb))
		seqf = d + '/' + score + '/' + head + '_' + base + '.seq'
		nhit = len(open(seqf).readlines())
		record = [resnum, nhit, ncon, nres]
		info.append(record)
	info_sort = sorted(info, key = lambda x :int(x[0]))
	for i_s in info_sort:
		outstr = '\t'.join(map(str, i_s)) + '\n'
		outfh.write(outstr)
	outfh.close()

# 需要导入模块: import General [as 别名]
# 或者: from General import getBase [as 别名]
    hits_all = []

    local_envf = args.local + '/' + d + '/' + args.localh + '_' + d + '.' + args.env[1]
    hits_local = []
    if os.path.isfile(local_envf):
        for l in open(local_envf).readlines():
            hit = ' '.join(l.strip().split()[0:2])
            hits_local.append(hit)

    for i in range(0, len(envfs)):
        hits = []
        for l in open(envfs[i]).readlines():
            hit = ' '.join(l.strip().split()[0:2])
            hits.append(hit)
        hits_all.append(hits)

    for i in range(0, len(envfs)):
        outstr = [General.getBase(envfs[i]).replace(args.head + '_', ''), ':']
        overlapToLocal = 0
        if len(hits_local) > 0:
            overlapToLocal = calculateOverlap(hits_all[i], hits_local)
        outstr.extend(['local', format(overlapToLocal, '.2f')])
        for j in range(0, len(envfs)):
            if i == j:
                continue
            overlap = calculateOverlap(hits_all[i], hits_all[j])
            outstr.extend([General.getBase(envfs[j]).replace(args.head + '_', ''), format(overlap, '.2f')])
        out.write('\t'.join(outstr) + '\n')

    out.close()

# 需要导入模块: import General [as 别名]
# 或者: from General import getBase [as 别名]
import os, sys
import General

inp = sys.argv[1]

# read fasta of casp files
fa = '/home/anthill/fzheng/home/termanal_upgrade/caspModels.fa'
fals = open(fa).readlines()
length = {}

db1 = {}
db1list = '/home/anthill/fzheng/ironfs/designScore/bc-30-sc-20140815/list'
for l in open(db1list):
	pid = General.getBase( General.removePath(l) )
	db1[pid] = 1

db2 = {}
db2list = '/home/anthill/fzheng/home/searchDB/bc-30-sc-20141022-newpds/list'
for l in open(db2list):
	pid = General.getBase( General.removePath(l) )
	db2[pid] = 1

for i in range(len(fals)):
	if fals[i].startswith('>'):
		qid = fals[i][1:].split('.')[0]
		length[qid] = len(fals[i+1].strip())

for inpl in open(inp):
	items = inpl.strip().split()
	qid = items[0].split('.')[0]
	items[0] = qid

# 需要导入模块: import General [as 别名]
# 或者: from General import getBase [as 别名]
import sys, glob
import General

if len(sys.argv) - 1 != 3:
	print '<usage> [protein] [extension] [output file]'
	exit(0)
pro, ext, outf = sys.argv[1:]

files = glob.glob(pro + '*' + ext)
files.sort()
scores = []
for f in files:
	score, n = 0.0, 0
	lines = open(f).readlines()
	for l in lines:
		score += float(l.strip().split()[-1])
		n += 1
	score /= n
	scores.append([General.getBase(f), score])
# scores = sorted(scores, reverse = True, key = lambda x : x[1])
 
out = open(outf, 'w')
for s in scores:
	out.write(s[0] + '\t'+format(s[1], '.3f')+'\n')

# 需要导入模块: import General [as 别名]
# 或者: from General import getBase [as 别名]
dirs = [x for x in os.listdir('.') if os.path.isdir(x)]
dirs.sort()

for d in dirs:
	if args.cont:
		confiles = [d + '/' + x for x in os.listdir(d) if x.endswith(args.ext)]
	else:
		confiles = [d +'/' + [x for x in os.listdir(d) if x.endswith(args.ext)][0] ]# side chain dependent info is not used here
	pid = d.split('_')[0].lower()
	pdb = args.sdir + '/' + pid + '.clean.pdb'

	for confile in confiles:
		conlst = open(confile)
		for con in conlst:
			info = con.strip().split()

			idx, res1, res2, cond = info[0], info[2], info[3], float(info[6])
			if args.cont:
				cond = max(float(info[4]), float(info[5]))
			seeds = []

			if cond > args.cut:
				seeds.extend([res1, res2])
			else:
				continue
			outpdb = General.getBase(confile)[:-1] + '_' + res2.replace(',', '') + '.pdb'
			if os.path.isfile(outpdb):
				continue
			fragmade = Fragment.makeFragments(pdb, seeds, outFile = outpdb, flank = args.flank)
			if fragmade == -1:
				print outpdb, 'not successfully created'

# 需要导入模块: import General [as 别名]
# 或者: from General import getBase [as 别名]
selfbin = General.selfbin(sys.argv[0])

par = argparse.ArgumentParser()
par.add_argument('--l', required = True, help = 'a list file')
par.add_argument('--ext', required = True, help = 'the extension of output file')
args = par.parse_args()

batchsize = 30 # will run 30 command in a single job

count = 0
pdbfs = []

for l in open(args.l):
    info = l.strip().split('/')
    subdir = info[-2]
    name = General.getBase(info[-1])

    pdbf = subdir + '/' + name + '.pdb'
    outf = subdir + '/' + name + '.' + args.ext
    if (not os.path.isfile(pdbf)) or (os.path.isfile(outf)):
        continue

    pdbfs.append(pdbf)
    count += 1

    if count == 30:
        cmd = ['python', selfbin+'/extractSCcontacts.py', '--ext', args.ext, '--p'] + pdbfs
        cmd = ' '.join(cmd)
        General.qsub([cmd])
        pdbfs = []
        count = 0

# 需要导入模块: import General [as 别名]
# 或者: from General import getBase [as 别名]
# output file names
nr_matchf = args.outh + '_' + args.m
nr_seqf = General.changeExt(nr_matchf, 'seq')
nr_env = None
oenv = General.changeExt(args.m, args.env)
if os.path.isfile(oenv):
	nr_env = General.changeExt(nr_matchf, args.env)


# write a custom .fasta file
for match in matches:
	match_region_indices = Analyze.index_from_match(match)
	central_index = match_region_indices[args.cres - 1]

	match_id = General.getBase( General.removePath( match.split()[1] ) )
	fullsequence = database[match_id]
	if central_index - args.wd < 1:
		seqcontext = fullsequence[0:(2 * args.wd + 1)]
	elif central_index + args.wd > len(fullsequence):
		seqcontext = fullsequence[-(2 * args.wd + 1):]
	else:
		seqcontext = fullsequence[(central_index - args.wd - 1):(central_index + args.wd)]
	tempfh.write('>match:'+str(matchind)+'\n'+seqcontext+'\n')

	if args.conres != None:
		con_index = match_region_indices[args.conres -1]
		if con_index - args.wd < 1:
			seqcontext = fullsequence[0:(2 * args.wd + 1)]
		elif central_index + args.wd > len(fullsequence):
			seqcontext = fullsequence[-(2 * args.wd + 1):]

# 需要导入模块: import General [as 别名]
# 或者: from General import getBase [as 别名]
	os.chdir(odir)
	os.chdir(d)
	pdbs = glob.glob('*.pdb')
	cmds = []
	resn = int(d.split('_')[1][2:])
	for pdb in pdbs:
		matchf = args.head + '_' + General.changeExt(pdb, 'match')
		if not os.path.isfile(matchf):
			continue

		pos = PDB.findPositionInPDB(pdb, resn)
		# if output file is already there, skip the job
		if os.path.isfile('nr'+args.id +'_'+matchf):
			continue

		cmd = ['python', selfbin + '/removeLocalRedundancy.py', '--m', matchf, '--cres', str(pos), '--id', args.id, '--outh', 'nr'+args.id]
		if not args.db == None:
			cmd.extend(['--db', args.db])
		if args.conR:
			conresn = General.getBase(pdb).split('_')[2][1:]
			conpos = PDB.findPositionInPDB(pdb, conresn)
			cmd.extend(['--conres', str(conpos)])
		if args.env != None:
			cmd.extend(['--env', args.env])
		cmd = ' '.join(cmd)
		cmds.append(cmd)

	if len(cmds) > 0:
		General.qsub(cmds)