Missing genotypes and normalize

2 files changed, 32 insertions, 2 deletions

diff --git a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py b/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
index 8c6e30e8..00bbf144 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
@@ -726,6 +726,11 @@ def gn2_redis(key,species):
     params = json.loads(json_params)
     
     tempdata = temp_data.TempData(params['temp_uuid'])
+
+    print('kinship', np.array(params['kinship_matrix'][0:10][0:10]))
+    print('pheno', params['pheno_vector'][0:10])
+    print('geno', params['genotype_matrix'][0:10][0:10])
+    
     if species == "human" :
         ps, ts = run_human(pheno_vector = np.array(params['pheno_vector']),
                   covariate_matrix = np.array(params['covariate_matrix']),
@@ -748,6 +753,7 @@ def gn2_redis(key,species):
     #Pushing json_results into a list where it is the only item because blpop needs a list
     Redis.rpush(results_key, json_results)
     Redis.expire(results_key, 60*60)
+    return ps, ts
 
 # This is the main function used by Genenetwork2 (with environment)
 def gn2_main():
@@ -766,6 +772,7 @@ def gn2_load_redis(key,species,kinship,pheno,geno):
     print("Loading Redis from parsed data")
     params = dict(pheno_vector = pheno.tolist(),
                   genotype_matrix = geno.tolist(),
+                  kinship_matrix= kinship.tolist(),
                   restricted_max_likelihood = True,
                   refit = False,
                   temp_uuid = "testrun_temp_uuid",
@@ -778,7 +785,7 @@ def gn2_load_redis(key,species,kinship,pheno,geno):
     Redis.set(key, json_params)
     Redis.expire(key, 60*60)
 
-    gn2_redis(key,species)
+    return gn2_redis(key,species)
     
 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 525c43fc..6db1bdbc 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
@@ -18,7 +18,9 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 from optparse import OptionParser
+import sys
 import tsvreader
+import numpy as np
 from lmm import gn2_load_redis
 
 usage = """
@@ -70,5 +72,26 @@ if options.geno:
     g = tsvreader.geno(options.geno)
     print g.shape
 
+def normalizeGenotype(G):
+    x = True - np.isnan(G)
+    m = G[x].mean()
+    s = np.sqrt(G[x].var())
+    G[np.isnan(G)] = m
+    if s == 0: G = G - m
+    else: G = (G - m) / s
+    return G
+
+gT = normalizeGenotype(g.T)
+
+# Remove individuals with missing phenotypes
+v = np.isnan(y)
+keep = True - v
+if v.sum():
+   if options.verbose: sys.stderr.write("Cleaning the phenotype vector by removing %d individuals...\n" % (v.sum()))
+   y  = y[keep]
+   gT = gT[keep,:]
+   k = k[keep,:][:,keep]
+
 if cmd == 'redis':
-    gn2_load_redis('testrun','other',k,y,g.T)
+    ps, ts = gn2_load_redis('testrun','other',k,y,gT)
+    print ps[0:10],ps[-10:]