Cleaning up some of the code, updating output and adding/changing comments

1 files changed, 18 insertions, 30 deletions

diff --git a/wqflask/wqflask/marker_regression/marker_regression.py b/wqflask/wqflask/marker_regression/marker_regression.py
index 64d3ef3d..5ddae0a1 100755
--- a/wqflask/wqflask/marker_regression/marker_regression.py
+++ b/wqflask/wqflask/marker_regression/marker_regression.py
@@ -228,11 +228,8 @@ 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("""
            getGenoCode <- function(header, name = 'unk'){
@@ -259,24 +256,22 @@ class MarkerRegression(object):
         """)
     
     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
-        
+
         print(r_library("qtl"))                     # Load R/qtl
-  
+
         ## Get pointers to some R/qtl functions
         scanone         = ro.r["scanone"]               # Map the scanone 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,9 +284,8 @@ 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")
 
@@ -302,26 +296,22 @@ class MarkerRegression(object):
         else:
             result_data_frame = scanone(cross_object, pheno = "the_pheno")
 
-        if int(self.num_perm) > 0:
-	     # Do permutation
+        if int(self.num_perm) > 0:                                                                  # Do permutation (if requested by user)
             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))
             else:
                 perm_data_frame = scanone(cross_object, pheno_col = "the_pheno", n_perm=int(self.num_perm))
 
-        self.suggestive, self.significant = self.process_rqtl_perm_results(perm_data_frame)
-
-        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
@@ -350,35 +340,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):