From ece60a2d25ed71980121ac1c1854c3ab3514f453 Mon Sep 17 00:00:00 2001 From: zsloan Date: Wed, 24 Feb 2016 23:36:12 +0000 Subject: GEMMA mostly working with gn1 mapping, though there may be some issue with the y scale --- wqflask/wqflask/marker_regression/gemma_mapping.py | 13 +++++++++---- wqflask/wqflask/marker_regression/marker_regression.py | 9 +++++---- 2 files changed, 14 insertions(+), 8 deletions(-) (limited to 'wqflask') diff --git a/wqflask/wqflask/marker_regression/gemma_mapping.py b/wqflask/wqflask/marker_regression/gemma_mapping.py index b1ab780c..997b692d 100644 --- a/wqflask/wqflask/marker_regression/gemma_mapping.py +++ b/wqflask/wqflask/marker_regression/gemma_mapping.py @@ -1,6 +1,9 @@ import os from base import webqtlConfig +from utility.tools import gemma_command + +GEMMA_PATH,GEMMA_COMMAND = gemma_command() def run_gemma(this_dataset, samples, vals): """Generates p-values for each marker using GEMMA""" @@ -9,9 +12,11 @@ def run_gemma(this_dataset, samples, vals): gen_pheno_txt_file(this_dataset, samples, vals) - os.chdir("{}gemma".format(webqtlConfig.HTMLPATH)) + os.chdir(GEMMA_PATH) - gemma_command = './gemma -bfile %s -k output_%s.cXX.txt -lmm 1 -o output/%s_output' % (this_dataset.group.name, + gemma_command = GEMMA_COMMAND + ' -bfile %s/%s -k %s/output/%s.cXX.txt -lmm 1 -o %s_output' % (GEMMA_PATH, + this_dataset.group.name, + GEMMA_PATH, this_dataset.group.name, this_dataset.group.name) print("gemma_command:" + gemma_command) @@ -25,14 +30,14 @@ def run_gemma(this_dataset, samples, vals): def gen_pheno_txt_file(this_dataset, samples, vals): """Generates phenotype file for GEMMA""" - with open("{}gemma/{}.fam".format(webqtlConfig.HTMLPATH, this_dataset.group.name), "w") as outfile: + with open("{}/{}.fam".format(GEMMA_PATH, this_dataset.group.name), "w") as outfile: for i, sample in enumerate(samples): outfile.write(str(sample) + " " + str(sample) + " 0 0 0 " + str(vals[i]) + "\n") def parse_gemma_output(this_dataset): included_markers = [] p_values = [] - with open("{}gemma/output/{}_output.assoc.txt".format(webqtlConfig.HTMLPATH, this_dataset.group.name)) as output_file: + with open("{}/output/{}_output.assoc.txt".format(GEMMA_PATH, this_dataset.group.name)) as output_file: for line in output_file: if line.startswith("chr"): continue diff --git a/wqflask/wqflask/marker_regression/marker_regression.py b/wqflask/wqflask/marker_regression/marker_regression.py index c16c885e..5377201b 100755 --- a/wqflask/wqflask/marker_regression/marker_regression.py +++ b/wqflask/wqflask/marker_regression/marker_regression.py @@ -39,10 +39,11 @@ from utility import helper_functions from utility import Plot, Bunch from utility import temp_data from utility.benchmark import Bench -from utility.tools import pylmm_command, plink_command +from utility.tools import pylmm_command, plink_command, gemma_command PYLMM_PATH,PYLMM_COMMAND = pylmm_command() PLINK_PATH,PLINK_COMMAND = plink_command() +GEMMA_PATH,GEMMA_COMMAND = gemma_command() class MarkerRegression(object): @@ -86,7 +87,7 @@ class MarkerRegression(object): self.dataset.group.get_markers() if self.mapping_method == "gemma": - self.score_type = "LOD" + self.score_type = "LRS" included_markers, p_values = gemma_mapping.run_gemma(self.dataset, self.samples, self.vals) self.dataset.group.get_specified_markers(markers = included_markers) self.dataset.group.markers.add_pvalues(p_values) @@ -219,9 +220,9 @@ class MarkerRegression(object): self.gen_pheno_txt_file() - os.chdir("/home/zas1024/gene/web/gemma") + #os.chdir("/home/zas1024/gene/web/gemma") - gemma_command = './gemma -bfile %s -k output_%s.cXX.txt -lmm 1 -o %s_output' % ( + gemma_command = GEMMA_COMMAND + ' -bfile %s -k output_%s.cXX.txt -lmm 1 -o %s_output' % ( self.dataset.group.name, self.dataset.group.name, self.dataset.group.name) -- cgit v1.2.3