daily pastebin goal
16%
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
  3. from cogent3 import LoadSeqs, LoadTree, DNA, LoadTable
  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.  
  24.     aln = LoadSeqs(path, moltype=DNA)
  25.     tree = LoadTree(tip_names=aln.names)
  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']
  64. tbl = LoadTable(header= headers, rows= [result])
  65. tbl
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top