diff options
| -rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/kinship.py | 33 | ||||
| -rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/runlmm.py | 5 |
2 files changed, 20 insertions, 18 deletions
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/kinship.py b/wqflask/wqflask/my_pylmm/pyLMM/kinship.py index 8caa753b..a5497a77 100644 --- a/wqflask/wqflask/my_pylmm/pyLMM/kinship.py +++ b/wqflask/wqflask/my_pylmm/pyLMM/kinship.py @@ -26,20 +26,31 @@ import Queue from optmatrix import matrix_initialize, matrixMultT +def kinship_full(G): + """ + 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 def compute_W(job,G,n,snps,compute_size): """ Read 1000 SNPs at a time into matrix and return the result """ - W = np.ones((n,compute_size)) * np.nan # W matrix has dimensions individuals x SNPs (initially all NaNs) + m = compute_size + W = np.ones((n,m)) * np.nan # W matrix has dimensions individuals x SNPs (initially all NaNs) for j in range(0,compute_size): - row = job*compute_size + j - if row >= compute_size or row>=snps: + pos = job*m + j # real position + if pos >= snps: W = W[:,range(0,j)] break - # print job,compute_size,j snp = G[job*compute_size+j] - # print snp.shape,snp if snp.var() == 0: continue W[:,j] = snp # set row to list of SNPs @@ -59,18 +70,6 @@ def compute_matrixMult(job,W,q = None): def f_init(q): compute_matrixMult.q = q -def kinship_full(G): - """ - 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 # diff --git a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py index 547ac7f1..a322f988 100644 --- a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py +++ b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py @@ -70,6 +70,9 @@ parser.add_option("--saveKvaKve", parser.add_option("--test", action="store_true", dest="testing", default=False, help="Testing mode") +parser.add_option("--test-kinship", + action="store_true", dest="test_kinship", default=False, + help="Testing mode for Kinship calculation") (options, args) = parser.parse_args() @@ -140,7 +143,7 @@ elif cmd == 'kinship': if not options.skip_genotype_normalization: G = np.apply_along_axis( genotype.normalize, axis=1, arr=G) - if False: + if options.test_kinship: K = kinship_full(G) print "Genotype",G.shape, "\n", G print "first Kinship method",K.shape,"\n",K |