diff options
| -rwxr-xr-x | wqflask/wqflask/marker_regression/marker_regression.py | 161 |
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()) |