gen220project/treeanalysis.py at master · erg55/gen220project · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
from Bio import Phylo
import os
import sys
import glob
import getpass
#filepath input
if len(sys.argv) == 3:
	files = glob.glob(sys.argv[1]+"/*.phylip")
	#files = ['./original_data/T58_L1.phylip']
	bootreps = " -N "+sys.argv[2]
else:
	print "FORMAT: python treeanalysis.py [folder with phylip] [bootstrap replicates]"
	print "EXAMPLE: python treeanalysis.py ./phylip 1000"
	sys.exit()
if len(files) == 0:
	print "no phylip files in the directory"
#files = ['./original_data/T58_L1.phylip']
#starting to process files
progbarc = 0
result = {}
#test
for f in files:
	print f
 	raxmlf = " -s "+f
 	#check for raxml leftover files
 	testcheck = glob.glob("./*.TEST")
 	for test in testcheck:
 		if test.split("/")[1] in os.listdir("./"):
			os.remove(test)
	#run RAxML
	user = getpass.getuser()
	if user == "aknys001":
		execfile('/usr/share/Modules/init/python.py')
		module('load','RAxML/8.2.3')
		os.system("cd ~/gen220project")
		os.system("raxmlHPC-PTHREADS-SSE3 -T 16 --silent -f a -c 25 -p 12345 -x 12345 -m GTRCAT -n TEST"+raxmlf+bootreps)
	else:
		os.system("raxmlHPC --silent -f a -c 25 -p 12345 -x 12345 -m GTRCAT -n TEST"+raxmlf+bootreps)
	tree = Phylo.read("RAxML_bipartitions.TEST", "newick")
	#tree = Phylo.read("./../testtree.tre", "newick")
	if tree.is_bifurcating():
		print "bifurcating"
	else:
		print "not bifurcating"
	counter = 0
	ct = 0
	test = tree.as_phyloxml()
	list4 = test.find_clades()
	totconf = 0
	lowconf = 0
	for l1 in list4:
		counter +=1
		if l1.is_bifurcating():
			ct +=1
		l2 = str(l1._get_confidence())
		if l2.find("Confidence") == 0:
			totconf +=1
			l3=l2.split(",")
			l4 = l3[1][:-1].split("=")
			if int(l4[1]) <70:
				lowconf +=1
	print "below 70 per cent", float(lowconf)/(totconf)
	print "clades #", counter, ct
	fname = f.split("/")[2]
	tr = fname.split(".")[0]
	result[tr] = float(lowconf)/(totconf)
	#clean up
	if user == "aknys001":
		os.remove("RAxML_bestTree.TEST")
		os.remove("RAxML_info.TEST")
		os.remove("RAxML_bipartitions.TEST")
 		os.remove("RAxML_bipartitionsBranchLabels.TEST")
 		os.remove("RAxML_bootstrap.TEST")
	#progress bar
	progbarc +=1
	progbar = int(round(float(progbarc)/len(files)*100, 0))
	hashes = '#' * int(progbar * 0.2)
	spaces = ' ' * (20 - len(hashes))
	print "-------------------------------------"
	print "\rProgress: [{0}] {1}%".format(hashes + spaces, progbar)
	print "-------------------------------------"
#output the final table
with open("treesfilter.tab", "w") as outfile:
	for tc, tcv in result.items():
		print >> outfile, tc, "\t", tcv
print "Done"