Merge branch 'master' of https://github.com/zsloan/genenetwork2 into zsloan

4 files changed, 123 insertions, 89 deletions

diff --git a/wqflask/wqflask/marker_regression/marker_regression.py b/wqflask/wqflask/marker_regression/marker_regression.py
index 64d3ef3d..db92e14c 100755
--- a/wqflask/wqflask/marker_regression/marker_regression.py
+++ b/wqflask/wqflask/marker_regression/marker_regression.py
@@ -65,6 +65,7 @@ class MarkerRegression(object):
             self.manhattan_plot = True
         else:
             self.manhattan_plot = False
+
         self.maf = start_vars['maf'] # Minor allele frequency
         self.suggestive = ""
         self.significant = ""
@@ -76,9 +77,18 @@ class MarkerRegression(object):
         elif self.mapping_method == "rqtl_plink":
             qtl_results = self.run_rqtl_plink()
         elif self.mapping_method == "rqtl_geno":
-            self.num_perm = start_vars['num_perm']
+            if start_vars['num_perm'] == "":
+                self.num_perm = 0
+            else:
+                self.num_perm = start_vars['num_perm']
             self.control = start_vars['control_marker']
 
+            if start_vars['pair_scan'] == "true":
+                self.pair_scan = True
+            else:
+                self.pair_scan = False
+            print("pair scan:", self.pair_scan)
+
             print("DOING RQTL GENO")
             qtl_results = self.run_rqtl_geno()
             print("qtl_results:", qtl_results)
@@ -228,25 +238,25 @@ class MarkerRegression(object):
         os.system(rqtl_command)
         
         count, p_values = self.parse_rqtl_output(plink_output_filename)
-    
-    def geno_to_rqtl_function(self):
-
-        #TODO: Need to figure out why some genofiles have the wrong format and don't convert properly
 
+    def geno_to_rqtl_function(self):        # TODO: Need to figure out why some genofiles have the wrong format and don't convert properly
         print("Adding a function to the R environment")
         ro.r("""
+           trim <- function( x ) { gsub("(^[[:space:]]+|[[:space:]]+$)", "", x) }
+
            getGenoCode <- function(header, name = 'unk'){
              mat = which(unlist(lapply(header,function(x){ length(grep(paste('@',name,sep=''), x)) })) == 1)
-             return(strsplit(header[mat],'')[[1]][6])
+             return(trim(strsplit(header[mat],':')[[1]][2]))
            }
 
            GENOtoCSVR <- function(genotypes = 'BXD.geno', out = 'cross.csvr', phenotype = NULL, sex = NULL, verbose = FALSE){
              header = readLines(genotypes, 40)
              toskip = which(unlist(lapply(header, function(x){ length(grep("Chr\t", x)) })) == 1)-1                            # Major hack to skip the geno headers
+             
              genocodes <- c(getGenoCode(header, 'mat'), getGenoCode(header, 'het'), getGenoCode(header, 'pat'))                # Get the genotype codes 
-
+             type <- getGenoCode(header, 'type')
              genodata <- read.csv(genotypes, sep='\t', skip=toskip, header=TRUE, na.strings=getGenoCode(header,'unk'), colClasses='character', comment.char = '#')
-             cat('Genodata:', toskip, " ", dim(genodata), '\n')
+             cat('Genodata:', toskip, " ", dim(genodata), genocodes, '\n')
              if(is.null(phenotype)) phenotype <- runif((ncol(genodata)-4))                                                     # If there isn't a phenotype, generate a random one
              if(is.null(sex)) sex <- rep('m', (ncol(genodata)-4))                                                              # If there isn't a sex phenotype, treat all as males
              outCSVR <- rbind(c('Pheno', '', '', phenotype),                                                                   # Phenotype
@@ -254,29 +264,31 @@ class MarkerRegression(object):
                               cbind(genodata[,c('Locus','Chr', 'cM')], genodata[, 5:ncol(genodata)]))                          # Genotypes
              write.table(outCSVR, file = out, row.names=FALSE, col.names=FALSE,quote=FALSE, sep=',')                           # Save it to a file
              require(qtl)
-             return(read.cross(file=out, 'csvr', genotypes=genocodes))                                                  # Load it using R/qtl read.cross  
+             cross = read.cross(file=out, 'csvr', genotypes=genocodes)
+             if(type == 'riset') cross <- convert2riself(cross)
+             return(cross)                                                  # Load it using R/qtl read.cross  
           }
         """)
     
     def run_rqtl_geno(self):
-
         print("Calling R/qtl from python")
 
-        #TODO: Need to get this working for other groups/inbred sets, calculating file on the fly
         self.geno_to_rqtl_function()
 
         ## 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
-        
+        r_sum           = ro.r["sum"]             # Map the ncol function
+
         print(r_library("qtl"))                     # Load R/qtl
-  
+
         ## Get pointers to some R/qtl functions
         scanone         = ro.r["scanone"]               # Map the scanone function
+        scantwo         = ro.r["scantwo"]               # Map the scantwo function
         calc_genoprob   = ro.r["calc.genoprob"]         # Map the calc.genoprob function
         read_cross      = ro.r["read.cross"]            # Map the read.cross function
         write_cross     = ro.r["write.cross"]           # Map the write.cross function
-        GENOtoCSVR      = ro.r["GENOtoCSVR"]            # Map the write.cross function
+        GENOtoCSVR      = ro.r["GENOtoCSVR"]            # Map the GENOtoCSVR function
 
         genofilelocation  = webqtlConfig.HTMLPATH + "genotypes/" + self.dataset.group.name + ".geno"
         crossfilelocation = webqtlConfig.HTMLPATH + "genotypes/" + self.dataset.group.name + ".cross"
@@ -289,55 +301,62 @@ class MarkerRegression(object):
             cross_object = calc_genoprob(cross_object)
         else:
             cross_object = calc_genoprob(cross_object, step=1, stepwidth="max")
-       
-        # Add the phenotype
-        cross_object = self.add_phenotype(cross_object, self.sanitize_rqtl_phenotype())
+
+        cross_object = self.add_phenotype(cross_object, self.sanitize_rqtl_phenotype())                 # Add the phenotype
 
         # for debug: write_cross(cross_object, "csvr", "test.csvr")
 
         # Scan for QTLs
-        if(self.control.replace(" ", "") != ""):
-            covar = self.create_covariates(cross_object)
-            result_data_frame = scanone(cross_object, pheno = "the_pheno", addcovar = covar)
-        else:
-            result_data_frame = scanone(cross_object, pheno = "the_pheno")
+        covar = self.create_covariates(cross_object)
 
-        if int(self.num_perm) > 0:
-	     # Do permutation
-            if(self.control.replace(" ", "") != ""):
-                covar = self.create_covariates(cross_object)
-                perm_data_frame = scanone(cross_object, pheno_col = "the_pheno", addcovar = covar, n_perm=int(self.num_perm))
+        if self.pair_scan:
+            if(r_sum(covar)[0] > 0):
+                print("Using covariate"); result_data_frame = scantwo(cross_object, pheno = "the_pheno", addcovar = covar)
+            else:
+                print("No covariates"); result_data_frame = scantwo(cross_object, pheno = "the_pheno")
+ 
+            print("pair scan results:", result_data_frame)
+
+            return 0
+        else:
+            if(r_sum(covar)[0] > 0):
+                print("Using covariate"); result_data_frame = scanone(cross_object, pheno = "the_pheno", addcovar = covar)
             else:
-                perm_data_frame = scanone(cross_object, pheno_col = "the_pheno", n_perm=int(self.num_perm))
+                print("No covariates"); result_data_frame = scanone(cross_object, pheno = "the_pheno")
 
-        self.suggestive, self.significant = self.process_rqtl_perm_results(perm_data_frame)
+            if int(self.num_perm) > 0:                                                                      # Do permutation (if requested by user)
+                if(r_sum(covar)[0] > 0):
+                    perm_data_frame = scanone(cross_object, pheno_col = "the_pheno", addcovar = covar, n_perm = int(self.num_perm))
+                else:
+                    perm_data_frame = scanone(cross_object, pheno_col = "the_pheno", n_perm = int(self.num_perm))
 
-        qtl_results = self.process_rqtl_results(result_data_frame)
+                self.process_rqtl_perm_results(perm_data_frame)                                             # Functions that sets the thresholds for the webinterface
 
-        return qtl_results
+            return self.process_rqtl_results(result_data_frame)
 
     def add_phenotype(self, 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(self, cross):
         ro.globalenv["the_cross"] = cross
-        ro.r('genotypes <- pull.geno(the_cross)')       # Get genotype matrix
-        userinputS = self.control.replace(" ", "").split(",")     # TODO sanitize user input !!! never trust a user
-        covariate_names = ', '.join('"{0}"'.format(w) for w in userinputS) 
-        print(covariate_names)
-        ro.r('covariates <- genotypes[,c(' + covariate_names + ')]')         # get covariate matrix, 
-        print("COVARIATES:", ro.r["covariates"])
+        ro.r('genotypes <- pull.geno(the_cross)')                             # Get the genotype matrix
+        userinputS = self.control.replace(" ", "").split(",")                 # TODO sanitize user input, Never Ever trust a user
+        covariate_names = ', '.join('"{0}"'.format(w) for w in userinputS)
+        print("Marker names of selected covariates:", covariate_names)
+        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
+        print("R/qtl matrix of covariates:", ro.r["covariates"])
         return ro.r["covariates"]
 
     def sanitize_rqtl_phenotype(self):
         pheno_as_string = "c("
-        null_pos = []
         for i, val in enumerate(self.vals):
             if val == "x":
-                null_pos.append(i)
                 if i < (len(self.vals) - 1):
                     pheno_as_string +=  "NA,"
                 else:
@@ -350,35 +369,33 @@ class MarkerRegression(object):
         pheno_as_string += ")"
         return pheno_as_string
 
-    def process_rqtl_results(self, result):
-        #TODO how to make this a one liner and not copy the stuff in a loop
+    def process_rqtl_results(self, result):        # TODO: how to make this a one liner and not copy the stuff in a loop
         qtl_results = []
-        
+
         output = [tuple([result[j][i] for j in range(result.ncol)]) for i in range(result.nrow)]
-        print("output", output)
-        
+        print("R/qtl scanone output:", output)
+
         for i, line in enumerate(result.iter_row()):
             marker = {}
             marker['name'] = result.rownames[i]
             marker['chr'] = output[i][0]
             marker['Mb'] = output[i][1]
             marker['lod_score'] = output[i][2]
-            
             qtl_results.append(marker)
-            
+
         return qtl_results
-    
+
     def process_rqtl_perm_results(self, results):
         perm_vals = []
         for line in str(results).split("\n")[1:(int(self.num_perm)+1)]:
-            print("line:", line.split())
+            print("R/qtl permutation line:", line.split())
             perm_vals.append(float(line.split()[1]))
-        
+
         self.suggestive = np.percentile(np.array(perm_vals), 67)
         self.significant = np.percentile(np.array(perm_vals), 95)
-        
+
         return self.suggestive, self.significant
-            
+
 
     def run_plink(self):
     
diff --git a/wqflask/wqflask/show_trait/show_trait.py b/wqflask/wqflask/show_trait/show_trait.py
index 62411640..11fd2644 100755
--- a/wqflask/wqflask/show_trait/show_trait.py
+++ b/wqflask/wqflask/show_trait/show_trait.py
@@ -1310,6 +1310,9 @@ def get_nearest_marker(this_trait, this_db):
     result = g.db.execute(query).fetchall()
     print("result:", result)
 
-    return result[0][0], result[1][0]
+    if result == []:
+        return "", ""
+    else:
+        return result[0][0], result[1][0]
     
     
diff --git a/wqflask/wqflask/templates/show_trait_mapping_tools.html b/wqflask/wqflask/templates/show_trait_mapping_tools.html
index d6c7b18e..946dab49 100755
--- a/wqflask/wqflask/templates/show_trait_mapping_tools.html
+++ b/wqflask/wqflask/templates/show_trait_mapping_tools.html
@@ -80,19 +80,31 @@
                 </div>
                 
                 <div class="tab-pane" id="pylmm">
+                    
                     <div style="padding: 20px" class="form-horizontal">
                         <div class="mapping_method_fields form-group">
-                            <label for="mapping_permutations" class="col-xs-1 control-label">Permutations</label>
-                            <div style="margin-left: 20px;" class="col-xs-2 controls">
-                                <input name="num_perm_pylmm" value="0" type="text" class="form-control">
+                            <label for="mapping_permutations" class="col-xs-2 control-label">Permutations</label>
+                            <div style="margin-left: 20px;" class="col-xs-4 controls">
+                                <input name="num_pylmm" value="2000" type="text" class="form-control">
                             </div>
                         </div>
                         <div id="permutations_alert" class="alert alert-error alert-warning" style="display:none;">
                             Please be aware that permutations can take a very long time (~20 minutes for 500 permutations)
                         </div>
                         <div class="mapping_method_fields form-group">
-                            <label class="col-xs-1 control-label">Manhattan Plot</label>
-                            <div style="margin-left:20px;" class="col-xs-4 controls">                      
+                            <label for="control_for" class="col-xs-2 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_pylmm" value="{{ nearest_marker1+","+nearest_marker2 }}" type="text" />
+                                {% else %}
+                                <input name="control_pylmm" value="" type="text" />
+                                {% endif %}
+                            </div>
+                        </div>
+
+                        <div class="mapping_method_fields form-group">
+                            <label style="text-align:left;" class="col-xs-12 control-label">Manhattan Plot</label>
+                            <div class="col-xs-12 controls">                      
                                 <label class="radio-inline">
                                     <input type="radio" name="manhattan_plot_pylmm" value="true">
                                     Yes
@@ -103,39 +115,30 @@
                                </label>
                             </div>
                         </div>
-                        <!--<div class="form-group" id="suggestive">
-                            <label for="suggestive_reaper" class="control-label" title="Control Locus">LOD score greater than: </label>
-                            <div class="controls">
-                                <input name="suggestive_reaper" type="text" value="0" class="form-control">
+                        <div class="form-group">
+                            <div style="padding-left:15px;" class="controls">
+                                <button id="pylmm_compute" class="btn submit_special btn-primary" data-url="/marker_regression" title="Compute Marker Regression">
+                                    <i class="icon-ok-circle icon-white"></i> Open Mapping Tool
+                                </button>
                             </div>
-                        </div>-->
-                    </div>
-                    
-                    <div class="form-group">
-                        <label for="marker_regression_submit" class="col-xs-1 control-label"></label>
-                        <div style="margin-left:20px;" class="col-xs-4 controls">
-                            <button id="pylmm_compute" class="btn submit_special btn-primary" data-url="/marker_regression" title="Compute Marker Regression">
-                                Compute
-                            </button>
                         </div>
                     </div>
-
                 </div>
                 <div class="tab-pane" id="rqtl_geno">
                     
                     <div style="padding: 20px" class="form-horizontal">
                         <div class="mapping_method_fields form-group">
-                            <label for="mapping_permutations" class="col-xs-1 control-label">Permutations</label>
-                            <div style="margin-left: 20px;" class="col-xs-2 controls">
-                                <input name="num_perm_rqtl_geno" value="500" type="text" class="form-control">
+                            <label for="mapping_permutations" class="col-xs-2 control-label">Permutations</label>
+                            <div style="margin-left: 20px;" class="col-xs-4 controls">
+                                <input name="num_perm_rqtl_geno" value="2000" type="text" class="form-control">
                             </div>
                         </div>
                         <div id="permutations_alert" class="alert alert-error alert-warning" style="display:none;">
                             Please be aware that permutations can take a very long time (~20 minutes for 500 permutations)
                         </div>
                         <div class="mapping_method_fields form-group">
-                            <label for="mapping_permutations" class="col-xs-1 control-label">Control&nbsp;for</label>
-                            <div style="margin-left: 20px;" class="col-xs-2 controls">
+                            <label for="control_for" class="col-xs-2 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_rqtl_geno" value="{{ nearest_marker1+","+nearest_marker2 }}" type="text" />
                                 {% else %}
@@ -144,8 +147,8 @@
                             </div>
                         </div>
                         <div class="mapping_method_fields form-group">
-                            <label class="col-xs-1 control-label">Manhattan&nbsp;Plot</label>
-                            <div style="margin-left:35px;" class="col-xs-4 controls">                      
+                            <label style="text-align:left;" class="col-xs-12 control-label">Manhattan Plot</label>
+                            <div class="col-xs-12 controls">                      
                                 <label class="radio-inline">
                                     <input type="radio" name="manhattan_plot_rqtl" value="true">
                                     Yes
@@ -156,17 +159,27 @@
                                </label>
                             </div>
                         </div>
-                    </div>
-                    
-                    <div class="form-group">
-                        <label for="marker_regression_submit" class="col-xs-1 control-label"></label>
-                        <div style="margin-left:20px;" class="col-xs-4 controls">
-                            <button id="rqtl_geno_compute" class="btn submit_special btn-primary" data-url="/marker_regression" title="Compute Marker Regression">
-                                <i class="icon-ok-circle icon-white"></i> Compute
-                            </button>
+                        <div class="mapping_method_fields form-group">
+                            <label style="text-align:left;" class="col-xs-12 control-label">Pair Scan</label>
+                            <div class="col-xs-12 controls">                      
+                                <label class="radio-inline">
+                                    <input type="radio" name="pair_scan" value="true">
+                                    Yes
+                                </label>
+                                <label class="radio-inline">
+                                    <input type="radio" name="pair_scan" value="false" checked="">
+                                    No
+                               </label>
+                            </div>
+                        </div>
+                        <div class="form-group">
+                            <div style="padding-left:15px;" class="controls">
+                                <button id="rqtl_geno_compute" class="btn submit_special btn-primary" data-url="/marker_regression" title="Compute Marker Regression">
+                                    <i class="icon-ok-circle icon-white"></i> Open Mapping Tool
+                                </button>
+                            </div>
                         </div>
                     </div>
-       
                 </div>
                 
                 {% if dataset.group.species == 'human' %}
diff --git a/wqflask/wqflask/views.py b/wqflask/wqflask/views.py
index 9110c5a1..e60852d5 100755
--- a/wqflask/wqflask/views.py
+++ b/wqflask/wqflask/views.py
@@ -285,7 +285,8 @@ def marker_regression_page():
         'maf',
         'manhattan_plot',
         'control_marker',
-        'control_marker_db'
+        'control_marker_db',
+        'pair_scan'
     )
 
     start_vars = {}