Regular and single chr views work for qtl reaper now

Single chr now works for R/qtl (didn't before)

Merged qtl reaper code into marker_regression.py

7 files changed, 135 insertions, 225 deletions

diff --git a/wqflask/wqflask/interval_mapping/interval_mapping.py b/wqflask/wqflask/interval_mapping/interval_mapping.py
index 46aac86c..e8ac99d7 100755
--- a/wqflask/wqflask/interval_mapping/interval_mapping.py
+++ b/wqflask/wqflask/interval_mapping/interval_mapping.py
@@ -58,7 +58,7 @@ class IntervalMapping(object):
 
         self.json_data = {}
         self.json_data['lodnames'] = ['lod.hk']
-        self.gen_reaper_results(tempdata)
+        self.gen_reaper_results()
 
         #Get chromosome lengths for drawing the interval map plot
         chromosome_mb_lengths = {}
@@ -102,76 +102,8 @@ class IntervalMapping(object):
         else:
             self.control_locus = None
 
-    def gen_qtl_results(self, tempdata):
-        """Generates qtl results for plotting interval map"""
 
-        self.dataset.group.get_markers()
-        genotype = self.dataset.group.read_genotype_file()
-
-        samples, values, variances = self.this_trait.export_informative()
-
-        trimmed_samples = []
-        trimmed_values = []
-        for i in range(0, len(samples)):
-            if samples[i] in self.dataset.group.samplelist:
-                trimmed_samples.append(samples[i])
-                trimmed_values.append(values[i])
-
-        #if self.weighted_regression:
-        #    self.lrs_array = self.genotype.permutation(strains = trimmed_samples,
-        #                                                         trait = trimmed_values, 
-        #                                                         variance = variances,
-        #                                                         nperm=self.num_permutations)
-        #else:
-        self.lrs_array = genotype.permutation(strains = trimmed_samples,
-                                                   trait = trimmed_values, 
-                                                   nperm= self.num_permutations)
-        self.suggestive = self.lrs_array[int(self.num_permutations*0.37-1)]
-        self.significant = self.lrs_array[int(self.num_permutations*0.95-1)]
-
-        print("samples:", trimmed_samples)
-
-        if self.control_locus:
-            #if self.weighted_regression:
-            #    self.qtl_results = self.dataset.genotype.regression(strains = samples,
-            #                                                  trait = values,
-            #                                                  variance = variances,
-            #                                                  control = self.control_locus)
-            #else:
-            reaper_results = genotype.regression(strains = trimmed_samples,
-                                                          trait = trimmed_values,
-                                                          control = self.control_locus)
-        else:
-            #if self.weighted_regression:
-            #    self.qtl_results = self.dataset.genotype.regression(strains = samples,
-            #                                                  trait = values,
-            #                                                  variance = variances)
-            #else:
-            print("strains:", trimmed_samples)
-            print("trait:", trimmed_values)
-            reaper_results = genotype.regression(strains = trimmed_samples,
-                                                          trait = trimmed_values)
-
-        #Need to convert the QTL objects that qtl reaper returns into a json serializable dictionary
-        self.qtl_results = []
-        self.json_data['chr'] = []
-        self.json_data['pos'] = []
-        self.json_data['lod.hk'] = []
-        self.json_data['markernames'] = []
-        for qtl in reaper_results:
-            reaper_locus = qtl.locus
-            self.json_data['pos'].append(reaper_locus.Mb)
-            self.json_data['lod.hk'].append(qtl.lrs)
-            self.json_data['markernames'].append(reaper_locus.name)
-            if self.additive:
-                self.json_data['additive'] = []
-                self.json_data['additive'].append(qtl.additive)
-            locus = {"name":reaper_locus.name, "chr":reaper_locus.chr, "cM":reaper_locus.cM, "Mb":reaper_locus.Mb}
-            qtl = {"lrs": qtl.lrs, "locus": locus, "additive": qtl.additive}
-            self.qtl_results.append(qtl)
-
-
-    def gen_reaper_results(self, tempdata):
+    def gen_reaper_results(self):
         genotype = self.dataset.group.read_genotype_file()
 
         samples, values, variances = self.this_trait.export_informative()
@@ -249,127 +181,3 @@ class IntervalMapping(object):
             trimmed_genotype_data.append(new_genotypes)
         return trimmed_genotype_data
 
-    #def get_qtl_results(self):
-    #    """Gets the LOD (or LRS) score at each marker in order do the qtl mapping"""
-    #
-    #    
-    #    
-    #    #self.genotype = self.genotype.addinterval()
-    #    #resultSlice = []
-    #    #controlGeno = []
-    #    
-    #    #if self.multipleInterval:
-    #    #    self.suggestive = 0
-    #    #    self.significance = 0
-    #    #    if self.selectedChr > -1:
-    #    #        self.genotype.chromosome = [self.genotype[self.selectedChr]]
-    #    #else:
-    #    #single interval mapping
-    #    #try:
-    #    #    self.suggestive = float(fd.formdata.getvalue('permSuggestive'))
-    #    #    self.significance = float(fd.formdata.getvalue('permSignificance'))
-    #    #except:
-    #    #    self.suggestive = None
-    #    #    self.significance = None
-    #    
-    #    #NOT DOING MULTIPLE TRAITS YET, BUT WILL NEED TO LATER
-    #    #_strains, _vals, _vars = self.traitList[0].exportInformative(weightedRegression)
-    #    
-    #    #if webqtlUtil.ListNotNull(_vars):
-    #    #    pass
-    #    #else:
-    #    #    weightedRegression = 0
-    #    #    _strains, _vals, _vars = self.traitList[0].exportInformative()
-    #        
-    #    ##locate genotype of control Locus
-    #    #if self.controlLocus:
-    #    #    controlGeno2 = []
-    #    #    _FIND = 0
-    #    #    for _chr in self.genotype:
-    #    #        for _locus in _chr:
-    #    #            if _locus.name == self.controlLocus:
-    #    #                controlGeno2 = _locus.genotype
-    #    #                _FIND = 1
-    #    #                break
-    #    #        if _FIND:
-    #    #            break
-    #    #    if controlGeno2:
-    #    #        _prgy = list(self.genotype.prgy)
-    #    #        for _strain in _strains:
-    #    #            _idx = _prgy.index(_strain)
-    #    #            controlGeno.append(controlGeno2[_idx])
-    #    #    else:
-    #    #        return "The control marker you selected is not in the genofile."
-    #    #
-    #    #    
-    #    #    if self.significance and self.suggestive:
-    #    #        pass
-    #    #    else:
-    #    #        if self.permChecked:
-    #    #            if weightedRegression:
-    #    #                self.LRSArray = self.genotype.permutation(strains = _strains, trait = _vals, 
-    #    #                    variance = _vars, nperm=fd.nperm)
-    #    #            else:
-    #    #                self.LRSArray = self.genotype.permutation(strains = _strains, trait = _vals, 
-    #    #                    nperm=fd.nperm)
-    #    #            self.suggestive = self.LRSArray[int(fd.nperm*0.37-1)]
-    #    #            self.significance = self.LRSArray[int(fd.nperm*0.95-1)]
-    #    #
-    #    #        else:
-    #    #            self.suggestive = 9.2
-    #    #            self.significance = 16.1
-    #    #
-    #    #    #calculating bootstrap
-    #    #    #from now on, genotype could only contain a single chromosome 
-    #    #    #permutation need to be performed genome wide, this is not the case for bootstrap
-    #    #    
-    #    #    #due to the design of qtlreaper, composite regression need to be performed genome wide
-    #    #    if not self.controlLocus and self.selectedChr > -1:
-    #    #        self.genotype.chromosome = [self.genotype[self.selectedChr]]
-    #    #    elif self.selectedChr > -1: #self.controlLocus and self.selectedChr > -1
-    #    #        lociPerChr = map(len, self.genotype)
-    #    #        resultSlice = reduce(lambda X, Y: X+Y, lociPerChr[:self.selectedChr], 0)
-    #    #        resultSlice = [resultSlice,resultSlice+lociPerChr[self.selectedChr]]
-    #    #    else:
-    #    #        pass
-    #        
-    #    #calculate QTL for each trait
-    #    self.qtl_results = []
-    #
-    #    #for thisTrait in self.traitList:
-    #    _strains, _vals, _vars = thisTrait.exportInformative(weightedRegression)
-    #    if self.controlLocus:
-    #        if weightedRegression:
-    #            qtlresult = self.genotype.regression(strains = _strains, trait = _vals, 
-    #                    variance = _vars, control = self.controlLocus)
-    #        else:
-    #            qtlresult = self.genotype.regression(strains = _strains, trait = _vals,
-    #                control = self.controlLocus)
-    #        if resultSlice:
-    #            qtlresult = qtlresult[resultSlice[0]:resultSlice[1]]
-    #    else:
-    #        if weightedRegression:
-    #            qtlresult = self.genotype.regression(strains = _strains, trait = _vals, 
-    #                    variance = _vars)
-    #        else:
-    #            qtlresult = self.genotype.regression(strains = _strains, trait = _vals)
-    #
-    #    self.qtlresults.append(qtlresult)
-    #    
-    #    if not self.multipleInterval:
-    #        if self.controlLocus and self.selectedChr > -1:
-    #            self.genotype.chromosome = [self.genotype[self.selectedChr]]
-    #            
-    #        if self.bootChecked:
-    #            if controlGeno:
-    #                self.bootResult = self.genotype.bootstrap(strains = _strains, trait = _vals,
-    #                    control = controlGeno, nboot=fd.nboot)
-    #            elif weightedRegression:
-    #                self.bootResult = self.genotype.bootstrap(strains = _strains, trait = _vals, 
-    #                    variance = _vars, nboot=fd.nboot)
-    #            else:
-    #                self.bootResult = self.genotype.bootstrap(strains = _strains, trait = _vals, 
-    #                    nboot=fd.nboot)
-    #        else:
-    #            self.bootResult = []
-
diff --git a/wqflask/wqflask/marker_regression/marker_regression.py b/wqflask/wqflask/marker_regression/marker_regression.py
index c396b14b..605fde63 100755
--- a/wqflask/wqflask/marker_regression/marker_regression.py
+++ b/wqflask/wqflask/marker_regression/marker_regression.py
@@ -110,6 +110,15 @@ class MarkerRegression(object):
 
             results = self.run_rqtl_geno()
             #print("qtl_results:", results)
+        elif self.mapping_method == "reaper":
+            if start_vars['num_perm'] == "":
+                self.num_perm = 0
+            else:
+                self.num_perm = int(start_vars['num_perm'])
+            self.additive = False
+            self.control = start_vars['control_marker']
+            self.do_control = start_vars['do_control']
+            results = self.gen_reaper_results()
         elif self.mapping_method == "plink":
             results = self.run_plink()
             #print("qtl_results:", pf(results))
@@ -130,7 +139,7 @@ class MarkerRegression(object):
                 if marker['chr1'] > 0 or marker['chr1'] == "X" or marker['chr1'] == "X/Y":
                     if marker['chr1'] > highest_chr or marker['chr1'] == "X" or marker['chr1'] == "X/Y":
                         highest_chr = marker['chr1']
-                    if 'lod_score' in marker:
+                    if 'lod_score' in marker.keys():
                         self.qtl_results.append(marker)
 
             for qtl in enumerate(self.qtl_results):
@@ -157,7 +166,7 @@ class MarkerRegression(object):
                 if marker['chr'] > 0 or marker['chr'] == "X" or marker['chr'] == "X/Y":
                     if marker['chr'] > highest_chr or marker['chr'] == "X" or marker['chr'] == "X/Y":
                         highest_chr = marker['chr']
-                    if 'lod_score' in marker:
+                    if ('lod_score' in marker.keys()) or ('lrs_value' in marker.keys()):
                         self.qtl_results.append(marker)
 
             self.json_data['chr'] = []
@@ -178,7 +187,7 @@ class MarkerRegression(object):
                 else:
                     self.json_data['chr'].append(str(qtl['chr']))
                 self.json_data['pos'].append(qtl['Mb'])
-                if 'lrs_value' in qtl:
+                if 'lrs_value' in qtl.keys():
                     self.json_data['lod.hk'].append(str(qtl['lrs_value']))
                 else:
                     self.json_data['lod.hk'].append(str(qtl['lod_score']))
@@ -588,6 +597,66 @@ class MarkerRegression(object):
         
         return sample_list
     
+    def gen_reaper_results(self):
+        genotype = self.dataset.group.read_genotype_file()
+
+        samples, values, variances = self.this_trait.export_informative()
+
+        trimmed_samples = []
+        trimmed_values = []
+        for i in range(0, len(samples)):
+            if samples[i] in self.dataset.group.samplelist:
+                trimmed_samples.append(samples[i])
+                trimmed_values.append(values[i])
+                
+        self.lrs_array = genotype.permutation(strains = trimmed_samples,
+                                                   trait = trimmed_values, 
+                                                   nperm= self.num_perm)
+        
+        self.suggestive = self.lrs_array[int(self.num_perm*0.37-1)]
+        self.significant = self.lrs_array[int(self.num_perm*0.95-1)]
+        self.json_data['suggestive'] = self.suggestive
+        self.json_data['significant'] = self.significant
+
+        print("samples:", trimmed_samples)
+
+        if self.control != "" and self.do_control == "true":
+            reaper_results = genotype.regression(strains = trimmed_samples,
+                                                          trait = trimmed_values,
+                                                          control = self.control)
+        else:
+            reaper_results = genotype.regression(strains = trimmed_samples,
+                                                          trait = trimmed_values)
+
+        self.json_data['chr'] = []
+        self.json_data['pos'] = []
+        self.json_data['lod.hk'] = []
+        self.json_data['markernames'] = []
+        if self.additive:
+            self.json_data['additive'] = []
+
+        #Need to convert the QTL objects that qtl reaper returns into a json serializable dictionary
+        qtl_results = []
+        for qtl in reaper_results:
+            reaper_locus = qtl.locus
+            #ZS: Convert chr to int
+            converted_chr = reaper_locus.chr
+            if reaper_locus.chr != "X" and reaper_locus.chr != "X/Y":
+                converted_chr = int(reaper_locus.chr)
+            self.json_data['chr'].append(converted_chr)
+            self.json_data['pos'].append(reaper_locus.Mb)
+            self.json_data['lod.hk'].append(qtl.lrs)
+            self.json_data['markernames'].append(reaper_locus.name)
+            if self.additive:
+                self.json_data['additive'].append(qtl.additive)
+            locus = {"name":reaper_locus.name, "chr":reaper_locus.chr, "cM":reaper_locus.cM, "Mb":reaper_locus.Mb}
+            qtl = {"lrs_value": qtl.lrs, "chr":converted_chr, "Mb":reaper_locus.Mb,
+                   "cM":reaper_locus.cM, "name":reaper_locus.name, "additive":qtl.additive, "dominance":qtl.dominance}
+            qtl_results.append(qtl)
+
+        return qtl_results
+
+
     def parse_plink_output(self, output_filename):
         plink_results={}
     
diff --git a/wqflask/wqflask/marker_regression/marker_regression_gn1.py b/wqflask/wqflask/marker_regression/marker_regression_gn1.py
index 23d5a631..a7b5fad3 100644
--- a/wqflask/wqflask/marker_regression/marker_regression_gn1.py
+++ b/wqflask/wqflask/marker_regression/marker_regression_gn1.py
@@ -178,9 +178,13 @@ class MarkerRegression(object):
         self.this_trait = start_vars['this_trait']
         self.species = start_vars['species']
 
-        self.vals = start_vars['vals'] #Needed to put into form for when it is resubmitted for single chr views or remapping
- 
-        self.mapping_method = start_vars['mapping_method'] #Needing for form submission when doing single chr mapping or remapping after changing options
+        #Needing for form submission when doing single chr mapping or remapping after changing options
+        self.vals = start_vars['vals'] 
+        self.mapping_method = start_vars['mapping_method'] 
+        if self.mapping_method == "rqtl_geno":
+            self.mapmethod_rqtl_geno = start_vars['method']
+            self.mapmodel_rqtl_geno = start_vars['model']
+            self.pair_scan = start_vars['pair_scan']
 
         self.js_data = start_vars['js_data']
 
@@ -225,10 +229,15 @@ class MarkerRegression(object):
         #self.permChecked = fd.formdata.getvalue('permCheck', True)
         self.bootChecked = False #ZS: For now setting to False, I'll add this option later once rest of figure works
         #self.bootChecked = fd.formdata.getvalue('bootCheck', '')
+        if 'do_control' in start_vars.keys():
+            self.doControl = start_vars['do_control']
+        else:
+            self.doControl = "false"
         if 'control' in start_vars.keys():
             self.controlLocus = start_vars['control']
         else:
             self.controlLocus = ""
+        
         #self.controlLocus = fd.formdata.getvalue('controlLocus', '')
 
         #try:
@@ -1760,6 +1769,7 @@ class MarkerRegression(object):
             lineColor = pid.lightblue
             startPosX = xLeftOffset
             for j, ChrInfo in enumerate(ChrAInfo):
+              if ChrInfo == self.selectedChr:
                 preLpos = -1
                 for i, item in enumerate(ChrInfo):
                     Lname,Lpos = item
@@ -1932,8 +1942,6 @@ class MarkerRegression(object):
                 previous_chr = qtlresult['chr']
                 previous_chr_as_int += 1
 
-                print("ChrLengthDistList:", self.ChrLengthDistList)
-
                 newStartPosX = (self.ChrLengthDistList[previous_chr_as_int - 1]+self.GraphInterval)*plotXScale
                 if newStartPosX != oldStartPosX:
                     startPosX += newStartPosX
diff --git a/wqflask/wqflask/static/new/javascript/show_trait_mapping_tools.js b/wqflask/wqflask/static/new/javascript/show_trait_mapping_tools.js
index 16eeb8df..a48ac2b7 100755
--- a/wqflask/wqflask/static/new/javascript/show_trait_mapping_tools.js
+++ b/wqflask/wqflask/static/new/javascript/show_trait_mapping_tools.js
@@ -236,9 +236,11 @@
       var form_data, url;
       console.log("In interval mapping");
       //$("#progress_bar_container").modal();
-      url = "/interval_mapping";
+      url = "/marker_regression";
       $('input[name=method]').val("reaper");
       $('input[name=manhattan_plot]').val($('input[name=manhattan_plot_reaper]:checked').val());
+      $('input[name=control_marker]').val($('input[name=control_reaper]').val());
+      $('input[name=do_control]').val($('input[name=do_control_reaper]:checked').val());
       $('input[name=mapping_display_all]').val($('input[name=display_all_reaper]'));
       $('input[name=suggestive]').val($('input[name=suggestive_reaper]'));
       form_data = $('#trait_data_form').serialize();
diff --git a/wqflask/wqflask/templates/marker_regression_gn1.html b/wqflask/wqflask/templates/marker_regression_gn1.html
index 05623462..9f7dd41f 100644
--- a/wqflask/wqflask/templates/marker_regression_gn1.html
+++ b/wqflask/wqflask/templates/marker_regression_gn1.html
@@ -22,6 +22,11 @@
         <input type="hidden" name="selected_chr" value="{{ selectedChr }}">
         <input type="hidden" name="manhattan_plot">
         <input type="hidden" name="num_perm">
+        <input type="hidden" name="do_control" value="{{ doControl }}">
+        <input type="hidden" name="control_marker" value="{{ controlLocus }}">
+        <input type="hidden" name="mapmethod_rqtl_geno" value="{{ mapmethod_rqtl_geno }}">
+        <input type="hidden" name="mapmodel_rqtl_geno" value="{{ mapmodel_rqtl_geno }}">
+        <input type="hidden" name="pair_scan" value="{{ pair_scan }}">
         <div id="chart_container">
             <div class="qtlcharts" id="topchart">
                 {{ gifmap|safe }}
diff --git a/wqflask/wqflask/templates/show_trait_mapping_tools.html b/wqflask/wqflask/templates/show_trait_mapping_tools.html
index de25f4ee..5b0c74c3 100755
--- a/wqflask/wqflask/templates/show_trait_mapping_tools.html
+++ b/wqflask/wqflask/templates/show_trait_mapping_tools.html
@@ -187,6 +187,24 @@
                                 <input name="num_perm_reaper" value="2000" type="text" class="form-control">
                             </div>
                         </div>
+                        <div class="mapping_method_fields form-group">
+                            <label for="control_for" class="col-xs-3 control-label">Control&nbsp;for</label>
+                            <div style="margin-left: 20px;" class="col-xs-4 controls">
+                                {% if dataset.type == 'ProbeSet' and this_trait.locus_chr != "" %}
+                                <input name="control_reaper" value="{{ nearest_marker }}" type="text" style="width: 160px;" class="form-control" />
+                                {% else %}
+                                <input name="control_reaper" value="" type="text" class="form-control" />
+                                {% endif %}
+                                <label class="radio-inline">
+                                    <input type="radio" name="do_control_reaper" value="true">
+                                    Yes
+                                </label>
+                                <label class="radio-inline">
+                                    <input type="radio" name="do_control_reaper" value="false" checked="">
+                                    No
+                                </label>
+                            </div>
+                        </div>
 <!-- 
                         <div class="mapping_method_fields form-group">
                             <label for="mapping_bootstraps" class="col-xs-3 control-label" title="Bootstrapping Resamples">Bootstrap</label>
diff --git a/wqflask/wqflask/views.py b/wqflask/wqflask/views.py
index 9c5a3c92..73d9ddab 100755
--- a/wqflask/wqflask/views.py
+++ b/wqflask/wqflask/views.py
@@ -411,28 +411,6 @@ def marker_regression_page():
 
     return rendered_template
 
-@app.route("/export", methods = ('POST',))
-def export():
-    print("request.form:", request.form)
-    svg_xml = request.form.get("data", "Invalid data")
-    filename = request.form.get("filename", "manhattan_plot_snp")
-    response = Response(svg_xml, mimetype="image/svg+xml")
-    response.headers["Content-Disposition"] = "attachment; filename=%s"%filename
-    return response
-
-@app.route("/export_pdf", methods = ('POST',))
-def export_pdf():
-    import cairosvg
-    print("request.form:", request.form)
-    svg_xml = request.form.get("data", "Invalid data")
-    print("svg_xml:", svg_xml)
-    filename = request.form.get("filename", "interval_map_pdf")
-    filepath = "/home/zas1024/gene/wqflask/output/"+filename
-    pdf_file = cairosvg.svg2pdf(bytestring=svg_xml)
-    response = Response(pdf_file, mimetype="application/pdf")
-    response.headers["Content-Disposition"] = "attachment; filename=%s"%filename
-    return response
-
 @app.route("/interval_mapping", methods=('POST',))
 def interval_mapping_page():
     initial_start_vars = request.form
@@ -487,6 +465,28 @@ def interval_mapping_page():
 
     return rendered_template
 
+@app.route("/export", methods = ('POST',))
+def export():
+    print("request.form:", request.form)
+    svg_xml = request.form.get("data", "Invalid data")
+    filename = request.form.get("filename", "manhattan_plot_snp")
+    response = Response(svg_xml, mimetype="image/svg+xml")
+    response.headers["Content-Disposition"] = "attachment; filename=%s"%filename
+    return response
+
+@app.route("/export_pdf", methods = ('POST',))
+def export_pdf():
+    import cairosvg
+    print("request.form:", request.form)
+    svg_xml = request.form.get("data", "Invalid data")
+    print("svg_xml:", svg_xml)
+    filename = request.form.get("filename", "interval_map_pdf")
+    filepath = "/home/zas1024/gene/wqflask/output/"+filename
+    pdf_file = cairosvg.svg2pdf(bytestring=svg_xml)
+    response = Response(pdf_file, mimetype="application/pdf")
+    response.headers["Content-Disposition"] = "attachment; filename=%s"%filename
+    return response
+
 @app.route("/corr_compute", methods=('POST',))
 def corr_compute_page():
     print("In corr_compute, request.form is:", pf(request.form))