From a5d0d42f6977c4473c28f12999a62ad360750041 Mon Sep 17 00:00:00 2001
From: Pjotr Prins
Date: Wed, 11 Mar 2015 12:07:54 +0300
Subject: & Added testing option, now the results match

---
 wqflask/wqflask/my_pylmm/pyLMM/lmm.py    | 16 ++++++++++------
 wqflask/wqflask/my_pylmm/pyLMM/runlmm.py | 11 ++++++++---
 2 files changed, 18 insertions(+), 9 deletions(-)

(limited to 'wqflask')

diff --git a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py b/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
index ab19bf08..5ffb9a4c 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
@@ -267,6 +267,7 @@ def run_other(pheno_vector,
     """
     
     print("In run_other")
+    print("REML=",restricted_max_likelihood," REFIT=",refit, "TESTING=",is_testing)
     with Bench("Calculate Kinship"):
         kinship_matrix = calculate_kinship(genotype_matrix, tempdata, is_testing)
     
@@ -280,6 +281,7 @@ def run_other(pheno_vector,
         lmm_ob.fit()
 
     print("genotype_matrix: ", genotype_matrix.shape)
+    print(genotype_matrix)
 
     with Bench("Doing GWAS"):
         t_stats, p_values = GWAS(pheno_vector,
@@ -720,7 +722,7 @@ class LMM:
        pl.title(title)
 
 
-def gn2_redis(key,species):
+def gn2_redis(key,species,is_testing=False):
     json_params = Redis.get(key)
     
     params = json.loads(json_params)
@@ -729,7 +731,8 @@ def gn2_redis(key,species):
 
     print('kinship', np.array(params['kinship_matrix']))
     print('pheno', np.array(params['pheno_vector']))
-    print('geno', np.array(params['genotype_matrix']))
+    geno = np.array(params['genotype_matrix'])
+    print('geno', geno.shape, geno)
     # sys.exit(1)
     
     if species == "human" :
@@ -741,10 +744,11 @@ def gn2_redis(key,species):
                   tempdata = tempdata)
     else:
         ps, ts = run_other(pheno_vector = np.array(params['pheno_vector']),
-                  genotype_matrix = np.array(params['genotype_matrix']),
+                  genotype_matrix = geno,
                   restricted_max_likelihood = params['restricted_max_likelihood'],
                   refit = params['refit'],
-                  tempdata = tempdata, is_testing=False)
+                  tempdata = tempdata,
+                  is_testing=is_testing)
         
     results_key = "pylmm:results:" + params['temp_uuid']
 
@@ -769,7 +773,7 @@ def gn2_main():
 
     gn2_redis(key,species)
 
-def gn2_load_redis(key,species,kinship,pheno,geno):
+def gn2_load_redis(key,species,kinship,pheno,geno,is_testing):
     print("Loading Redis from parsed data")
     params = dict(pheno_vector = pheno.tolist(),
                   genotype_matrix = geno.tolist(),
@@ -786,7 +790,7 @@ def gn2_load_redis(key,species,kinship,pheno,geno):
     Redis.set(key, json_params)
     Redis.expire(key, 60*60)
 
-    return gn2_redis(key,species)
+    return gn2_redis(key,species,is_testing)
     
 if __name__ == '__main__':
     print("WARNING: Calling pylmm from lmm.py will become OBSOLETE, use runlmm.py instead!")
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
index 907a6835..738268be 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
@@ -50,6 +50,9 @@ parser.add_option("--geno",dest="geno",
 parser.add_option("-q", "--quiet",
                   action="store_false", dest="verbose", default=True,
                   help="don't print status messages to stdout")
+parser.add_option("--test",
+                  action="store_true", dest="testing", default=False,
+                  help="Testing mode")
 
 (options, args) = parser.parse_args()
 
@@ -76,9 +79,9 @@ def normalizeGenotype(G):
     # Run for every SNP list (num individuals)
     x = True - np.isnan(G)  # Matrix of True/False
     m = G[x].mean()         # Global mean value
-    print m
+    # print m
     s = np.sqrt(G[x].var()) # Global stddev
-    print s
+    # print s
     G[np.isnan(G)] = m      # Plug-in mean values for missing
     if s == 0:
         G = G - m           # Subtract the mean
@@ -92,6 +95,7 @@ gn = []
 for snp in g:
     gn.append( normalizeGenotype(snp) )
 
+gn = g
 gn = np.array(gn)
 print("After1",gn)
 gnT = gn.T
@@ -99,6 +103,7 @@ print("After",gnT)
 # G = gnT
 G = gnT
 print "G shape",G.shape
+# sys.exit(1)
 # assert(G[0,0]==-2.25726341)
 
 # Remove individuals with missing phenotypes
@@ -111,5 +116,5 @@ if v.sum():
    k = k[keep,:][:,keep]
 
 if cmd == 'redis':
-    ps, ts = gn2_load_redis('testrun','other',k,y,G)
+    ps, ts = gn2_load_redis('testrun','other',np.array(k),y,G,options.testing)
     print np.array(ps)
-- 
cgit v1.2.3