aboutsummaryrefslogtreecommitdiff
path: root/wqflask
diff options
context:
space:
mode:
Diffstat (limited to 'wqflask')
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/kinship.py159
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/runlmm.py13
2 files changed, 89 insertions, 83 deletions
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/kinship.py b/wqflask/wqflask/my_pylmm/pyLMM/kinship.py
index 1dca1c57..8caa753b 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/kinship.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/kinship.py
@@ -57,91 +57,94 @@ def compute_matrixMult(job,W,q = None):
return job
def f_init(q):
- compute_matrixMult.q = q
+ compute_matrixMult.q = q
def kinship_full(G):
- print G.shape
- m = G.shape[0] # snps
- n = G.shape[1] # inds
- sys.stderr.write(str(m)+" SNPs\n")
- assert m>n, "n should be larger than m (snps>inds)"
- m = np.dot(G.T,G)
- m = m/G.shape[0]
- return m
+ """
+ Calculate the Kinship matrix using a full dot multiplication
+ """
+ print G.shape
+ m = G.shape[0] # snps
+ n = G.shape[1] # inds
+ sys.stderr.write(str(m)+" SNPs\n")
+ assert m>n, "n should be larger than m (snps>inds)"
+ 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):
- numThreads = None
- if options.numThreads:
- numThreads = int(options.numThreads)
- options.computeSize = 1000
- matrix_initialize(options.useBLAS)
-
- sys.stderr.write(str(G.shape)+"\n")
- n = G.shape[1] # inds
- inds = n
- m = G.shape[0] # snps
- snps = m
- sys.stderr.write(str(m)+" SNPs\n")
- assert m>n, "n should be larger than m (snps>inds)"
-
- q = mp.Queue()
- p = mp.Pool(numThreads, f_init, [q])
- iterations = snps/options.computeSize+1
- if options.testing:
+ numThreads = None
+ if options.numThreads:
+ numThreads = int(options.numThreads)
+ options.computeSize = 1000
+ matrix_initialize(options.useBLAS)
+
+ sys.stderr.write(str(G.shape)+"\n")
+ n = G.shape[1] # inds
+ inds = n
+ m = G.shape[0] # snps
+ snps = m
+ sys.stderr.write(str(m)+" SNPs\n")
+ assert m>n, "n should be larger than m (snps>inds)"
+
+ q = mp.Queue()
+ p = mp.Pool(numThreads, f_init, [q])
+ iterations = snps/options.computeSize+1
+ if options.testing:
iterations = 8
- # jobs = range(0,8) # range(0,iterations)
-
- results = []
-
- K = np.zeros((n,n)) # The Kinship matrix has dimension individuals x individuals
-
- 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 numThreads == 1:
- compute_matrixMult(job,W,q)
- j,x = q.get()
- if options.verbose: sys.stderr.write("Job "+str(j)+" finished\n")
- K_j = x
- # print j,K_j[:,0]
- K = K + K_j
- else:
- results.append(p.apply_async(compute_matrixMult, (job,W)))
- # Do we have a result?
- try:
- j,x = q.get_nowait()
- if options.verbose: sys.stderr.write("Job "+str(j)+" finished\n")
- K_j = x
- # print j,K_j[:,0]
- K = K + K_j
- completed += 1
- except Queue.Empty:
- pass
-
- if numThreads == None or numThreads > 1:
- for job in range(len(results)-completed):
- j,x = q.get(True,15)
- if options.verbose: sys.stderr.write("Job "+str(j)+" finished\n")
- K_j = x
- # print j,K_j[:,0]
- K = K + K_j
-
- K = K / float(snps)
- outFile = 'runtest.kin'
- if options.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")
- Kva,Kve = linalg.eigh(K)
- 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
+ # jobs = range(0,8) # range(0,iterations)
+
+ results = []
+
+ K = np.zeros((n,n)) # The Kinship matrix has dimension individuals x individuals
+
+ 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 numThreads == 1:
+ compute_matrixMult(job,W,q)
+ j,x = q.get()
+ if options.verbose: sys.stderr.write("Job "+str(j)+" finished\n")
+ K_j = x
+ # print j,K_j[:,0]
+ K = K + K_j
+ else:
+ results.append(p.apply_async(compute_matrixMult, (job,W)))
+ # Do we have a result?
+ try:
+ j,x = q.get_nowait()
+ if options.verbose: sys.stderr.write("Job "+str(j)+" finished\n")
+ K_j = x
+ # print j,K_j[:,0]
+ K = K + K_j
+ completed += 1
+ except Queue.Empty:
+ pass
+
+ if numThreads == None or numThreads > 1:
+ for job in range(len(results)-completed):
+ j,x = q.get(True,15)
+ if options.verbose: sys.stderr.write("Job "+str(j)+" finished\n")
+ K_j = x
+ # print j,K_j[:,0]
+ K = K + K_j
+
+ K = K / float(snps)
+ outFile = 'runtest.kin'
+ if options.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")
+ Kva,Kve = linalg.eigh(K)
+ 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/runlmm.py b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
index 26485136..547ac7f1 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
@@ -140,7 +140,7 @@ elif cmd == 'kinship':
if not options.skip_genotype_normalization:
G = np.apply_along_axis( genotype.normalize, axis=1, arr=G)
- if True:
+ if False:
K = kinship_full(G)
print "Genotype",G.shape, "\n", G
print "first Kinship method",K.shape,"\n",K
@@ -157,11 +157,14 @@ elif cmd == 'kinship':
k3 = round(K3[0][0],4)
if options.geno == 'data/small.geno':
assert k1==0.7939, "k1=%f" % k1
- assert k2==0.7939, "k1=%f" % k1
- assert k3==0.7939, "k1=%f" % k1
+ assert k2==0.7939, "k2=%f" % k2
+ assert k3==0.7939, "k3=%f" % k3
if options.geno == 'data/small_na.geno':
assert k1==0.7172, "k1=%f" % k1
- assert k2==0.7172, "k1=%f" % k1
- assert k3==0.7172, "k1=%f" % k1
+ assert k2==0.7172, "k2=%f" % k2
+ assert k3==0.7172, "k3=%f" % k3
+ if options.geno == 'data/test8000.geno':
+ assert k3==1.4352, "k3=%f" % k3
+
else:
print "Doing nothing"