From 6a0a3626baad96deb1e8dc7d27fe1fa15e8c5b98 Mon Sep 17 00:00:00 2001 From: zsloan Date: Mon, 3 Oct 2016 19:16:56 +0000 Subject: Fixed problem that caused R/qtl with permutations to not work correctly --- .../wqflask/marker_regression/marker_regression.py | 5 +++- wqflask/wqflask/marker_regression/rqtl_mapping.py | 35 ++++++++++++---------- 2 files changed, 23 insertions(+), 17 deletions(-) diff --git a/wqflask/wqflask/marker_regression/marker_regression.py b/wqflask/wqflask/marker_regression/marker_regression.py index bb964961..f2a2eb8c 100644 --- a/wqflask/wqflask/marker_regression/marker_regression.py +++ b/wqflask/wqflask/marker_regression/marker_regression.py @@ -166,7 +166,10 @@ class MarkerRegression(object): self.model = start_vars['mapmodel_rqtl_geno'] if start_vars['pair_scan'] == "true": self.pair_scan = True - results = rqtl_mapping.run_rqtl_geno(self.vals, self.dataset, self.method, self.model, self.permCheck, self.num_perm, self.do_control, self.control_marker, self.manhattan_plot, self.pair_scan) + if self.permCheck and self.num_perm > 0: + perm_output, self.suggestive, self.significant, results = rqtl_mapping.run_rqtl_geno(self.vals, self.dataset, self.method, self.model, self.permCheck, self.num_perm, self.do_control, self.control_marker, self.manhattan_plot, self.pair_scan) + else: + results = rqtl_mapping.run_rqtl_geno(self.vals, self.dataset, self.method, self.model, self.permCheck, self.num_perm, self.do_control, self.control_marker, self.manhattan_plot, self.pair_scan) elif self.mapping_method == "reaper": if "startMb" in start_vars: #ZS: Check if first time page loaded, so it can default to ON if "additiveCheck" in start_vars: diff --git a/wqflask/wqflask/marker_regression/rqtl_mapping.py b/wqflask/wqflask/marker_regression/rqtl_mapping.py index 53ea6053..93bf717c 100644 --- a/wqflask/wqflask/marker_regression/rqtl_mapping.py +++ b/wqflask/wqflask/marker_regression/rqtl_mapping.py @@ -1,4 +1,5 @@ import rpy2.robjects as ro +import numpy as np from base.webqtlConfig import TMPDIR from utility import webqtlUtil @@ -72,9 +73,10 @@ def run_rqtl_geno(vals, dataset, method, model, permCheck, num_perm, do_control, else: perm_data_frame = scanone(cross_object, pheno_col = "the_pheno", n_perm = num_perm, model=model, method=method) - process_rqtl_perm_results(num_perm, perm_data_frame) # Functions that sets the thresholds for the webinterface - - return process_rqtl_results(result_data_frame) + perm_output, suggestive, significant = process_rqtl_perm_results(num_perm, perm_data_frame) # Functions that sets the thresholds for the webinterface + return perm_output, suggestive, significant, process_rqtl_results(result_data_frame) + else: + return process_rqtl_results(result_data_frame) def geno_to_rqtl_function(dataset): # TODO: Need to figure out why some genofiles have the wrong format and don't convert properly @@ -161,6 +163,19 @@ def process_pair_scan_results(result): return pair_scan_results +def process_rqtl_perm_results(num_perm, results): + perm_vals = [] + for line in str(results).split("\n")[1:(num_perm+1)]: + #print("R/qtl permutation line:", line.split()) + perm_vals.append(float(line.split()[1])) + + perm_output = perm_vals + suggestive = np.percentile(np.array(perm_vals), 67) + significant = np.percentile(np.array(perm_vals), 95) + print("SIGNIFICANT:", significant) + + return perm_output, suggestive, significant + def process_rqtl_results(result): # TODO: how to make this a one liner and not copy the stuff in a loop qtl_results = [] @@ -175,16 +190,4 @@ def process_rqtl_results(result): # TODO: how to make this a one liner an marker['lod_score'] = output[i][2] qtl_results.append(marker) - return qtl_results - -def process_rqtl_perm_results(num_perm, results): - perm_vals = [] - for line in str(results).split("\n")[1:(num_perm+1)]: - #print("R/qtl permutation line:", line.split()) - perm_vals.append(float(line.split()[1])) - - perm_output = perm_vals - suggestive = np.percentile(np.array(perm_vals), 67) - significant = np.percentile(np.array(perm_vals), 95) - - return perm_output, suggestive, significant \ No newline at end of file + return qtl_results \ No newline at end of file -- cgit v1.2.3