diff options
Diffstat (limited to 'wqflask')
| -rwxr-xr-x | wqflask/wqflask/marker_regression/marker_regression.py | 121 | ||||
| -rwxr-xr-x | wqflask/wqflask/show_trait/show_trait.py | 5 | ||||
| -rwxr-xr-x | wqflask/wqflask/templates/show_trait_mapping_tools.html | 83 | ||||
| -rwxr-xr-x | wqflask/wqflask/views.py | 3 |
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 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 for</label> - <div style="margin-left: 20px;" class="col-xs-2 controls"> + <label for="control_for" class="col-xs-2 control-label">Control 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 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 = {} |