aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--wqflask/wqflask/marker_regression/display_mapping_results.py3
-rw-r--r--wqflask/wqflask/marker_regression/qtlreaper_mapping.py116
-rw-r--r--wqflask/wqflask/marker_regression/run_mapping.py49
-rw-r--r--wqflask/wqflask/templates/show_trait_mapping_tools.html17
-rw-r--r--wqflask/wqflask/views.py3
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
@@ -116,6 +116,15 @@
<div class="tab-pane" id="interval_mapping">
<div style="margin-top: 20px" class="form-horizontal">
<div class="mapping_method_fields form-group">
+ <label for="reaper_version" style="text-align: right;" class="col-xs-3 control-label">Version<sup><a title="'New' is the new qtlreaper implementation written in RUST by Christian Fischer. 'Original' corresponds to the original version written in C.">?</a></sup></label>
+ <div style="margin-left:20px;" class="col-xs-3 controls">
+ <select id="reaper_version" class="form-control">
+ <option value="new">New</option>
+ <option value="original">Original</option>
+ </select>
+ </div>
+ </div>
+ <div class="mapping_method_fields form-group">
<label for="chr_select" style="text-align: right;" class="col-xs-3 control-label">Chromosome</label>
<div style="margin-left:20px;" class="col-xs-2 controls">
<select id="chr_reaper" class="form-control">
@@ -347,9 +356,11 @@
<dt>GEMMA</dt>
<dd>GEMMA is software implementing the Genome-wide Efficient Mixed Model Association algorithm for a standard linear mixed model for genome-wide association studies (GWAS).</dd>
{% if dataset.group.mapping_id == "1" %}
- <dt>Interval Mapping</dt>
- <dd>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.</dd>
- <dt>R/qtl</dt>
+ <dt style="margin-top: 20px;">Interval Mapping</dt>
+ <dd>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.<br><br>
+ 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.<br>
+ <i>Please let us know if you encounter any issues while using the new version.</i></dd>
+ <dt style="margin-top: 20px;">R/qtl</dt>
<dd>R/qtl is an extensible, interactive environment for mapping quantitative trait loci (QTL) in experimental crosses.</dd>
{% endif %}
</dl>
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():