diff options
| -rw-r--r-- | wqflask/wqflask/my_pylmm/pyLMM/runlmm.py | 74 |
1 files changed, 30 insertions, 44 deletions
diff --git a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py index f095bb73..2d02e195 100644 --- a/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py +++ b/wqflask/wqflask/my_pylmm/pyLMM/runlmm.py @@ -27,6 +27,8 @@ import genotype import phenotype from standalone import uses +progress,mprint,debug,info,fatal = uses('progress','mprint','debug','info','fatal') + usage = """ python runlmm.py [options] command @@ -103,6 +105,29 @@ if options.geno and cmd != 'iterator': g = tsvreader.geno(options.geno) print g.shape +def check_results(ps,ts): + print np.array(ps) + print len(ps),sum(ps) + # Test results + p1 = round(ps[0],4) + p2 = round(ps[-1],4) + # sys.stderr.write(options.geno+"\n") + if options.geno == 'data/small.geno': + info("Validating results for "+options.geno) + assert p1==0.0708, "p1=%f" % p1 + assert p2==0.1417, "p2=%f" % p2 + if options.geno == 'data/small_na.geno': + info("Validating results for "+options.geno) + assert p1==0.0897, "p1=%f" % p1 + assert p2==0.0405, "p2=%f" % p2 + if options.geno == 'data/test8000.geno': + info("Validating results for "+options.geno) + # assert p1==0.8984, "p1=%f" % p1 + # assert p2==0.9621, "p2=%f" % p2 + assert round(sum(ps)) == 4070 + assert len(ps) == 8000 + info("Run completed") + if cmd == 'run': if options.remove_missing_phenotypes: raise Exception('Can not use --remove-missing-phenotypes with LMM2') @@ -110,22 +135,13 @@ if cmd == 'run': n = len(y) m = g.shape[1] ps, ts = run_gwas('other',n,m,k,y,g.T) - print np.array(ps) - print len(ps),sum(ps) - # Test results - p1 = round(ps[0],4) - p2 = round(ps[-1],4) - + check_results(ps,ts) elif cmd == 'iterator': if options.remove_missing_phenotypes: raise Exception('Can not use --remove-missing-phenotypes with LMM2') snp_iterator = tsvreader.geno_iter(options.geno) ps, ts = gn2_load_redis_iter('testrun_iter','other',k,y,snp_iterator) - print np.array(ps) - print len(ps),sum(ps) - # Test results - p1 = round(ps[0],4) - p2 = round(ps[-1],4) + check_results(ps,ts) elif cmd == 'redis_new': # The main difference between redis_new and redis is that missing # phenotypes are handled by the first @@ -138,23 +154,7 @@ elif cmd == 'redis_new': gt = G.T G = None ps, ts = gn2_load_redis('testrun','other',k,Y,gt,new_code=True) - print np.array(ps) - print len(ps),sum(ps) - # Test results - p1 = round(ps[0],4) - p2 = round(ps[-1],4) - sys.stderr.write(options.geno+"\n") - if options.geno == 'data/small.geno': - assert p1==0.0708, "p1=%f" % p1 - assert p2==0.1417, "p2=%f" % p2 - if options.geno == 'data/small_na.geno': - assert p1==0.0897, "p1=%f" % p1 - assert p2==0.0405, "p2=%f" % p2 - if options.geno == 'data/test8000.geno': - # assert p1==0.8984, "p1=%f" % p1 - # assert p2==0.9621, "p2=%f" % p2 - assert round(sum(ps)) == 4070 - assert len(ps) == 8000 + check_results(ps,ts) elif cmd == 'redis': # Emulating the redis setup of GN2 G = g @@ -177,21 +177,7 @@ elif cmd == 'redis': gt = G.T G = None ps, ts = gn2_load_redis('testrun','other',k,Y,gt, new_code=False) - print np.array(ps) - print len(ps),sum(ps) - # Test results 4070.02346579 - p1 = round(ps[0],4) - p2 = round(ps[-1],4) - sys.stderr.write(options.geno+"\n") - if options.geno == 'data/small.geno': - assert p1==0.0708, "p1=%f" % p1 - assert p2==0.1417, "p2=%f" % p2 - if options.geno == 'data/small_na.geno': - assert p1==0.0897, "p1=%f" % p1 - assert p2==0.0405, "p2=%f" % p2 - if options.geno == 'data/test8000.geno': - assert int(sum(ps)) == 4070 - assert len(ps) == 8000 + check_results(ps,ts) elif cmd == 'kinship': G = g print "Original G",G.shape, "\n", G @@ -235,4 +221,4 @@ elif cmd == 'kinship': assert k3==1.4352, "k3=%f" % k3 else: - print "Doing nothing" + fatal("Doing nothing") |