diff options
| -rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/lmm.py | 9 | ||||
| -rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/runlmm.py | 25 |
2 files changed, 32 insertions, 2 deletions
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py b/wqflask/wqflask/my_pylmm/pyLMM/lmm.py index 8c6e30e8..00bbf144 100644 --- a/wqflask/wqflask/my_pylmm/pyLMM/lmm.py +++ b/wqflask/wqflask/my_pylmm/pyLMM/lmm.py @@ -726,6 +726,11 @@ def gn2_redis(key,species): params = json.loads(json_params) tempdata = temp_data.TempData(params['temp_uuid']) + + print('kinship', np.array(params['kinship_matrix'][0:10][0:10])) + print('pheno', params['pheno_vector'][0:10]) + print('geno', params['genotype_matrix'][0:10][0:10]) + if species == "human" : ps, ts = run_human(pheno_vector = np.array(params['pheno_vector']), covariate_matrix = np.array(params['covariate_matrix']), @@ -748,6 +753,7 @@ def gn2_redis(key,species): #Pushing json_results into a list where it is the only item because blpop needs a list Redis.rpush(results_key, json_results) Redis.expire(results_key, 60*60) + return ps, ts # This is the main function used by Genenetwork2 (with environment) def gn2_main(): @@ -766,6 +772,7 @@ def gn2_load_redis(key,species,kinship,pheno,geno): print("Loading Redis from parsed data") params = dict(pheno_vector = pheno.tolist(), genotype_matrix = geno.tolist(), + kinship_matrix= kinship.tolist(), restricted_max_likelihood = True, refit = False, temp_uuid = "testrun_temp_uuid", @@ -778,7 +785,7 @@ def gn2_load_redis(key,species,kinship,pheno,geno): Redis.set(key, json_params) Redis.expire(key, 60*60) - gn2_redis(key,species) + return gn2_redis(key,species) 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 525c43fc..6db1bdbc 100644 --- a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py +++ b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py @@ -18,7 +18,9 @@ # along with this program. If not, see <http://www.gnu.org/licenses/>. from optparse import OptionParser +import sys import tsvreader +import numpy as np from lmm import gn2_load_redis usage = """ @@ -70,5 +72,26 @@ if options.geno: g = tsvreader.geno(options.geno) print g.shape +def normalizeGenotype(G): + x = True - np.isnan(G) + m = G[x].mean() + s = np.sqrt(G[x].var()) + G[np.isnan(G)] = m + if s == 0: G = G - m + else: G = (G - m) / s + return G + +gT = normalizeGenotype(g.T) + +# Remove individuals with missing phenotypes +v = np.isnan(y) +keep = True - v +if v.sum(): + if options.verbose: sys.stderr.write("Cleaning the phenotype vector by removing %d individuals...\n" % (v.sum())) + y = y[keep] + gT = gT[keep,:] + k = k[keep,:][:,keep] + if cmd == 'redis': - gn2_load_redis('testrun','other',k,y,g.T) + ps, ts = gn2_load_redis('testrun','other',k,y,gT) + print ps[0:10],ps[-10:] |