aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rwxr-xr-xwqflask/wqflask/marker_regression/marker_regression.py161
1 files changed, 17 insertions, 144 deletions
diff --git a/wqflask/wqflask/marker_regression/marker_regression.py b/wqflask/wqflask/marker_regression/marker_regression.py
index 99f23676..140da0c5 100755
--- a/wqflask/wqflask/marker_regression/marker_regression.py
+++ b/wqflask/wqflask/marker_regression/marker_regression.py
@@ -231,7 +231,7 @@ class MarkerRegression(object):
def run_rqtl_geno(self):
- print("CORRECT FUNCTION BEING CALLED")
+ print("Calling R/qtl from python")
## Get pointers to some common R functions
r_library = ro.r["library"] # Map the library function
@@ -255,14 +255,26 @@ class MarkerRegression(object):
# Add the phenotype
cross_object = self.add_phenotype(cross_object, self.sanitize_rqtl_phenotype())
+
# for debug: write_cross(cross_object, "csvr", "test.csvr")
- # Scan for QTLs
- # TODO: use addcovar, and permutations
+
+ # Scan for QTLs
if(self.control.replace(" ", "") != ""):
- covar = self.create_covariatesShort(cross_object)
+ 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")
+
+ if int(self.num_perm) > 0:
+ # Do permutation
+ if(self.control.replace(" ", "") != ""):
+ covar = self.create_covariatesShort(cross_object)
+ perm_data_frame = scanone(cross_object, pheno = "the_pheno", addcovar = covar, n_perm=int(self.num_perm))
+ else:
+ perm_data_frame = scanone(cross_object, pheno = "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)
return qtl_results
@@ -273,7 +285,7 @@ class MarkerRegression(object):
return ro.r["the_cross"]
- def create_covariatesShort(self, 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
@@ -283,46 +295,6 @@ class MarkerRegression(object):
print("COVARIATES:", ro.r["covariates"])
return ro.r["covariates"]
-
- def create_covariates(self):
- if self.control != "":
- control_markers = self.control.split(",")
- control_string = ""
- for i, control in enumerate(control_markers):
- control_trait = GeneralTrait(name=str(control), dataset_name=str(self.dataset.group.name + "Geno"))
- control_vals = []
- for sample in self.dataset.group.samplelist:
- if sample in control_trait.data:
- control_vals.append(control_trait.data[sample].value)
- else:
- control_vals.append("x")
- print("control_vals:", control_vals)
- control_as_string = "c("
- for j, val2 in enumerate(control_vals):
- if val2 == "x":
- if j < (len(control_vals) - 1):
- control_as_string += "NA,"
- else:
- control_as_string += "NA"
- else:
- if j < (len(control_vals) - 1):
- control_as_string += str(val2) + ","
- else:
- control_as_string += str(val2)
- #if i < (len(control_vals) - 1):
- # control_as_string += str(new_val2) + ","
- #else:
- # control_as_string += str(new_val2)
- control_as_string += ")"
- print("control_as_string:", control_as_string)
- if i < (len(control_markers)-1):
- control_string += control_as_string + ","
- else:
- control_string += control_as_string
-
- robjects.r('covariates <- cbind( '+ control_string +')')
-
-
def sanitize_rqtl_phenotype(self):
pheno_as_string = "c("
null_pos = []
@@ -358,107 +330,8 @@ class MarkerRegression(object):
qtl_results.append(marker)
return qtl_results
-
- def run_rqtl_geno_old(self):
- robjects.packages.importr("qtl")
- robjects.r('the_cross <- read.cross(format="csvr", dir="/home/zas1024/PLINK2RQTL/test", file="BXD.csvr")')
- if self.manhattan_plot:
- robjects.r('the_cross <- calc.genoprob(the_cross)')
- else:
- robjects.r('the_cross <- calc.genoprob(the_cross, step=1, stepwidth="max")')
- 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:
- pheno_as_string += "NA"
- else:
- if i < (len(self.vals) - 1):
- pheno_as_string += str(val) + ","
- else:
- pheno_as_string += str(val)
-
- pheno_as_string += ")"
-
- robjects.r('the_cross$pheno <- cbind(pull.pheno(the_cross), the_pheno = '+ pheno_as_string +')')
-
- print("self.control:", self.control)
-
- if self.control != "":
- control_markers = self.control.split(",")
- control_string = ""
- for i, control in enumerate(control_markers):
- control_trait = GeneralTrait(name=str(control), dataset_name=str(self.dataset.group.name + "Geno"))
- control_vals = []
- for sample in self.dataset.group.samplelist:
- if sample in control_trait.data:
- control_vals.append(control_trait.data[sample].value)
- else:
- control_vals.append("x")
- print("control_vals:", control_vals)
- control_as_string = "c("
- for j, val2 in enumerate(control_vals):
- if val2 == "x":
- if j < (len(control_vals) - 1):
- control_as_string += "NA,"
- else:
- control_as_string += "NA"
- else:
- if j < (len(control_vals) - 1):
- control_as_string += str(val2) + ","
- else:
- control_as_string += str(val2)
- #if i < (len(control_vals) - 1):
- # control_as_string += str(new_val2) + ","
- #else:
- # control_as_string += str(new_val2)
- control_as_string += ")"
- print("control_as_string:", control_as_string)
- if i < (len(control_markers)-1):
- control_string += control_as_string + ","
- else:
- control_string += control_as_string
-
- robjects.r('covariates <- cbind( '+ control_string +')')
- #Scan for QTL
- r_string = 'scanone(the_cross, pheno.col="the_pheno", addcovar=covariates)'
- #r_string = 'scanone(the_cross, pheno.col="the_pheno", n.cluster=16, n.perm='+self.num_perm+', addcovar=covariates, intcovar=covariates[,'+ str(len(control_markers)) +'])'
- print("r_string:", r_string)
-
- # extract the profile
-
- if int(self.num_perm) > 0:
- # Do permutation
- r_string = 'scanone(the_cross, pheno.col="the_pheno", n.cluster=16, addcovar=covariates, intcovar=covariates[,'+ str(len(control_markers)) +'])'
- # Extract the tresholds
- thresholds = robjects.r(r_string)
- self.suggestive, self.significant = self.process_rqtl_perm_results(thresholds)
-
- #r_string = 'scanone(the_cross, pheno.col='+pheno_as_string+', addcovar='+control_as_string+')'
-
- else:
- #r_string = 'scanone(the_cross, pheno.col='+pheno_as_string+', n.perm='+self.num_perm+')'
- r_string = 'scanone(the_cross, pheno.col="the_pheno", n.cluster=16, n.perm='+self.num_perm+')'
- if self.num_perm.isdigit() and int(self.num_perm) > 0:
- results = robjects.r(r_string)
- self.suggestive, self.significant = self.process_rqtl_perm_results(results)
- r_string = 'scanone(the_cross, pheno.col="the_pheno", n.cluster=16)'
-
- print("r_string:", r_string)
- result_data_frame = robjects.r(r_string)
- #print("results:", result_data_frame)
-
- print("before process rqtl results")
-
- qtl_results = self.process_rqtl_results(result_data_frame)
-
- 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())