Added Christian's RUST qtlreaper as an option

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():