From 4773b58bd9fdf6b25bddd2dd75570d56e7d1325c Mon Sep 17 00:00:00 2001 From: zsloan Date: Mon, 29 Jul 2019 13:43:23 -0500 Subject: Added Christian's RUST qtlreaper as an option --- .../marker_regression/display_mapping_results.py | 3 + .../wqflask/marker_regression/qtlreaper_mapping.py | 116 ++++++++++++++++++++- wqflask/wqflask/marker_regression/run_mapping.py | 49 +++++---- .../templates/show_trait_mapping_tools.html | 17 ++- wqflask/wqflask/views.py | 3 +- 5 files changed, 160 insertions(+), 28 deletions(-) diff --git a/wqflask/wqflask/marker_regression/display_mapping_results.py b/wqflask/wqflask/marker_regression/display_mapping_results.py index d9601405..e6924f40 100644 --- a/wqflask/wqflask/marker_regression/display_mapping_results.py +++ b/wqflask/wqflask/marker_regression/display_mapping_results.py @@ -233,6 +233,9 @@ class DisplayMappingResults(object): if self.use_loco == "True": self.gwa_filename = start_vars['gwa_filename'] + if 'reaper_version' in start_vars.keys() and self.mapping_method == "reaper": + self.reaper_version = start_vars['reaper_version'] + self.selectedChr = int(start_vars['selected_chr']) self.strainlist = start_vars['samples'] diff --git a/wqflask/wqflask/marker_regression/qtlreaper_mapping.py b/wqflask/wqflask/marker_regression/qtlreaper_mapping.py index d58c59c8..b74a1c4d 100644 --- a/wqflask/wqflask/marker_regression/qtlreaper_mapping.py +++ b/wqflask/wqflask/marker_regression/qtlreaper_mapping.py @@ -1,7 +1,121 @@ +import os, math, string, random, json + +from base import webqtlConfig +from base.trait import GeneralTrait +from base.data_set import create_dataset +from utility.tools import flat_files, REAPER_COMMAND, TEMPDIR + import utility.logger logger = utility.logger.getLogger(__name__ ) -def gen_reaper_results(this_trait, dataset, samples_before, trait_vals, json_data, num_perm, bootCheck, num_bootstrap, do_control, control_marker, manhattan_plot): +def run_reaper(this_trait, this_dataset, samples, vals, json_data, num_perm, boot_check, num_bootstrap, do_control, control_marker, manhattan_plot): + """Generates p-values for each marker using qtlreaper""" + + if this_dataset.group.genofile != None: + genofile_name = this_dataset.group.genofile[:-5] + else: + genofile_name = this_dataset.group.name + + trait_filename = str(this_trait.name) + "_" + str(this_dataset.name) + "_pheno" + gen_pheno_txt_file(this_dataset, genofile_name, samples, vals, trait_filename) + + output_filename = this_dataset.group.name + "_GWA_" + ''.join(random.choice(string.ascii_uppercase + string.digits) for _ in range(6)) + bootstrap_filename = None + permu_filename = None + + opt_list = [] + if boot_check and num_bootstrap > 0: + bootstrap_filename = this_dataset.group.name + "_BOOTSTRAP_" + ''.join(random.choice(string.ascii_uppercase + string.digits) for _ in range(6)) + + opt_list.append("-b") + opt_list.append("--n_bootstrap " + str(num_bootstrap)) + opt_list.append("--bootstrap_output " + webqtlConfig.GENERATED_IMAGE_DIR + bootstrap_filename + ".txt") + if num_perm > 0: + permu_filename = this_dataset.group.name + "_PERM_" + ''.join(random.choice(string.ascii_uppercase + string.digits) for _ in range(6)) + opt_list.append("-n " + str(num_perm)) + opt_list.append("--permu_output " + webqtlConfig.GENERATED_IMAGE_DIR + permu_filename + ".txt") + if control_marker != "" and do_control == "true": + opt_list.append("-c " + control_marker) + + reaper_command = REAPER_COMMAND + ' --geno {0}/{1}.geno --traits {2}/gn2/{3}.txt {4} -o {5}{6}.txt'.format(flat_files('genotype'), + genofile_name, + TEMPDIR, + trait_filename, + " ".join(opt_list), + webqtlConfig.GENERATED_IMAGE_DIR, + output_filename) + + logger.debug("reaper_command:" + reaper_command) + os.system(reaper_command) + + marker_obs, permu_vals, bootstrap_vals = parse_reaper_output(output_filename, permu_filename, bootstrap_filename) + + suggestive = 0 + significant = 0 + if len(permu_vals) > 0: + suggestive = permu_vals[int(num_perm*0.37-1)] + significant = permu_vals[int(num_perm*0.95-1)] + + return marker_obs, permu_vals, suggestive, significant, bootstrap_vals + +def gen_pheno_txt_file(this_dataset, genofile_name, samples, vals, trait_filename): + """Generates phenotype file for GEMMA""" + + with open("{}/gn2/{}.txt".format(TEMPDIR, trait_filename), "w") as outfile: + outfile.write("Trait\t") + + filtered_sample_list = [] + filtered_vals_list = [] + for i, sample in enumerate(samples): + if vals[i] != "x": + filtered_sample_list.append(sample) + filtered_vals_list.append(vals[i]) + + samples_string = "\t".join(filtered_sample_list) + outfile.write(samples_string + "\n") + outfile.write("T1\t") + values_string = "\t".join(filtered_vals_list) + outfile.write(values_string) + +def parse_reaper_output(gwa_filename, permu_filename, bootstrap_filename): + included_markers = [] + p_values = [] + marker_obs = [] + + with open("{}{}.txt".format(webqtlConfig.GENERATED_IMAGE_DIR, gwa_filename)) as output_file: + for line in output_file: + if line.startswith("ID\t"): + continue + else: + marker = {} + marker['name'] = line.split("\t")[1] + try: + marker['chr'] = int(line.split("\t")[2]) + except: + marker['chr'] = line.split("\t")[2] + marker['cM'] = float(line.split("\t")[3]) + marker['Mb'] = float(line.split("\t")[4]) + marker['p_value'] = float(line.split("\t")[7]) + marker['lrs_value'] = float(line.split("\t")[5]) + marker['lod_score'] = marker['lrs_value'] / 4.61 + marker['additive'] = float(line.split("\t")[6]) + marker_obs.append(marker) + + permu_vals = [] + if permu_filename: + with open("{}{}.txt".format(webqtlConfig.GENERATED_IMAGE_DIR, permu_filename)) as permu_file: + for line in permu_file: + permu_vals.append(float(line)) + + bootstrap_vals = [] + if bootstrap_filename: + with open("{}{}.txt".format(webqtlConfig.GENERATED_IMAGE_DIR, bootstrap_filename)) as bootstrap_file: + for line in bootstrap_file: + bootstrap_vals.append(int(line)) + + return marker_obs, permu_vals, bootstrap_vals + +def run_original_reaper(this_trait, dataset, samples_before, trait_vals, json_data, num_perm, bootCheck, num_bootstrap, do_control, control_marker, manhattan_plot): genotype = dataset.group.read_genotype_file(use_reaper=True) if manhattan_plot != True: diff --git a/wqflask/wqflask/marker_regression/run_mapping.py b/wqflask/wqflask/marker_regression/run_mapping.py index 6e9fe85c..648d646c 100644 --- a/wqflask/wqflask/marker_regression/run_mapping.py +++ b/wqflask/wqflask/marker_regression/run_mapping.py @@ -239,33 +239,36 @@ class RunMapping(object): self.bootCheck = False self.num_bootstrap = 0 + self.reaper_version = start_vars['reaper_version'] + self.control_marker = start_vars['control_marker'] self.do_control = start_vars['do_control'] logger.info("Running qtlreaper") - results, self.perm_output, self.suggestive, self.significant, self.bootstrap_results = rust_reaper_mapping.run_reaper(self.this_trait, - self.dataset, - self.samples, - self.vals, - self.json_data, - self.num_perm, - self.bootCheck, - self.num_bootstrap, - self.do_control, - self.control_marker, - self.manhattan_plot) - - # results, self.json_data, self.perm_output, self.suggestive, self.significant, self.bootstrap_results = qtlreaper_mapping.gen_reaper_results(self.this_trait, - # self.dataset, - # self.samples, - # self.vals, - # self.json_data, - # self.num_perm, - # self.bootCheck, - # self.num_bootstrap, - # self.do_control, - # self.control_marker, - # self.manhattan_plot) + if self.reaper_version == "new": + results, self.perm_output, self.suggestive, self.significant, self.bootstrap_results = qtlreaper_mapping.run_reaper(self.this_trait, + self.dataset, + self.samples, + self.vals, + self.json_data, + self.num_perm, + self.bootCheck, + self.num_bootstrap, + self.do_control, + self.control_marker, + self.manhattan_plot) + else: + results, self.json_data, self.perm_output, self.suggestive, self.significant, self.bootstrap_results = qtlreaper_mapping.run_original_reaper(self.this_trait, + self.dataset, + self.samples, + self.vals, + self.json_data, + self.num_perm, + self.bootCheck, + self.num_bootstrap, + self.do_control, + self.control_marker, + self.manhattan_plot) elif self.mapping_method == "plink": self.score_type = "-log(p)" self.manhattan_plot = True diff --git a/wqflask/wqflask/templates/show_trait_mapping_tools.html b/wqflask/wqflask/templates/show_trait_mapping_tools.html index b5c37b8b..22707f41 100644 --- a/wqflask/wqflask/templates/show_trait_mapping_tools.html +++ b/wqflask/wqflask/templates/show_trait_mapping_tools.html @@ -115,6 +115,15 @@ {% if dataset.group.mapping_id == "1" %}
+
+ +
+ +
+
@@ -347,9 +356,11 @@
GEMMA
GEMMA is software implementing the Genome-wide Efficient Mixed Model Association algorithm for a standard linear mixed model for genome-wide association studies (GWAS).
{% if dataset.group.mapping_id == "1" %} -
Interval Mapping
-
Interval mapping is a process in which the statistical significance of a hypothetical QTL is evaluated at regular points across a chromosome, even in the absence of explicit genotype data at those points.
-
R/qtl
+
Interval Mapping
+
Interval mapping is a process in which the statistical significance of a hypothetical QTL is evaluated at regular points across a chromosome, even in the absence of explicit genotype data at those points.

+ The default option is a version of qtlrearper rewritten in RUST that should fix some of the issues the original had. You can still use the original with the Version option.
+ Please let us know if you encounter any issues while using the new version.
+
R/qtl
R/qtl is an extensible, interactive environment for mapping quantitative trait loci (QTL) in experimental crosses.
{% endif %} diff --git a/wqflask/wqflask/views.py b/wqflask/wqflask/views.py index 2ff6d4d7..fe858d52 100644 --- a/wqflask/wqflask/views.py +++ b/wqflask/wqflask/views.py @@ -628,7 +628,8 @@ def mapping_results_page(): 'haplotypeAnalystCheck', 'mapmethod_rqtl_geno', 'mapmodel_rqtl_geno', - 'temp_trait' + 'temp_trait', + 'reaper_version' ) start_vars = {} for key, value in initial_start_vars.iteritems(): -- cgit v1.2.3