about summary refs log tree commit diff
diff options
context:
space:
mode:
-rw-r--r--wqflask/utility/gen_geno_ob.py4
-rw-r--r--wqflask/wqflask/marker_regression/display_mapping_results.py36
-rw-r--r--wqflask/wqflask/marker_regression/rqtl_mapping.py116
-rw-r--r--wqflask/wqflask/marker_regression/run_mapping.py4
-rw-r--r--wqflask/wqflask/static/new/javascript/show_trait.js7
-rw-r--r--wqflask/wqflask/templates/mapping_results.html2
-rw-r--r--wqflask/wqflask/templates/show_trait_mapping_tools.html18
7 files changed, 141 insertions, 46 deletions
diff --git a/wqflask/utility/gen_geno_ob.py b/wqflask/utility/gen_geno_ob.py
index aa5b27c4..db40f6ea 100644
--- a/wqflask/utility/gen_geno_ob.py
+++ b/wqflask/utility/gen_geno_ob.py
@@ -45,10 +45,12 @@ class genotype(object):
         chr_ob = None
         for marker in qtl_results:
             locus = Locus(self)
-            if str(marker['chr']) != this_chr:
+            if (str(marker['chr']) != this_chr) and this_chr != "X": #ZS: This is really awkward but works as a temporary fix
                 if this_chr != "":
                     self.chromosomes.append(chr_ob)
                 this_chr = str(marker['chr'])
+                if this_chr == "20":
+                    this_chr = "X"
                 chr_ob = Chr(this_chr, self)
             if 'chr' in marker:
                 locus.chr = str(marker['chr'])
diff --git a/wqflask/wqflask/marker_regression/display_mapping_results.py b/wqflask/wqflask/marker_regression/display_mapping_results.py
index 302a4ff9..b8f84721 100644
--- a/wqflask/wqflask/marker_regression/display_mapping_results.py
+++ b/wqflask/wqflask/marker_regression/display_mapping_results.py
@@ -334,6 +334,10 @@ class DisplayMappingResults(object):
         ################################################################
         self.ChrList = [("All", -1)]
         for i, indChr in enumerate(self.genotype):
+            if self.dataset.group.species == "mouse" and indChr.name == "20":
+                self.ChrList.append(("X", i))
+            elif self.dataset.group.species == "rat" and indChr.name == "21":
+                self.ChrList.append(("X", i))
             self.ChrList.append((indChr.name, i))
 
         self.ChrLengthMbList = g.db.execute("""
@@ -356,10 +360,7 @@ class DisplayMappingResults(object):
 
         self.ChrLengthCMList = []
         for i, _chr in enumerate(self.genotype):
-            if self.mapping_method == "rqtl_geno" and self.genotype.filler == True:
-                self.ChrLengthCMList.append(_chr[-1].cM)
-            else:
-                self.ChrLengthCMList.append(_chr[-1].cM - _chr[0].cM)
+            self.ChrLengthCMList.append(_chr[-1].cM - _chr[0].cM)
 
         self.ChrLengthCMSum = reduce(lambda x, y:x+y, self.ChrLengthCMList, 0.0)
 
@@ -1568,10 +1569,7 @@ class DisplayMappingResults(object):
             if self.selectedChr == -1: #ZS: If viewing full genome/all chromosomes
                 for i, _chr in enumerate(self.genotype):
                     thisChr = []
-                    if self.mapping_method == "rqtl_geno" and self.genotype.filler == True:
-                        Locus0CM = 0
-                    else:
-                        Locus0CM = _chr[0].cM
+                    Locus0CM = _chr[0].cM
                     nLoci = len(_chr)
                     if  nLoci <= 8:
                         for _locus in _chr:
@@ -1593,10 +1591,7 @@ class DisplayMappingResults(object):
                 for i, _chr in enumerate(self.genotype):
                     if _chr.name == self.ChrList[self.selectedChr][0]:
                         thisChr = []
-                        if self.mapping_method == "rqtl_geno" and self.genotype.filler == True:
-                            Locus0CM = 0
-                        else:
-                            Locus0CM = _chr[0].cM
+                        Locus0CM = _chr[0].cM
                         for _locus in _chr:
                             if _locus.name != ' - ':
                                 if _locus.cM != preLpos:
@@ -1892,12 +1887,21 @@ class DisplayMappingResults(object):
                     oldStartPosX = newStartPosX
 
             #ZS: This is beause the chromosome value stored in qtlresult['chr'] can be (for example) either X or 20 depending upon the mapping method/scale used
-            if self.plotScale == "physic":
-                this_chr = str(self.ChrList[self.selectedChr][0])
-            else:
+            this_chr = str(self.ChrList[self.selectedChr][0])
+            if self.plotScale != "physic":
                 this_chr = str(self.ChrList[self.selectedChr][1]+1)
+
             if self.selectedChr == -1 or str(qtlresult['chr']) == this_chr:
-                Xc = startPosX + (qtlresult['Mb']-startMb)*plotXScale
+                if self.plotScale != "physic" and self.genotype.filler == True:
+                    if self.selectedChr != -1:
+                        start_cm = self.genotype[self.selectedChr - 1][0].cM
+                        Xc = startPosX + (qtlresult['Mb'] - start_cm)*plotXScale
+                    else:
+                        start_cm = self.genotype[previous_chr_as_int][0].cM
+                        Xc = startPosX + ((qtlresult['Mb']-start_cm-startMb)*plotXScale)*(((qtlresult['Mb']-start_cm-startMb)*plotXScale)/((qtlresult['Mb']-start_cm-startMb+self.GraphInterval)*plotXScale))
+                else:
+                    Xc = startPosX + (qtlresult['Mb']-startMb)*plotXScale
+
                 # updated by NL 06-18-2011:
                 # fix the over limit LRS graph issue since genotype trait may give infinite LRS;
                 # for any lrs is over than 460(LRS max in this system), it will be reset to 460
diff --git a/wqflask/wqflask/marker_regression/rqtl_mapping.py b/wqflask/wqflask/marker_regression/rqtl_mapping.py
index 2e3ea406..d76f3812 100644
--- a/wqflask/wqflask/marker_regression/rqtl_mapping.py
+++ b/wqflask/wqflask/marker_regression/rqtl_mapping.py
@@ -1,14 +1,17 @@
 import rpy2.robjects as ro
+import rpy2.robjects.numpy2ri as np2r
 import numpy as np
 
 from base.webqtlConfig import TMPDIR
+from base.trait import GeneralTrait
+from base.data_set import create_dataset
 from utility import webqtlUtil
 from utility.tools import locate, TEMPDIR
 
 import utility.logger
 logger = utility.logger.getLogger(__name__ )
 
-def run_rqtl_geno(vals, dataset, method, model, permCheck, num_perm, do_control, control_marker, manhattan_plot, pair_scan):
+def run_rqtl_geno(vals, dataset, method, model, permCheck, num_perm, do_control, control_marker, manhattan_plot, pair_scan, samples, cofactors):
     ## Get pointers to some common R functions
     r_library     = ro.r["library"]                 # Map the library function
     r_c           = ro.r["c"]                       # Map the c function
@@ -43,14 +46,22 @@ def run_rqtl_geno(vals, dataset, method, model, permCheck, num_perm, do_control,
     else:
         cross_object = calc_genoprob(cross_object, step=1, stepwidth="max")
 
-    cross_object = add_phenotype(cross_object, sanitize_rqtl_phenotype(vals))                 # Add the phenotype
+    cross_object = add_phenotype(cross_object, sanitize_rqtl_phenotype(vals), "the_pheno")                 # Add the phenotype
 
     # Scan for QTLs
-    covar = create_covariates(control_marker, cross_object)                                   # Create the additive covariate matrix
+    marker_covars = create_marker_covariates(control_marker, cross_object)  # Create the additive covariate markers
+
+    if cofactors != "":
+        cross_object, trait_covars = add_cofactors(cross_object, dataset, cofactors, samples)                            # Create the covariates from selected traits
+        ro.r('all_covars <- cbind(marker_covars, trait_covars)')
+    else:
+        ro.r('all_covars <- marker_covars')
+
+    covars = ro.r['all_covars']
 
     if pair_scan:
         if do_control == "true":
-            logger.info("Using covariate"); result_data_frame = scantwo(cross_object, pheno = "the_pheno", addcovar = covar, model=model, method=method, n_cluster = 16)
+            logger.info("Using covariate"); result_data_frame = scantwo(cross_object, pheno = "the_pheno", addcovar = covars, model=model, method=method, n_cluster = 16)
         else:
             logger.info("No covariates"); result_data_frame = scantwo(cross_object, pheno = "the_pheno", model=model, method=method, n_cluster = 16)
 
@@ -61,14 +72,14 @@ def run_rqtl_geno(vals, dataset, method, model, permCheck, num_perm, do_control,
 
         return process_pair_scan_results(result_data_frame)
     else:
-        if do_control == "true":
-            logger.info("Using covariate"); result_data_frame = scanone(cross_object, pheno = "the_pheno", addcovar = covar, model=model, method=method)
+        if do_control == "true" or cofactors != "":
+            logger.info("Using covariate"); result_data_frame = scanone(cross_object, pheno = "the_pheno", addcovar = covars, model=model, method=method)
         else:
             logger.info("No covariates"); result_data_frame = scanone(cross_object, pheno = "the_pheno", model=model, method=method)
 
         if num_perm > 0 and permCheck == "ON":                                                                   # Do permutation (if requested by user)
-            if do_control == "true":
-                perm_data_frame = scanone(cross_object, pheno_col = "the_pheno", addcovar = covar, n_perm = num_perm, model=model, method=method)
+            if do_control == "true" or cofactors != "":
+                perm_data_frame = scanone(cross_object, pheno_col = "the_pheno", addcovar = covars, n_perm = num_perm, model=model, method=method)
             else:
                 perm_data_frame = scanone(cross_object, pheno_col = "the_pheno", n_perm = num_perm, model=model, method=method)
 
@@ -118,23 +129,6 @@ def generate_cross_from_geno(dataset):        # TODO: Need to figure out why som
       }
     """ % (dataset.group.name + ".geno"))
 
-def add_phenotype(cross, pheno_as_string):
-    ro.globalenv["the_cross"] = cross
-    ro.r('the_cross$pheno <- cbind(pull.pheno(the_cross), the_pheno = '+ pheno_as_string +')')
-    return ro.r["the_cross"]
-
-def create_covariates(control_marker, cross):
-    ro.globalenv["the_cross"] = cross
-    ro.r('genotypes <- pull.geno(the_cross)')                             # Get the genotype matrix
-    userinputS = control_marker.replace(" ", "").split(",")               # TODO: sanitize user input, Never Ever trust a user
-    covariate_names = ', '.join('"{0}"'.format(w) for w in userinputS)
-    ro.r('covnames <- c(' + covariate_names + ')')
-    ro.r('covInGeno <- which(covnames %in% colnames(genotypes))')
-    ro.r('covnames <- covnames[covInGeno]')
-    ro.r("cat('covnames (purged): ', covnames,'\n')")
-    ro.r('covariates <- genotypes[,covnames]')                            # Get the covariate matrix by using the marker name as index to the genotype file
-    return ro.r["covariates"]
-
 def sanitize_rqtl_phenotype(vals):
     pheno_as_string = "c("
     for i, val in enumerate(vals):
@@ -151,6 +145,78 @@ def sanitize_rqtl_phenotype(vals):
     pheno_as_string += ")"
     return pheno_as_string
 
+def add_phenotype(cross, pheno_as_string, col_name):
+    ro.globalenv["the_cross"] = cross
+    ro.r('the_cross$pheno <- cbind(pull.pheno(the_cross), ' + col_name + ' = '+ pheno_as_string +')')
+    return ro.r["the_cross"]
+
+def pull_covar(cross, covar_name_string):
+    ro.globalenv["the_cross"] = cross
+    ro.r('trait_covars <- pull.pheno(the_cross, ' + covar_name_string + ')')
+
+    return ro.r["trait_covars"]
+
+def add_cofactors(cross, this_dataset, covariates, samples):
+    ro.numpy2ri.activate()
+
+    covariate_list = covariates.split(",")
+    covar_name_string = "c("
+    for i, covariate in enumerate(covariate_list):
+        this_covar_data = []
+        covar_as_string = "c("
+        trait_name = covariate.split(":")[0]
+        dataset_ob = create_dataset(covariate.split(":")[1])
+        trait_ob = GeneralTrait(dataset=dataset_ob,
+                                name=trait_name,
+                                cellid=None)
+
+        this_dataset.group.get_samplelist()
+        trait_samples = this_dataset.group.samplelist
+        trait_sample_data = trait_ob.data
+        for index, sample in enumerate(trait_samples):
+            if sample in samples:
+                if sample in trait_sample_data:
+                    sample_value = trait_sample_data[sample].value
+                    this_covar_data.append(sample_value)
+                else:
+                    this_covar_data.append("NA")
+
+        for j, item in enumerate(this_covar_data):
+            if j < (len(this_covar_data) - 1):
+                covar_as_string += str(item) + ","
+            else:
+                covar_as_string += str(item)
+
+        covar_as_string += ")"
+
+        col_name = "covar_" + str(i)
+
+        if i < (len(covariate_list) - 1):
+            covar_name_string += '"' + col_name + '", '
+        else:
+            covar_name_string += '"' + col_name + '"'
+
+        cross = add_phenotype(cross, covar_as_string, col_name)
+
+    covar_name_string += ")"
+
+    covars = pull_covar(cross, covar_name_string)
+
+    return cross, covars
+
+def create_marker_covariates(control_marker, cross):
+    ro.globalenv["the_cross"] = cross
+    ro.r('genotypes <- pull.geno(the_cross)')                             # Get the genotype matrix
+    userinputS = control_marker.replace(" ", "").split(",")               # TODO: sanitize user input, Never Ever trust a user
+    covariate_names = ', '.join('"{0}"'.format(w) for w in userinputS)
+    ro.r('covnames <- c(' + covariate_names + ')')
+    ro.r('covInGeno <- which(covnames %in% colnames(genotypes))')
+    ro.r('covnames <- covnames[covInGeno]')
+    ro.r("cat('covnames (purged): ', covnames,'\n')")
+    ro.r('marker_covars <- genotypes[,covnames]')                            # Get the covariate matrix by using the marker name as index to the genotype file
+
+    return ro.r["marker_covars"]
+
 def process_pair_scan_results(result):
     pair_scan_results = []
 
diff --git a/wqflask/wqflask/marker_regression/run_mapping.py b/wqflask/wqflask/marker_regression/run_mapping.py
index f03b046e..3006c4ff 100644
--- a/wqflask/wqflask/marker_regression/run_mapping.py
+++ b/wqflask/wqflask/marker_regression/run_mapping.py
@@ -216,9 +216,9 @@ class RunMapping(object):
             #if start_vars['pair_scan'] == "true":
             #    self.pair_scan = True
             if self.permCheck and self.num_perm > 0:
-                self.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)
+                self.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, self.samples, self.covariates)
             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)
+                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, self.samples, self.covariates)
         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/static/new/javascript/show_trait.js b/wqflask/wqflask/static/new/javascript/show_trait.js
index 0162f858..a2fa37e0 100644
--- a/wqflask/wqflask/static/new/javascript/show_trait.js
+++ b/wqflask/wqflask/static/new/javascript/show_trait.js
@@ -107,6 +107,13 @@ $("#remove_covariates").click(function () {
     $("input[name=covariates]").val("")
     $(".selected_covariates").val("")
 });
+$(".select_covariates").click(function () {
+  open_covariate_selection();
+});
+$(".remove_covariates").click(function () {
+  $("input[name=covariates]").val("")
+  $(".selected_covariates").val("")
+});
 d3.select("#clear_compare_trait").on("click", (function(_this) {
   return function() {
     return $('.scatter-matrix-container').remove();
diff --git a/wqflask/wqflask/templates/mapping_results.html b/wqflask/wqflask/templates/mapping_results.html
index 80928541..86816d17 100644
--- a/wqflask/wqflask/templates/mapping_results.html
+++ b/wqflask/wqflask/templates/mapping_results.html
@@ -478,7 +478,7 @@
           title: "<b>Histogram of Permutation Test</b>",
           xaxis: {
                   autorange: true,
-                  title: "<b>LRS</b>",
+                  title: "<b>{{ LRS_LOD }}</b>",
                   titlefont: {
                     family: "arial",
                     size: 20
diff --git a/wqflask/wqflask/templates/show_trait_mapping_tools.html b/wqflask/wqflask/templates/show_trait_mapping_tools.html
index 6da25a5d..66903c31 100644
--- a/wqflask/wqflask/templates/show_trait_mapping_tools.html
+++ b/wqflask/wqflask/templates/show_trait_mapping_tools.html
@@ -325,7 +325,23 @@
                                </label>
                             </div>
                         </div>
-
+                        <div class="mapping_method_fields form-group">
+                            <div class="col-xs-3" style="text-align: right;">
+                            <label class="control-label">Covariates</label>
+                            <font size="2">Select covariate(s) from a collection</font>
+                            </div>
+                            <div style="margin-left:20px;" class="col-xs-7">
+                              {% if g.user_session.num_collections < 1 %}
+                              No collections available. Please add traits to a collection to use them as covariates.
+                              {% else %}
+                              <div style="margin-bottom: 10px;">
+                                <button type="button" class="btn btn-default select_covariates" style="width: 80px; padding-right: 10px;">Select</button>
+                                <button type="button" class="btn btn-default remove_covariates" style="width: 80px;">Remove</button>
+                              </div>
+                              <textarea rows="3" cols="50" readonly placeholder="No covariates selected" style="overflow-y: scroll; resize: none; width: 200px;" class="selected_covariates"></textarea>
+                              {% endif %}
+                            </div>
+                        </div>
                         <div class="mapping_method_fields form-group">
                             <label class="col-xs-3 control-label"></label>
                             <div style="margin-left:20px;" class="col-xs-6">