-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathmake_run_thinned.py
70 lines (65 loc) · 3.25 KB
/
make_run_thinned.py
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
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Mon Jul 1 15:29:33 2019
@author: lpsmith
"""
rootdir = "tsinfer_1xr/"
reactionRate = "2.679985"
ntips = 1000
x2tips = ntips*2
runall = open("run_all.bat", "w")
for i in ["0"]:
for j in ["0", "1", "2", "3", "4", "5", "6", "7", "8", "9"]:
for k in ["0", "1", "2", "3", "4", "5", "6", "7", "8", "9"]:
ijk = i+j+k
simout = open("sim_thinned_" + ijk + ".sge", "w")
runall.write("qsub -l mfree=4G sim_thinned_" + ijk + ".sge\n")
simout.write("module load modules modules-init modules-gs gcc/latest python/latest numpy/1.16.0 scipy/1.2.0 glibc/2.14\n")
simout.write("cd " + rootdir + str(ntips) + "tip_thinned/\n")
simout.write("rm reccount\n")
simout.write("rm migcount\n")
simout.write("rm totreccount\n")
simout.write("rm RFresults\n")
simout.write("# ./ms " + str(x2tips) + " 1 -T -t 550.000000 -r " + reactionRate + " 100000 | tail -n +4 | grep -v // > truetree" + ijk + "\n")
simout.write("python3 truncate_ms_tree.py truetree" + ijk + " treefile" + ijk + "\n")
simout.write("python thinned.py treefile" + ijk + " thinnedtree" + ijk + " " + str(ntips) + "\n")
simout.write("./seq-gen -mF84 -l 100000 -s 0.005500 -p 1000 thinnedtree" + ijk + " > dnafile" + ijk + "\n")
simout.write("python3 dna_to_bool.py dnafile" + ijk + " translated_seqgen_dna" + ijk + "\n")
simout.write("python3 parse_seqgen.py translated_seqgen_dna" + ijk + " treedump" + ijk + " inferredtree" + ijk + "\n")
#simout.write("mv timefile timefile" + ijk + "\n")
simout.write("# compare converted original tree to inferred tree\n")
simout.write("python ts_argcompare.py treefile" + ijk + " inferredtree" + ijk + "\n")
simout.write("rm seqgen_dna translated_seqgen_dna treedump" + ijk + "\n")
simout.write("rm reccount migcount totreccount\n")
simout.close()
runall.close()
initout = open("sim_init_trees.sge", "w")
initout.write("module load modules modules-init modules-gs gcc/latest python/latest numpy/1.16.0 scipy/1.2.0 glibc/2.14\n")
initout.write("cd " + rootdir + str(ntips) + "tip_thinned/\n")
initout.write("for i in 0\n")
initout.write("do\n")
initout.write("for j in 0 1 2 3 4 5 6 7 8 9\n")
initout.write("do\n")
initout.write("for k in 0 1 2 3 4 5 6 7 8 9\n")
initout.write("do\n")
initout.write(" ./ms " + str(x2tips) + " 1 -T -t 550.000000 -r " + reactionRate + " 100000 | tail -n +4 | grep -v // > truetree$i$j$k\n")
initout.write("done\n")
initout.write("done\n")
initout.write("done\n")
initout.close()
argout = open("sim_argcompare.sge", "w")
argout.write("module load modules modules-init modules-gs gcc/latest python/latest numpy/1.16.0 scipy/1.2.0 glibc/2.14\n")
argout.write("cd " + rootdir + str(ntips) + "tip_thinned/\n")
argout.write("rm RFresults\n")
argout.write("for i in 0\n")
argout.write("do\n")
argout.write("for j in 0 1 2 3 4 5 6 7 8 9\n")
argout.write("do\n")
argout.write("for k in 0 1 2 3 4 5 6 7 8 9\n")
argout.write("do\n")
argout.write(" python ts_argcompare.py thinnedtree$i$j$k inferredtree$i$j$k\n")
argout.write("done\n")
argout.write("done\n")
argout.write("done\n")
argout.close()