aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorPjotr Prins2015-03-11 12:07:54 +0300
committerPjotr Prins2015-03-11 12:07:54 +0300
commita5d0d42f6977c4473c28f12999a62ad360750041 (patch)
treef637882f89d5007acc5b50179edbefa9963bcb11
parent6c165a4b9ec848eae0baeb4e4e3382bde4e9c1be (diff)
downloadgenenetwork2-a5d0d42f6977c4473c28f12999a62ad360750041.tar.gz
& Added testing option, now the results match
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/lmm.py16
-rw-r--r--wqflask/wqflask/my_pylmm/pyLMM/runlmm.py11
2 files changed, 18 insertions, 9 deletions
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py b/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
index ab19bf08..5ffb9a4c 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/lmm.py
@@ -267,6 +267,7 @@ def run_other(pheno_vector,
"""
print("In run_other")
+ print("REML=",restricted_max_likelihood," REFIT=",refit, "TESTING=",is_testing)
with Bench("Calculate Kinship"):
kinship_matrix = calculate_kinship(genotype_matrix, tempdata, is_testing)
@@ -280,6 +281,7 @@ def run_other(pheno_vector,
lmm_ob.fit()
print("genotype_matrix: ", genotype_matrix.shape)
+ print(genotype_matrix)
with Bench("Doing GWAS"):
t_stats, p_values = GWAS(pheno_vector,
@@ -720,7 +722,7 @@ class LMM:
pl.title(title)
-def gn2_redis(key,species):
+def gn2_redis(key,species,is_testing=False):
json_params = Redis.get(key)
params = json.loads(json_params)
@@ -729,7 +731,8 @@ def gn2_redis(key,species):
print('kinship', np.array(params['kinship_matrix']))
print('pheno', np.array(params['pheno_vector']))
- print('geno', np.array(params['genotype_matrix']))
+ geno = np.array(params['genotype_matrix'])
+ print('geno', geno.shape, geno)
# sys.exit(1)
if species == "human" :
@@ -741,10 +744,11 @@ def gn2_redis(key,species):
tempdata = tempdata)
else:
ps, ts = run_other(pheno_vector = np.array(params['pheno_vector']),
- genotype_matrix = np.array(params['genotype_matrix']),
+ genotype_matrix = geno,
restricted_max_likelihood = params['restricted_max_likelihood'],
refit = params['refit'],
- tempdata = tempdata, is_testing=False)
+ tempdata = tempdata,
+ is_testing=is_testing)
results_key = "pylmm:results:" + params['temp_uuid']
@@ -769,7 +773,7 @@ def gn2_main():
gn2_redis(key,species)
-def gn2_load_redis(key,species,kinship,pheno,geno):
+def gn2_load_redis(key,species,kinship,pheno,geno,is_testing):
print("Loading Redis from parsed data")
params = dict(pheno_vector = pheno.tolist(),
genotype_matrix = geno.tolist(),
@@ -786,7 +790,7 @@ def gn2_load_redis(key,species,kinship,pheno,geno):
Redis.set(key, json_params)
Redis.expire(key, 60*60)
- return gn2_redis(key,species)
+ return gn2_redis(key,species,is_testing)
if __name__ == '__main__':
print("WARNING: Calling pylmm from lmm.py will become OBSOLETE, use runlmm.py instead!")
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
index 907a6835..738268be 100644
--- a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
+++ b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py
@@ -50,6 +50,9 @@ parser.add_option("--geno",dest="geno",
parser.add_option("-q", "--quiet",
action="store_false", dest="verbose", default=True,
help="don't print status messages to stdout")
+parser.add_option("--test",
+ action="store_true", dest="testing", default=False,
+ help="Testing mode")
(options, args) = parser.parse_args()
@@ -76,9 +79,9 @@ def normalizeGenotype(G):
# Run for every SNP list (num individuals)
x = True - np.isnan(G) # Matrix of True/False
m = G[x].mean() # Global mean value
- print m
+ # print m
s = np.sqrt(G[x].var()) # Global stddev
- print s
+ # print s
G[np.isnan(G)] = m # Plug-in mean values for missing
if s == 0:
G = G - m # Subtract the mean
@@ -92,6 +95,7 @@ gn = []
for snp in g:
gn.append( normalizeGenotype(snp) )
+gn = g
gn = np.array(gn)
print("After1",gn)
gnT = gn.T
@@ -99,6 +103,7 @@ print("After",gnT)
# G = gnT
G = gnT
print "G shape",G.shape
+# sys.exit(1)
# assert(G[0,0]==-2.25726341)
# Remove individuals with missing phenotypes
@@ -111,5 +116,5 @@ if v.sum():
k = k[keep,:][:,keep]
if cmd == 'redis':
- ps, ts = gn2_load_redis('testrun','other',k,y,G)
+ ps, ts = gn2_load_redis('testrun','other',np.array(k),y,G,options.testing)
print np.array(ps)