From e9865707ef447b8bc23eb8c872703f156936499d Mon Sep 17 00:00:00 2001
From: Pjotr Prins
Date: Fri, 3 Apr 2015 14:03:32 +0200
Subject: - Calculate n,m from the start - added test function to runlmm.py to
 run without Redis (25% faster)

---
 wqflask/wqflask/my_pylmm/pyLMM/lmm.py    | 17 ++++++++++-------
 wqflask/wqflask/my_pylmm/pyLMM/runlmm.py | 17 +++++++++++++++--
 2 files changed, 25 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 ad6375e9..e51742c4 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
@@ -795,9 +795,9 @@ class LMM:
        pl.title(title)
 
 
-def run_gwas(species,k,y,geno,cov,reml,refit,inputfn,new_code):
+def run_gwas(species,n,m,k,y,geno,cov=None,reml=True,refit=False,inputfn=None,new_code=True):
     """
-    Invoke pylmm using a genotype (SNP) iterator
+    Invoke pylmm using genotype as a matrix or as a (SNP) iterator.
     """
     info("gwas_without_redis")
     print('pheno', y)
@@ -848,8 +848,11 @@ def gwas_with_redis(key,species,new_code=True):
         if v is not None:
             v = np.array(v)
         return v
-    
-    ps,ts = run_gwas(species,narray('kinship_matrix'),narray('pheno_vector'),narray('genotype_matrix'),narray('covariate_matrix'),params['restricted_max_likelihood'],params['refit'],params['input_file_name'],new_code)
+
+    y = narray('pheno_vector')
+    n = len(y)
+    m = params['num_genotypes']
+    ps,ts = run_gwas(species,n,m,narray('kinship_matrix'),y,narray('genotype_matrix'),narray('covariate_matrix'),params['restricted_max_likelihood'],params['refit'],params['input_file_name'],new_code)
         
     results_key = "pylmm:results:" + params['temp_uuid']
 
@@ -873,6 +876,7 @@ def gn2_load_redis(key,species,kinship,pheno,geno,new_code=True):
         k = kinship.tolist()
     params = dict(pheno_vector = pheno.tolist(),
                   genotype_matrix = geno.tolist(),
+                  num_genotypes = geno.shape[1],
                   kinship_matrix = k,
                   covariate_matrix = None,
                   input_file_name = None,
@@ -881,8 +885,7 @@ def gn2_load_redis(key,species,kinship,pheno,geno,new_code=True):
                   temp_uuid = "testrun_temp_uuid",
                         
                   # meta data
-                  timestamp = datetime.datetime.now().isoformat(),
-    )
+                  timestamp = datetime.datetime.now().isoformat())
             
     json_params = json.dumps(params)
     Redis.set(key, json_params)
@@ -907,7 +910,7 @@ def gn2_load_redis_iter(key,species,kinship,pheno,geno_iterator):
         k = kinship.tolist()
     params = dict(pheno_vector = pheno.tolist(),
                   genotype_matrix = "iterator",
-                  genotypes = i,
+                  num_genotypes = i,
                   kinship_matrix = k,
                   covariate_matrix = None,
                   input_file_name = None,
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
index 3801529e..f095bb73 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
@@ -21,7 +21,7 @@ from optparse import OptionParser
 import sys
 import tsvreader
 import numpy as np
-from lmm import gn2_load_redis, gn2_load_redis_iter, calculate_kinship_new
+from lmm import gn2_load_redis, gn2_load_redis_iter, calculate_kinship_new, run_gwas
 from kinship import kinship, kinship_full
 import genotype
 import phenotype
@@ -103,7 +103,20 @@ if options.geno and cmd != 'iterator':
     g = tsvreader.geno(options.geno)
     print g.shape
 
-if cmd == 'iterator':
+if cmd == 'run':
+    if options.remove_missing_phenotypes:
+        raise Exception('Can not use --remove-missing-phenotypes with LMM2')
+    snp_iterator =  tsvreader.geno_iter(options.geno)
+    n = len(y)
+    m = g.shape[1]
+    ps, ts = run_gwas('other',n,m,k,y,g.T)
+    print np.array(ps)
+    print len(ps),sum(ps)
+    # Test results
+    p1 = round(ps[0],4)
+    p2 = round(ps[-1],4)
+
+elif cmd == 'iterator':
     if options.remove_missing_phenotypes:
         raise Exception('Can not use --remove-missing-phenotypes with LMM2')
     snp_iterator =  tsvreader.geno_iter(options.geno)
-- 
cgit v1.2.3