pylmm: Edit conflict

6 files changed, 72 insertions, 39 deletions

diff --git a/wqflask/runserver.py b/wqflask/runserver.py
index 9d5686a9..fadae6bf 100755
--- a/wqflask/runserver.py
+++ b/wqflask/runserver.py
@@ -20,9 +20,9 @@ from wqflask import app
 
 import logging
 #from themodule import TheHandlerYouWant
-file_handler = logging.FileHandler("/tmp/flask_gn_log_danny_unsecure")
-file_handler.setLevel(logging.DEBUG)
-app.logger.addHandler(file_handler)
+# file_handler = logging.FileHandler("/tmp/flask_gn_log_danny_unsecure")
+# file_handler.setLevel(logging.DEBUG)
+# app.logger.addHandler(file_handler)
 
 import logging_tree
 logging_tree.printout()
diff --git a/wqflask/wqflask/marker_regression/marker_regression.py b/wqflask/wqflask/marker_regression/marker_regression.py
index 7708356b..fba34b99 100755
--- a/wqflask/wqflask/marker_regression/marker_regression.py
+++ b/wqflask/wqflask/marker_regression/marker_regression.py
@@ -40,6 +40,7 @@ from utility import temp_data
 
 from utility.benchmark import Bench
 
+PYLMM_COMMAND= 'python /home/pjotr/izip/git/opensource/python/gn2/wqflask/wqflask/my_pylmm/pyLMM/lmm.py'
 
 class MarkerRegression(object):
 
@@ -126,8 +127,9 @@ class MarkerRegression(object):
 
         #Need to convert the QTL objects that qtl reaper returns into a json serializable dictionary
         self.qtl_results = []
-        for qtl in self.filtered_markers:
-            print("lod score is:", qtl['lod_score'])
+        for index,qtl in enumerate(self.filtered_markers):
+            if index<40:
+                print("lod score is:", qtl['lod_score'])
             if qtl['chr'] == highest_chr and highest_chr != "X" and highest_chr != "X/Y":
                 print("changing to X")
                 self.json_data['chr'].append("X")
@@ -144,7 +146,7 @@ class MarkerRegression(object):
             self.json_data['chrnames'].append([self.species.chromosomes.chromosomes[key].name, self.species.chromosomes.chromosomes[key].mb_length])
             chromosome_mb_lengths[key] = self.species.chromosomes.chromosomes[key].mb_length
         
-        print("json_data:", self.json_data)
+        # print("json_data:", self.json_data)
         
 
         self.js_data = dict(
@@ -272,7 +274,7 @@ class MarkerRegression(object):
         """)
     
     def run_rqtl_geno(self):
-        print("Calling R/qtl from python")
+        print("Calling R/qtl")
 
         self.geno_to_rqtl_function()
 
@@ -655,8 +657,7 @@ class MarkerRegression(object):
                 Redis.set(key, json_params)
                 Redis.expire(key, 60*60)
     
-                command = 'python /home/zas1024/gene/wqflask/wqflask/my_pylmm/pyLMM/lmm.py --key {} --species {}'.format(key,
-                                                                                                                        "other")
+                command = PYLMM_COMMAND+' --key {} --species {}'.format(key,"other")
     
                 os.system(command)
     
@@ -713,8 +714,8 @@ class MarkerRegression(object):
             #            "refit": False,
             #            "temp_data": tempdata}
             
-            print("genotype_matrix:", str(genotype_matrix.tolist()))
-            print("pheno_vector:", str(pheno_vector.tolist()))
+            # print("genotype_matrix:", str(genotype_matrix.tolist()))
+            # print("pheno_vector:", str(pheno_vector.tolist()))
             
             params = dict(pheno_vector = pheno_vector.tolist(),
                         genotype_matrix = genotype_matrix.tolist(),
@@ -732,7 +733,7 @@ class MarkerRegression(object):
             Redis.expire(key, 60*60)
             print("before printing command")
 
-            command = 'python /home/zas1024/gene/wqflask/wqflask/my_pylmm/pyLMM/lmm.py --key {} --species {}'.format(key,
+            command = PYLMM_COMMAND + ' --key {} --species {}'.format(key,
                                                                                                                     "other")
             print("command is:", command)
             print("after printing command")
@@ -745,7 +746,7 @@ class MarkerRegression(object):
             json_results = Redis.blpop("pylmm:results:" + temp_uuid, 45*60)
             results = json.loads(json_results[1])
             p_values = [float(result) for result in results['p_values']]
-            print("p_values:", p_values)
+            print("p_values:", p_values[:10])
             #p_values = self.trim_results(p_values)
             t_stats = results['t_stats']
             
@@ -806,7 +807,7 @@ class MarkerRegression(object):
 
         print("Before creating the command")
 
-        command = 'python /home/zas1024/gene/wqflask/wqflask/my_pylmm/pyLMM/lmm.py --key {} --species {}'.format(key,
+        command = PYLMM_COMMAND+' --key {} --species {}'.format(key,
                                                                                                                 "human")
         
         print("command is:", command)
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/gn2.py b/wqflask/wqflask/my_pylmm/pyLMM/gn2.py
index 40b2021d..821195c8 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/gn2.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/gn2.py
@@ -17,7 +17,7 @@ logger = logging.getLogger('lmm2')
 logging.basicConfig(level=logging.DEBUG)
 np.set_printoptions(precision=3,suppress=True)
 
-progress_location = None 
+progress_location = None
 progress_current  = None
 progress_prev_perc     = None
 
@@ -25,17 +25,17 @@ def progress_default_func(location,count,total):
     global progress_current
     value = round(count*100.0/total)
     progress_current = value
-    
+
 progress_func = progress_default_func
 
 def progress_set_func(func):
     global progress_func
     progress_func = func
-    
+
 def progress(location, count, total):
     global progress_location
     global progress_prev_perc
-    
+
     perc = round(count*100.0/total)
     if perc != progress_prev_perc and (location != progress_location or perc > 98 or perc > progress_prev_perc + 5):
         progress_func(location, count, total)
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py b/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
index 98bbead8..5b06c9ae 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
@@ -42,25 +42,27 @@ from redis import Redis
 Redis = Redis()
 
 import sys
-sys.path.append("/home/zas1024/gene/wqflask/")
-
-has_gn2=True
 
 from utility.benchmark import Bench
 from utility import temp_data
 
-sys.path.append("/home/zas1024/gene/wqflask/wqflask/my_pylmm/pyLMM/")
-
 from kinship import kinship, kinship_full, kvakve
 import genotype
 import phenotype
 import gwas
 
+has_gn2=True
+sys.stderr.write("INFO: pylmm system path is "+":".join(sys.path)+"\n")
+sys.stderr.write("INFO: pylmm file is "+__file__+"\n")
+
 # ---- A trick to decide on the environment:
 try:
-    from wqflask.my_pylmm.pyLMM import chunks
+    sys.stderr.write("INFO: trying loading module\n")
+    import utility.formatting # this is never used, just to check the environment
+    sys.stderr.write("INFO: This is a genenetwork2 environment\n")
     from gn2 import uses, progress_set_func
 except ImportError:
+    # Failed to load gn2
     has_gn2=False
     import standalone as handlers
     from standalone import uses, progress_set_func
@@ -852,30 +854,49 @@ def gwas_with_redis(key,species,new_code=True):
         debug("Updating REDIS percent_complete=%d" % (round(i*100.0/total)))
     progress_set_func(update_tempdata)
 
-    def narray(key):
-        print(key)
-        v = params[key]
+    def narray(t):
+        info("Type is "+t)
+        v = params.get(t)
         if v is not None:
-            v = np.array(v).astype(np.float)
+            # Note input values can be array of string or float
+            v1 = [x if x != 'NA' else 'nan' for x in v]
+            v = np.array(v1).astype(np.float)
         return v
 
-    def narrayT(key):
-        m = narray(key)
+    def marray(t):
+        info("Type is "+t)
+        v = params.get(t)
+        if v is not None:
+            m = []
+            for r in v:
+              # Note input values can be array of string or float
+              r1 = [x if x != 'NA' else 'nan' for x in r]
+              m.append(np.array(r1).astype(np.float))
+            return np.array(m)
+        return np.array(v)
+
+    def marrayT(t):
+        m = marray(t)
         if m is not None:
             return m.T
         return m
     
     # We are transposing before we enter run_gwas - this should happen on the webserver
     # side (or when reading data from file)
-    k = narray('kinship_matrix')
-    g = narrayT('genotype_matrix')
+    k = marray('kinship_matrix')
+    g = marrayT('genotype_matrix')
+    mprint("geno",g)
     y = narray('pheno_vector')
     n = len(y)
-    m = params['num_genotypes']
-    ps,ts = run_gwas(species,n,m,k,y,g,narray('covariate_matrix'),params['restricted_max_likelihood'],params['refit'],params['input_file_name'],new_code)
+    m = params.get('num_genotypes')
+    if m is None:
+      m = g.shape[0]
+    info("m=%d,n=%d" % (m,n))
+    ps,ts = run_gwas(species,n,m,k,y,g,narray('covariate_matrix'),params['restricted_max_likelihood'],params['refit'],params.get('input_file_name'),new_code)
         
     results_key = "pylmm:results:" + params['temp_uuid']
 
+    # fatal(results_key)
     json_results = json.dumps(dict(p_values = ps,
                                    t_stats = ts))
     
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/lmm2.py b/wqflask/wqflask/my_pylmm/pyLMM/lmm2.py
index 358bf27e..c65843ec 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/lmm2.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/lmm2.py
@@ -24,14 +24,24 @@ from scipy import optimize
 from optmatrix import matrixMult
 import kinship
 
+sys.stderr.write("INFO: pylmm (lmm2) system path is "+":".join(sys.path)+"\n")
+sys.stderr.write("INFO: pylmm (lmm2) file is "+__file__+"\n")
+
 # ---- A trick to decide on the environment:
 try:
-    from wqflask.my_pylmm.pyLMM import chunks
-    from gn2 import uses
+    sys.stderr.write("INFO: trying loading module\n")
+    import utility.formatting # this is never used, just to check the environment
+    sys.stderr.write("INFO: This is a genenetwork2 environment (lmm2)\n")
+    from gn2 import uses, progress_set_func
 except ImportError:
-    sys.stderr.write("WARNING: LMM2 standalone version missing the Genenetwork2 environment\n")
+    # Failed to load gn2
     has_gn2=False
-    from standalone import uses
+    import standalone as handlers
+    from standalone import uses, progress_set_func
+    sys.stderr.write("WARNING: LMM2 standalone version missing the Genenetwork2 environment\n")
+
+progress,mprint,debug,info,fatal = uses('progress','mprint','debug','info','fatal')
+
 
 def calculateKinship(W,center=False):
       """
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/tsvreader.py b/wqflask/wqflask/my_pylmm/pyLMM/tsvreader.py
index 27daf43f..66b34ee2 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/tsvreader.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/tsvreader.py
@@ -56,7 +56,8 @@ def geno(fn):
 
     print fn
     with open(fn,'r') as tsvin:
-        assert(tsvin.readline().strip() == "# Genotype format version 1.0")
+        line = tsvin.readline().strip()
+        assert line == "# Genotype format version 1.0", line
         tsvin.readline()
         tsvin.readline()
         tsvin.readline()