about summary refs log tree commit diff
diff options
context:
space:
mode:
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/lmm.py7
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/runlmm.py28
2 files changed, 22 insertions, 13 deletions
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py b/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
index 87b999e8..d6830787 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
@@ -338,6 +338,7 @@ def calculate_kinship(genotype_matrix, temp_data=None):
     m = genotype_matrix.shape[1]
     print("genotype 2D matrix n (inds) is:", n)
     print("genotype 2D matrix m (snps) is:", m)
+    assert m>n, "n should be larger than m (snps>inds)"
     keep = []
     for counter in range(m):
         #print("type of genotype_matrix[:,counter]:", pf(genotype_matrix[:,counter]))
@@ -777,9 +778,13 @@ def gn2_main():
 
 def gn2_load_redis(key,species,kinship,pheno,geno):
     print("Loading Redis from parsed data")
+    if kinship == None:
+        k = None
+    else:
+        k = kinship.tolist()
     params = dict(pheno_vector = pheno.tolist(),
                   genotype_matrix = geno.tolist(),
-                  kinship_matrix= kinship.tolist(),
+                  kinship_matrix= k,
                   restricted_max_likelihood = True,
                   refit = False,
                   temp_uuid = "testrun_temp_uuid",
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
index 80478368..7a77ad0a 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
@@ -98,16 +98,20 @@ if options.geno:
 
 if cmd == 'redis':
     # Emulating the redis setup of GN2
-    gn = []
-    for ind_g in g:
-        gn.append( genotype.normalize(ind_g) )
-    gnt = np.array(gn).T
-    if y:
-        Y,G = phenotype.remove_missing(y,gnt,options.verbose)
-        print "G",G.shape,G
-    else:
-        G = gnt
-    ps, ts = gn2_load_redis('testrun','other',k,Y,G,options.testing)
+    G = g
+    print "Original G",G.shape, "\n", G
+    if y != None:
+        gnt = np.array(g).T
+        Y,g = phenotype.remove_missing(y,g.T,options.verbose)
+        G = g.T
+        print "Removed missing phenotypes",G.shape, "\n", G
+    if options.maf_normalization:
+        G = np.apply_along_axis( genotype.replace_missing_with_MAF, axis=0, arr=g )
+        print "MAF replacements: \n",G
+    if not options.skip_genotype_normalization:
+        G = np.apply_along_axis( genotype.normalize, axis=1, arr=G)
+
+    ps, ts = gn2_load_redis('testrun','other',k,Y,G.T)
     print np.array(ps)
     print round(ps[0],4)
     assert(options.testing and round(ps[0],4)==0.7262)
@@ -116,8 +120,8 @@ if cmd == 'redis':
 elif cmd == 'kinship':
     G = g
     print "Original G",G.shape, "\n", G
-    if y:
-        gnt = np.array(gn).T
+    if y != None:
+        gnt = np.array(g).T
         Y,g = phenotype.remove_missing(y,g.T,options.verbose)
         G = g.T
         print "Removed missing phenotypes",G.shape, "\n", G