diff options
author | Pjotr Prins | 2015-03-11 12:07:54 +0300 |
---|---|---|
committer | Pjotr Prins | 2015-03-11 12:07:54 +0300 |
commit | a5d0d42f6977c4473c28f12999a62ad360750041 (patch) | |
tree | f637882f89d5007acc5b50179edbefa9963bcb11 /wqflask | |
parent | 6c165a4b9ec848eae0baeb4e4e3382bde4e9c1be (diff) | |
download | genenetwork2-a5d0d42f6977c4473c28f12999a62ad360750041.tar.gz |
& Added testing option, now the results match
Diffstat (limited to 'wqflask')
-rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/lmm.py | 16 | ||||
-rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/runlmm.py | 11 |
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) |