diff options
| -rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/kinship.py | 159 | ||||
| -rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/runlmm.py | 13 |
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" |