• API
• FAQ
• Tools
• Archive
daily pastebin goal
56%
SHARE
TWEET

# Untitled

a guest Oct 12th, 2018 66 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. import os
2. from numpy import log10
4. from cogent3.evolve.models import CNFGTR
5. from cogent3.maths.stats import chisqprob
6.
7. submodel = CNFGTR()
8.
9. def do_lr_test(lnL_alt, lnL_null, nfp_alt, nfp_null):
10.     LR = 2 * (lnL_alt - lnL_null)
11.     LR = LR if LR >= 0 else 0
12.     df = nfp_alt - nfp_null
13.     pval = chisqprob(LR, df)
14.     return pval
15.
16. def eval_alignment(path):
17.     # append the gene name, the p-vals from the two LRTs,
18.     # omega from the Human edge and the length from the Human edge
19.
20.     results = []
21.     gene_name = os.path.basename(path).split(".")[0]
22.     results.append(gene_name)
23.
26.     lf = submodel.make_likelihood_function(tree)
27.
28.     lf.set_alignment(aln)
29.
30.     ##### First LRT calc
31.     #define the neutral null
32.     lf.set_param_rule('omega', value=1.0, is_constant=True)
33.     lf.optimise()
34.     neutral_null_nfp = lf.get_num_free_params()
35.     neutral_null_lnL = lf.get_log_likelihood()
36.
37.     #define the neutral alt
38.     lf.set_param_rule('omega', value=1.0, is_constant=False)
39.     lf.optimise()
40.     neutral_alt_nfp = lf.get_num_free_params()
41.     neutral_alt_lnL = lf.get_log_likelihood()
42.
43.     results.append(do_lr_test(neutral_alt_lnL, neutral_null_lnL, neutral_alt_nfp, neutral_null_nfp))
44.
45.     ##### Second LRT calc
46.     #define the Human alternate
47.     lf.set_param_rule('omega', is_constant=True, edge="Human")
48.     human_alt_nfp = lf.get_num_free_params()
49.     human_alt_lnL = lf.get_log_likelihood()
50.
51.     results.append(do_lr_test(human_alt_lnL, neutral_null_lnL, human_alt_nfp, neutral_null_nfp))
52.     print(do_lr_test(human_alt_lnL, neutral_null_lnL, human_alt_nfp, neutral_null_nfp))
53.
54.
55.     results.append(lf.get_param_value("omega", edge="Human"))
56.     results.append(lf.get_param_value("length", edge="Human"))
57.
58.     return results
59.
60. paths = get_one2one_align_paths()
61. result = eval_alignment(paths[0])
62.
63. headers = ['Name', 'pval(Neutral)', 'pval(Human diff)', 'Omega from Human edge', 'Length from Human edge']