@@ -24,6 +24,9 @@ def main(args):
24
24
save_file = './save/PSJS_HKY_' + '_' .join (paralog ) + '_' + IGC_pm .replace (' ' , '_' ) + '_Guess_' + str (args .guess ) + '_nonclock_save.txt'
25
25
log_file = './log/PSJS_HKY_' + '_' .join (paralog ) + '_' + IGC_pm .replace (' ' , '_' ) + '_Guess_' + str (args .guess ) + '_nonclock_log.txt'
26
26
summary_file = './summary/PSJS_HKY_' + '_' .join (paralog ) + '_' + IGC_pm .replace (' ' , '_' ) + '_Guess_' + str (args .guess ) + '_nonclock_summary.txt'
27
+ gradient_file = './summary/PSJS_HKY_' + '_' .join (paralog ) + '_' + IGC_pm .replace (' ' , '_' ) + '_Guess_' + str (args .guess ) + '_nonclock_gradient.txt'
28
+ hessian_file = './summary/PSJS_HKY_' + '_' .join (paralog ) + '_' + IGC_pm .replace (' ' , '_' ) + '_Guess_' + str (args .guess ) + '_nonclock_hessian.txt'
29
+ godambe_file = './summary/PSJS_HKY_' + '_' .join (paralog ) + '_' + IGC_pm .replace (' ' , '_' ) + '_Guess_' + str (args .guess ) + '_nonclock_godambe.txt'
27
30
28
31
initial_tract_length_list = np .log ([30.0 , 200.0 ])
29
32
guess_lnp = - initial_tract_length_list [args .guess - 1 ]
@@ -50,76 +53,80 @@ def main(args):
50
53
#test.unpack_x(x)
51
54
52
55
53
- test .get_mle ()
56
+ test .get_mle (stringent_level = 'high' )
54
57
test .get_individual_summary (summary_file )
55
58
59
+ Godambe_x = np .array ([test .psjsmodel .x_IGC [0 ] - test .psjsmodel .x_IGC [1 ], test .psjsmodel .x_IGC [1 ] - np .log (1.0 - np .exp (test .psjsmodel .x_IGC [1 ]))])
60
+ godambe = test .get_Godambe_matrix (Godambe_x )
61
+ np .savetxt (open (godambe_file , 'w+' ), np .array (godambe ))
62
+ test .get_gradient_hessian (Godambe_x , gradient_file , hessian_file )
56
63
57
64
58
65
if __name__ == '__main__' :
59
- ## parser = argparse.ArgumentParser()
60
- ## parser.add_argument('--paralog1', required = True, help = 'Name of the 1st paralog')
61
- ## parser.add_argument('--paralog2', required = True, help = 'Name of the 2nd paralog')
62
- ## parser.add_argument('--G', type = int, dest = 'guess', default = 1, help = 'Guess case')
63
- ## parser.add_argument('--heterogeneity', dest = 'rate_variation', action = 'store_true', help = 'rate heterogeneity control')
64
- ## parser.add_argument('--homogeneity', dest = 'rate_variation', action = 'store_false', help = 'rate heterogeneity control')
65
- ## parser.add_argument('--coding', dest = 'cdna', action = 'store_true', help = 'coding sequence control')
66
- ## parser.add_argument('--noncoding', dest = 'cdna', action = 'store_false', help = 'coding sequence control')
67
- ## parser.add_argument('--samecodon', dest = 'allow_same_codon', action = 'store_true', help = 'whether allow pair sites from same codon')
68
- ## parser.add_argument('--no-samecodon', dest = 'allow_same_codon', action = 'store_false', help = 'whether allow pair sites from same codon')
69
- ##
70
- ## main(parser.parse_args())
71
-
72
-
73
-
74
- MyStruct = namedtuple ('MyStruct' , 'paralog1 paralog2 tract_length dim cdna rate_variation allow_same_codon guess' )
75
- args = MyStruct (paralog1 = 'EDN' , paralog2 = 'ECP' , tract_length = 30.0 , dim = 1 ,
76
- cdna = True , rate_variation = True , allow_same_codon = True ,
77
- guess = 2 )
78
-
79
- paralog = [args .paralog1 , args .paralog2 ]
80
-
81
- gene_to_orlg_file = './' + '_' .join (paralog ) + '_GeneToOrlg.txt'
82
- alignment_file = './' + '_' .join (paralog ) + '_Cleaned_input.fasta'
83
-
84
- tree_newick = './input_tree.newick'
85
- DupLosList = './EDN_ECP_DupLost.txt'
86
- terminal_node_list = ['Tamarin' , 'Macaque' , 'Orangutan' , 'Gorilla' , 'Chimpanzee' ]
87
- node_to_pos = {'D1' :0 }
88
- pm_model = 'HKY'
89
-
90
- IGC_pm = 'One rate'
91
- seq_index_file = './' + '_' .join (paralog ) + '_seq_index.txt'
92
- averaged_x = np .loadtxt ('./EDN_ECP_JS_HKY_One_rate_nonclock_save.txt' )
93
- averaged_x_js , averaged_x_rates = averaged_x [:5 ], averaged_x [5 :]
94
-
95
- save_file = './save/PSJS_HKY_' + '_' .join (paralog ) + '_' + IGC_pm .replace (' ' , '_' ) + '_Guess_' + str (args .guess ) + '_nonclock_save.txt'
96
- log_file = './log/PSJS_HKY_' + '_' .join (paralog ) + '_' + IGC_pm .replace (' ' , '_' ) + '_Guess_' + str (args .guess ) + '_nonclock_log.txt'
97
- summary_file = './summary/PSJS_HKY_' + '_' .join (paralog ) + '_' + IGC_pm .replace (' ' , '_' ) + '_Guess_' + str (args .guess ) + '_nonclock_summary.txt'
98
-
99
- initial_tract_length_list = np .log ([30.0 , 200.0 ])
100
- guess_lnp = - initial_tract_length_list [args .guess - 1 ]
66
+ parser = argparse .ArgumentParser ()
67
+ parser .add_argument ('--paralog1' , required = True , help = 'Name of the 1st paralog' )
68
+ parser .add_argument ('--paralog2' , required = True , help = 'Name of the 2nd paralog' )
69
+ parser .add_argument ('--G' , type = int , dest = 'guess' , default = 1 , help = 'Guess case' )
70
+ parser .add_argument ('--heterogeneity' , dest = 'rate_variation' , action = 'store_true' , help = 'rate heterogeneity control' )
71
+ parser .add_argument ('--homogeneity' , dest = 'rate_variation' , action = 'store_false' , help = 'rate heterogeneity control' )
72
+ parser .add_argument ('--coding' , dest = 'cdna' , action = 'store_true' , help = 'coding sequence control' )
73
+ parser .add_argument ('--noncoding' , dest = 'cdna' , action = 'store_false' , help = 'coding sequence control' )
74
+ parser .add_argument ('--samecodon' , dest = 'allow_same_codon' , action = 'store_true' , help = 'whether allow pair sites from same codon' )
75
+ parser .add_argument ('--no-samecodon' , dest = 'allow_same_codon' , action = 'store_false' , help = 'whether allow pair sites from same codon' )
101
76
102
- if args .rate_variation :
103
- x_js = np .concatenate ((averaged_x_js [:- 1 ], np .log ([0.7 , 3.0 ]), [averaged_x_js [- 1 ] + guess_lnp , guess_lnp ]))
104
- if args .allow_same_codon :
105
- save_file = save_file .replace ('_nonclock' , '_rv_SCOK_nonclock' )
106
- log_file = log_file .replace ('_nonclock' , '_rv_SCOK_nonclock' )
107
- summary_file = summary_file .replace ('_nonclock' , '_rv_SCOK_nonclock' )
108
- else :
109
- save_file = save_file .replace ('_nonclock' , '_rv_NOSC_nonclock' )
110
- log_file = log_file .replace ('_nonclock' , '_rv_SCOK_nonclock' )
111
- summary_file = summary_file .replace ('_nonclock' , '_rv_NOSC_nonclock' )
112
- else :
113
- x_js = np .concatenate ((averaged_x_js [:- 1 ], [averaged_x_js [- 1 ] + guess_lnp , guess_lnp ]))
114
-
115
-
116
-
117
-
118
- test = PSJSGeneconv (alignment_file , gene_to_orlg_file , seq_index_file , args .cdna , args .allow_same_codon , tree_newick , DupLosList , x_js , pm_model , IGC_pm ,
119
- args .rate_variation , node_to_pos , terminal_node_list , save_file , log_file )
120
- #x = np.concatenate((x_js, averaged_x_rates))
121
- #test.unpack_x(x)
77
+ main (parser .parse_args ())
122
78
79
+
123
80
124
- test .get_mle ()
125
- test .get_individual_summary (summary_file )
81
+ ## MyStruct = namedtuple('MyStruct', 'paralog1 paralog2 tract_length dim cdna rate_variation allow_same_codon guess')
82
+ ## args = MyStruct(paralog1 = 'EDN', paralog2 = 'ECP', tract_length = 30.0, dim = 1,
83
+ ## cdna = True, rate_variation = True, allow_same_codon = True,
84
+ ## guess = 2)
85
+ ##
86
+ ## paralog = [args.paralog1, args.paralog2]
87
+ ##
88
+ ## gene_to_orlg_file = './' + '_'.join(paralog) +'_GeneToOrlg.txt'
89
+ ## alignment_file = './' + '_'.join(paralog) +'_Cleaned_input.fasta'
90
+ ##
91
+ ## tree_newick = './input_tree.newick'
92
+ ## DupLosList = './EDN_ECP_DupLost.txt'
93
+ ## terminal_node_list = ['Tamarin', 'Macaque', 'Orangutan', 'Gorilla', 'Chimpanzee']
94
+ ## node_to_pos = {'D1':0}
95
+ ## pm_model = 'HKY'
96
+ ##
97
+ ## IGC_pm = 'One rate'
98
+ ## seq_index_file = './' + '_'.join(paralog) +'_seq_index.txt'
99
+ ## averaged_x = np.loadtxt('./EDN_ECP_JS_HKY_One_rate_nonclock_save.txt')
100
+ ## averaged_x_js, averaged_x_rates = averaged_x[:5], averaged_x[5:]
101
+ ##
102
+ ## save_file = './save/PSJS_HKY_'+ '_'.join(paralog) + '_' + IGC_pm.replace(' ', '_') + '_Guess_' + str(args.guess) + '_nonclock_save.txt'
103
+ ## log_file = './log/PSJS_HKY_'+ '_'.join(paralog) + '_' + IGC_pm.replace(' ', '_') + '_Guess_' + str(args.guess) + '_nonclock_log.txt'
104
+ ## summary_file = './summary/PSJS_HKY_'+ '_'.join(paralog) + '_' + IGC_pm.replace(' ', '_') + '_Guess_' + str(args.guess) + '_nonclock_summary.txt'
105
+ ##
106
+ ## initial_tract_length_list = np.log([30.0, 200.0])
107
+ ## guess_lnp = -initial_tract_length_list[args.guess-1]
108
+ ##
109
+ ## if args.rate_variation:
110
+ ## x_js = np.concatenate((averaged_x_js[:-1], np.log([0.7, 3.0]), [averaged_x_js[-1] + guess_lnp, guess_lnp]))
111
+ ## if args.allow_same_codon:
112
+ ## save_file = save_file.replace('_nonclock', '_rv_SCOK_nonclock')
113
+ ## log_file = log_file.replace('_nonclock', '_rv_SCOK_nonclock')
114
+ ## summary_file = summary_file.replace('_nonclock', '_rv_SCOK_nonclock')
115
+ ## else:
116
+ ## save_file = save_file.replace('_nonclock', '_rv_NOSC_nonclock')
117
+ ## log_file = log_file.replace('_nonclock', '_rv_SCOK_nonclock')
118
+ ## summary_file = summary_file.replace('_nonclock', '_rv_NOSC_nonclock')
119
+ ## else:
120
+ ## x_js = np.concatenate((averaged_x_js[:-1], [averaged_x_js[-1] + guess_lnp, guess_lnp]))
121
+ ##
122
+ ##
123
+ ##
124
+ ##
125
+ ## test = PSJSGeneconv(alignment_file, gene_to_orlg_file, seq_index_file, args.cdna, args.allow_same_codon, tree_newick, DupLosList, x_js, pm_model, IGC_pm,
126
+ ## args.rate_variation, node_to_pos, terminal_node_list, save_file, log_file)
127
+ ## #x = np.concatenate((x_js, averaged_x_rates))
128
+ ## #test.unpack_x(x)
129
+ ##
130
+ ##
131
+ ## test.get_mle()
132
+ ## test.get_individual_summary(summary_file)
0 commit comments