EvolutionModels/create_branchandsite.py at master · leanngward/EvolutionModels · 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
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
import os
import sys
def null_codeml_branches_runner(codeml_prefix,seq_file,treefile,outfile,scriptfile_b,gene_name):
    import subprocess
    #create the control file for codeml to run
    control_file_name = codeml_prefix  + gene_name + '_null.ctl'
    with open(control_file_name, 'w') as fh:
        fh.write('seqfile = ' + seq_file +'\n')
        fh.write('treefile = '+ treefile +'\n')
        fh.write('outfile = ' + outfile +'\n')
        fh.write('noisy = 9 \n')
        fh.write('verbose = 0 \n')
        fh.write('runmode = 0 \n')
        fh.write('seqtype = 1 \n')
        fh.write('CodonFreq = 2 \n')
        fh.write('clock = 0 \n')
        fh.write('aaDist = 0 \n')
        fh.write('model = 2 \n')
        fh.write('NSsites = 2 \n')
        fh.write('Mgene = 0 \n')
        fh.write('fix_kappa = 0 \n')
        fh.write('kappa = 2 \n')
        fh.write('fix_omega = 1 \n')
        fh.write('omega = 1 \n')
        fh.write('fix_alpha = 1 \n')
        fh.write('alpha = 0 \n')
        fh.write('Malpha = 0 \n')
        fh.write('ncatG = 8  \n')
        fh.write('getSE = 0 \n')
        fh.write('RateAncestor = 1 \n')
        fh.write('Small_Diff = .5e-6 \n')
        fh.write('method = 0 \n')

    with open(scriptfile_b, 'w') as fh:
        fh.write('#!/bin/bash\n')
        fh.write('cd '+codeml_prefix+' \n')
        fh.write('/home/amanda/software/paml4.9j/bin/codeml ')
        fh.write(control_file_name)
        fh.write('\n')

    ## submit the pbs.sh
    subprocess.call(['qsub', scriptfile_b])


def codeml_branches_runner(codeml_prefix,seq_file,treefile,outfile,scriptfile_b, gene_name):
    import subprocess
    #create the control file for codeml to run
    control_file_name = codeml_prefix  + gene_name + '_alt.ctl'
    with open(control_file_name, 'w') as fh:
        fh.write('seqfile = ' + seq_file +'\n')
        fh.write('treefile = '+ treefile +'\n')
        fh.write('outfile = ' + outfile +'\n')
        fh.write('noisy = 9 \n')
        fh.write('verbose = 0 \n')
        fh.write('runmode = 0 \n')
        fh.write('seqtype = 1 \n')
        fh.write('CodonFreq = 2 \n')
        fh.write('clock = 0 \n')
        fh.write('aaDist = 0 \n')
        fh.write('model = 2 \n')
        fh.write('NSsites = 2 \n')
        fh.write('Mgene = 0 \n')
        fh.write('fix_kappa = 0 \n')
        fh.write('kappa = 2 \n')
        fh.write('fix_omega = 0 \n')
        fh.write('omega = 0.4 \n')
        fh.write('fix_alpha = 1 \n')
        fh.write('alpha = 0 \n')
        fh.write('Malpha = 0 \n')
        fh.write('ncatG = 8  \n')
        fh.write('getSE = 0 \n')
        fh.write('RateAncestor = 1 \n')
        fh.write('Small_Diff = .5e-6 \n')
        fh.write('method = 0 \n')

    with open(scriptfile_b, 'w') as fh:
        fh.write('#!/bin/bash\n')
        fh.write('cd '+codeml_prefix+' \n')
        fh.write('/home/amanda/software/paml4.9j/bin/codeml ')
        fh.write(control_file_name)
        fh.write('\n')

    ## submit the pbs.sh
    subprocess.call(['qsub', scriptfile_b])


#Provide directory of the gene name subdirectories
dirName = sys.argv[1]

#Provide location of edited tree files
treeDir = sys.argv[2]

#nullInFile = "/home/leann/lib/longevity_dnds/null.ctl"
#altInFile = "/home/leann/lib/longevity_dnds/alternative.ctl"
#nullInFile = input("Input Blank Null file name: ")
#altInFile = input("Input Blank Alt file name: ")

dirList = os.listdir(dirName)
geneDirName = ''
geneName = ''


for items in dirList:
	geneDirName = dirName + items
	if(os.path.isdir(geneDirName)):
		geneName = items
		seqFileName = dirName+geneName+"/"+"nuc_"+geneName+"_aligned.phy"
		treeFileName = treeDir+geneName+"_edited_tree.nwk" #LOCATION OF EDITED TREEFILES
		nulloutFileName = dirName+geneName+"/"+geneName+"_null_branchsite.out"
		altoutFileName = dirName+geneName+"/"+geneName+"_alt_branchsite.out"
		altScriptName = dirName+geneName+"_alt_bsite_runcod.sh"
		nullScriptName = dirName+geneName+"_null_bsite_runcod.sh"
		codeml_branches_runner(dirName,seqFileName,treeFileName,altoutFileName,altScriptName,geneName)
		null_codeml_branches_runner(dirName,seqFileName,treeFileName,nulloutFileName,nullScriptName,geneName)