From 7b47dd4a9ca1227a0df4f3e7bf2358633944de0f Mon Sep 17 00:00:00 2001
From: Pjotr Prins
Date: Sat, 14 Mar 2015 08:39:44 +0300
Subject: kinship: working on pylmm

---
 wqflask/wqflask/my_pylmm/pyLMM/kinship.py | 10 ++++++++--
 wqflask/wqflask/my_pylmm/pyLMM/lmm.py     |  5 +++--
 wqflask/wqflask/my_pylmm/pyLMM/runlmm.py  | 26 +++++++++++++++++++-------
 3 files changed, 30 insertions(+), 11 deletions(-)

diff --git a/wqflask/wqflask/my_pylmm/pyLMM/kinship.py b/wqflask/wqflask/my_pylmm/pyLMM/kinship.py
index dc2d717d..353784aa 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/kinship.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/kinship.py
@@ -58,6 +58,12 @@ def compute_matrixMult(job,W,q = None):
 def f_init(q):
     compute_matrixMult.q = q
 
+def kinship_full(G,options):
+    print G.shape
+    m = np.dot(G.T,G)
+    m = m/G.shape[0]
+    return m
+
 # Calculate the kinship matrix from G (SNPs as rows!), returns K
 #
 def kinship(G,options):
@@ -118,7 +124,7 @@ def kinship(G,options):
           # print j,K_j[:,0]
           K = K + K_j
 
-    print K.shape,K
+    print "kiship.kinship: ",K.shape,K
     K = K / float(snps)
     outFile = 'runtest.kin'
     if options.verbose: sys.stderr.write("Saving Kinship file to %s\n" % outFile)
@@ -130,7 +136,7 @@ def kinship(G,options):
        if options.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 48a931c5..36c3602f 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
@@ -326,7 +326,7 @@ 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, is_testing=False):
+def calculate_kinship(genotype_matrix, temp_data=None, is_testing=False):
     """
     genotype_matrix is an n x m matrix encoding SNP minor alleles.
     
@@ -361,7 +361,8 @@ def calculate_kinship(genotype_matrix, temp_data, is_testing=False):
         genotype_matrix[:,counter] = (genotype_matrix[:,counter] - values_mean) / np.sqrt(vr)
         
         percent_complete = int(round((counter/m)*45))
-        temp_data.store("percent_complete", percent_complete)
+        if temp_data != None:
+            temp_data.store("percent_complete", percent_complete)
         
     genotype_matrix = genotype_matrix[:,keep]
     print("genotype_matrix: ", pf(genotype_matrix))
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
index 69a80411..627cc7a4 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
@@ -21,8 +21,8 @@ from optparse import OptionParser
 import sys
 import tsvreader
 import numpy as np
-from lmm import gn2_load_redis
-from kinship import kinship
+from lmm import gn2_load_redis, calculate_kinship
+from kinship import kinship, kinship_full
 from genotype import normalizeGenotype
 from phenotype import removeMissingPhenotypes
 
@@ -106,10 +106,22 @@ if cmd == 'redis':
     assert(options.testing and round(ps[-1],4)==0.3461)
 elif cmd == 'kinship':
     gn = []
-    for snp in G.T:
-        gn.append( normalizeGenotype(snp) )
-    G2 = np.array(gn)
-    print G2.shape,G2
-    K = kinship(G2,options)
+    for ind_g in g:
+        if len(gn)>=8000: break
+        gn.append( normalizeGenotype(ind_g) )
+    K = kinship_full(np.array(gn),options)
+    print "first Kinship method",K.shape,K
+    K = kinship(np.array(gn),options)
+    print "second Kinship method",K.shape,K
+    sys.exit(1)
+    gnt = np.array(gn).T
+    Y,g = removeMissingPhenotypes(y,gnt,options.verbose)
+    G = g
+    print G.shape,G
+    K = calculate_kinship(np.copy(G),None,options)
+    print G.shape,G
+    print "first Kinship method",K.shape,K
+    K = kinship(G.T,options)
+    assert(K[0][0]==1.28)
 else:
     print "Doing nothing"
-- 
cgit v1.2.3