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