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(-)
(limited to 'wqflask')
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