aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--wqflask/wqflask/marker_regression/gemma_mapping.py46
1 files changed, 8 insertions, 38 deletions
diff --git a/wqflask/wqflask/marker_regression/gemma_mapping.py b/wqflask/wqflask/marker_regression/gemma_mapping.py
index 61e4897c..630a3afa 100644
--- a/wqflask/wqflask/marker_regression/gemma_mapping.py
+++ b/wqflask/wqflask/marker_regression/gemma_mapping.py
@@ -80,7 +80,7 @@ def run_gemma(this_trait, this_dataset, samples, vals, covariates, use_loco, maf
os.system(generate_k_command)
- gemma_command = GEMMA_WRAPPER_COMMAND + ' --json --input %s/gn2/%s.json -- ' % (TEMPDIR, k_output_filename) + GEMMAOPTS + ' -a %s/%s_snps.txt -lmm 2 -g %s/%s_geno.txt -p %s/gn2/%s.txt' % (flat_files('genotype/bimbam'),
+ gemma_command = GEMMA_WRAPPER_COMMAND + ' --json --input %s/gn2/%s.json -- ' % (TEMPDIR, k_output_filename) + GEMMAOPTS + ' -a %s/%s_snps.txt -lmm 9 -g %s/%s_geno.txt -p %s/gn2/%s.txt' % (flat_files('genotype/bimbam'),
genofile_name,
flat_files('genotype/bimbam'),
genofile_name,
@@ -101,7 +101,7 @@ def run_gemma(this_trait, this_dataset, samples, vals, covariates, use_loco, maf
marker_obs = parse_loco_output(this_dataset, gwa_output_filename)
return marker_obs, gwa_output_filename
else:
- marker_obs = parse_loco_output(this_dataset, gwa_output_filename)
+ marker_obs = parse_loco_output(this_dataset, gwa_output_filename, use_loco)
return marker_obs, gwa_output_filename
def gen_pheno_txt_file(this_dataset, genofile_name, vals, trait_filename):
@@ -145,41 +145,7 @@ def gen_covariates_file(this_dataset, covariates, samples):
outfile.write(str(this_covariate[i]) + "\t")
outfile.write("\n")
-def parse_gemma_output(genofile_name):
- included_markers = []
- p_values = []
- marker_obs = []
-
- with open("{}{}_output.assoc.txt".format(webqtlConfig.GENERATED_IMAGE_DIR, genofile_name)) as output_file:
- for line in output_file:
- if line.startswith("chr\t"):
- continue
- else:
- marker = {}
- marker['name'] = line.split("\t")[1]
- if line.split("\t")[0] != "X" and line.split("\t")[0] != "X/Y":
- if "chr" in line.split("\t")[0]:
- marker['chr'] = int(line.split("\t")[0][3:])
- else:
- marker['chr'] = int(line.split("\t")[0])
- else:
- marker['chr'] = line.split("\t")[0]
- marker['Mb'] = float(line.split("\t")[2]) / 1000000
- marker['p_value'] = float(line.split("\t")[9])
- if math.isnan(marker['p_value']) or (marker['p_value'] <= 0):
- marker['lod_score'] = 0
- #marker['lrs_value'] = 0
- else:
- marker['lod_score'] = -math.log10(marker['p_value'])
- #marker['lrs_value'] = -math.log10(marker['p_value']) * 4.61
- marker_obs.append(marker)
-
- included_markers.append(line.split("\t")[1])
- p_values.append(float(line.split("\t")[9]))
-
- return marker_obs
-
-def parse_loco_output(this_dataset, gwa_output_filename):
+def parse_loco_output(this_dataset, gwa_output_filename, loco="True"):
output_filelist = []
with open("{}/gn2/".format(TEMPDIR) + gwa_output_filename + ".json") as data_file:
@@ -218,7 +184,11 @@ def parse_loco_output(this_dataset, gwa_output_filename):
else:
marker['chr'] = line.split("\t")[0]
marker['Mb'] = float(line.split("\t")[2]) / 1000000
- marker['p_value'] = float(line.split("\t")[9])
+ if loco == "True":
+ marker['p_value'] = float(line.split("\t")[9])
+ else:
+ marker['p_value'] = float(line.split("\t")[10])
+ marker['additive'] = float(line.split("\t")[7])
if math.isnan(marker['p_value']) or (marker['p_value'] <= 0):
marker['lod_score'] = 0
#marker['lrs_value'] = 0