Merge branch 'master' of github.com:pjotrp/genenetwork2

3 files changed, 35 insertions, 24 deletions

diff --git a/wqflask/wqflask/my_pylmm/pyLMM/kinship.py b/wqflask/wqflask/my_pylmm/pyLMM/kinship.py
index 9ab48510..28f2042d 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/kinship.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/kinship.py
@@ -73,12 +73,11 @@ def f_init(q):
 
 # Calculate the kinship matrix from G (SNPs as rows!), returns K
 #
-def kinship(G,options):
+def kinship(G,computeSize=1000,numThreads=None,useBLAS=False,verbose=True):
    numThreads = None
-   if options.numThreads:
-      numThreads = int(options.numThreads)
-   options.computeSize = 1000
-   matrix_initialize(options.useBLAS)
+   if numThreads:
+      numThreads = int(numThreads)
+   matrix_initialize(useBLAS)
    
    sys.stderr.write(str(G.shape)+"\n")
    n = G.shape[1] # inds
@@ -92,9 +91,9 @@ def kinship(G,options):
    p = mp.Pool(numThreads, f_init, [q])
    cpu_num = mp.cpu_count()
    print "CPU cores:",cpu_num
-   iterations = snps/options.computeSize+1
-   if options.testing:
-      iterations = 8
+   iterations = snps/computeSize+1
+   # if testing:
+   #    iterations = 8
    # jobs = range(0,8) # range(0,iterations)
 
    results = []
@@ -103,14 +102,14 @@ def kinship(G,options):
 
    completed = 0
    for job in range(iterations):
-      if options.verbose:
-         sys.stderr.write("Processing job %d first %d SNPs\n" % (job, ((job+1)*options.computeSize)))
-      W = compute_W(job,G,n,snps,options.computeSize)
+      if verbose:
+         sys.stderr.write("Processing job %d first %d SNPs\n" % (job, ((job+1)*computeSize)))
+      W = compute_W(job,G,n,snps,computeSize)
       if numThreads == 1:
          # Single-core
          compute_matrixMult(job,W,q)
          j,x = q.get()
-         if options.verbose: sys.stderr.write("Job "+str(j)+" finished\n")
+         if verbose: sys.stderr.write("Job "+str(j)+" finished\n")
          K_j = x
          # print j,K_j[:,0]
          K = K + K_j
@@ -122,7 +121,7 @@ def kinship(G,options):
             time.sleep(0.1)
             try: 
                j,x = q.get_nowait()
-               if options.verbose: sys.stderr.write("Job "+str(j)+" finished\n")
+               if verbose: sys.stderr.write("Job "+str(j)+" finished\n")
                K_j = x
                # print j,K_j[:,0]
                K = K + K_j
@@ -134,7 +133,7 @@ def kinship(G,options):
       # results contains the growing result set
       for job in range(len(results)-completed):
          j,x = q.get(True,15)
-         if options.verbose: sys.stderr.write("Job "+str(j)+" finished\n")
+         if verbose: sys.stderr.write("Job "+str(j)+" finished\n")
          K_j = x
          # print j,K_j[:,0]
          K = K + K_j
@@ -143,13 +142,13 @@ def kinship(G,options):
    K = K / float(snps)
    
    # outFile = 'runtest.kin'
-   # if options.verbose: sys.stderr.write("Saving Kinship file to %s\n" % outFile)
+   # if verbose: sys.stderr.write("Saving Kinship file to %s\n" % outFile)
    # np.savetxt(outFile,K)
 
-   # if options.saveKvaKve:
-   #    if options.verbose: sys.stderr.write("Obtaining Eigendecomposition\n")
+   # if saveKvaKve:
+   #    if verbose: sys.stderr.write("Obtaining Eigendecomposition\n")
    #    Kva,Kve = linalg.eigh(K)
-   #    if options.verbose: sys.stderr.write("Saving eigendecomposition to %s.[kva | kve]\n" % outFile)
+   #    if verbose: sys.stderr.write("Saving eigendecomposition to %s.[kva | kve]\n" % outFile)
    #    np.savetxt(outFile+".kva",Kva)
    #    np.savetxt(outFile+".kve",Kve)
    return K      
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py b/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
index a9744e72..1383d3fc 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
@@ -49,6 +49,8 @@ has_gn2=True
 
 from utility.benchmark import Bench
 from utility import temp_data
+from kinship import kinship, kinship_full
+import genotype
 
 try:
     from wqflask.my_pylmm.pyLMM import chunks
@@ -270,7 +272,7 @@ def run_other(pheno_vector,
     print("In run_other")
     print("REML=",restricted_max_likelihood," REFIT=",refit)
     with Bench("Calculate Kinship"):
-        kinship_matrix = calculate_kinship(genotype_matrix, tempdata)
+        kinship_matrix,genotype_matrix = calculate_kinship(genotype_matrix, tempdata)
     
     print("kinship_matrix: ", pf(kinship_matrix))
     print("kinship_matrix.shape: ", pf(kinship_matrix.shape))
@@ -326,7 +328,15 @@ def matrixMult(A,B):
     return linalg.fblas.dgemm(alpha=1.,a=AA,b=BB,trans_a=transA,trans_b=transB)
 
 
-def calculate_kinship(genotype_matrix, temp_data=None):
+def calculate_kinship_new(genotype_matrix, temp_data=None):
+    """ 
+    Call the new kinship calculation where genotype_matrix contains
+    inds (columns) by snps (rows).
+    """
+    G = np.apply_along_axis( genotype.normalize, axis=0, arr=genotype_matrix)
+    return kinship(G.T),G
+
+def calculate_kinship_old(genotype_matrix, temp_data=None):
     """
     genotype_matrix is an n x m matrix encoding SNP minor alleles.
     
@@ -366,7 +376,9 @@ def calculate_kinship(genotype_matrix, temp_data=None):
     genotype_matrix = genotype_matrix[:,keep]
     print("genotype_matrix: ", pf(genotype_matrix))
     kinship_matrix = np.dot(genotype_matrix, genotype_matrix.T) * 1.0/float(m)
-    return kinship_matrix
+    return kinship_matrix,genotype_matrix
+
+calculate_kinship = calculate_kinship_new  # alias
 
 def GWAS(pheno_vector,
          genotype_matrix,
@@ -735,8 +747,6 @@ def gn2_redis(key,species):
 
     
     print('pheno', np.array(params['pheno_vector']))
-
-    # sys.exit(1)
     
     if species == "human" :
         print('kinship', np.array(params['kinship_matrix']))
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
index 469ba6c9..6bb79856 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
@@ -116,7 +116,9 @@ if cmd == 'redis':
     g = None
     gnt = None
 
-    ps, ts = gn2_load_redis('testrun','other',k,Y,G.T)
+    gt = G.T
+    G = None
+    ps, ts = gn2_load_redis('testrun','other',k,Y,gt)
     print np.array(ps)
     # Test results
     p1 = round(ps[0],4)