aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/lmm.py9
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/runlmm.py25
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:]