diff options
7 files changed, 177 insertions, 48 deletions
diff --git a/wqflask/wqflask/marker_regression/display_mapping_results.py b/wqflask/wqflask/marker_regression/display_mapping_results.py index 2a53b60e..f70bc555 100644 --- a/wqflask/wqflask/marker_regression/display_mapping_results.py +++ b/wqflask/wqflask/marker_regression/display_mapping_results.py @@ -265,14 +265,12 @@ class DisplayMappingResults(object): else: self.colorCollection = [self.LRS_COLOR] + self.dataset.group.genofile = self.genofile_string.split(":")[0] if self.mapping_method == "reaper" and self.manhattan_plot != True: self.genotype = self.dataset.group.read_genotype_file(use_reaper=True) else: self.genotype = self.dataset.group.read_genotype_file() - #if self.mapping_method == "rqtl_geno" and self.genotype.filler == True: - # self.genotype = self.genotype.read_rdata_output(self.qtlresults) - #Darwing Options try: if self.selectedChr > -1: @@ -1761,9 +1759,9 @@ class DisplayMappingResults(object): break if all_int: - max_lrs_width = canvas.stringWidth("%d" % LRS_LOD_Max, font=LRSScaleFont) + 30 + max_lrs_width = canvas.stringWidth("%d" % LRS_LOD_Max, font=LRSScaleFont) + 40 else: - max_lrs_width = canvas.stringWidth("%2.1f" % LRS_LOD_Max, font=LRSScaleFont) + 20 + max_lrs_width = canvas.stringWidth("%2.1f" % LRS_LOD_Max, font=LRSScaleFont) + 30 #draw the "LRS" or "LOD" string to the left of the axis canvas.drawString(self.LRS_LOD, xLeftOffset - max_lrs_width - 15*(zoom-1), \ @@ -1899,13 +1897,16 @@ class DisplayMappingResults(object): this_chr = str(self.ChrList[self.selectedChr][1]+1) if self.selectedChr == -1 or str(qtlresult['chr']) == this_chr: - if self.plotScale != "physic" and self.genotype.filler == True: - if self.selectedChr != -1: - start_cm = self.genotype[self.selectedChr - 1][0].cM - Xc = startPosX + (qtlresult['Mb'] - start_cm)*plotXScale - else: - start_cm = self.genotype[previous_chr_as_int][0].cM - Xc = startPosX + ((qtlresult['Mb']-start_cm-startMb)*plotXScale)*(((qtlresult['Mb']-start_cm-startMb)*plotXScale)/((qtlresult['Mb']-start_cm-startMb+self.GraphInterval)*plotXScale)) + if self.plotScale != "physic" and self.mapping_method == "reaper" and not self.manhattan_plot: + Xc = startPosX + (qtlresult['cM']-startMb)*plotXScale + if hasattr(self.genotype, "filler"): + if self.genotype.filler: + if self.selectedChr != -1: + start_cm = self.genotype[self.selectedChr - 1][0].cM + Xc = startPosX + (qtlresult['Mb'] - start_cm)*plotXScale + else: + start_cm = self.genotype[previous_chr_as_int][0].cM + Xc = startPosX + ((qtlresult['Mb']-start_cm-startMb)*plotXScale)*(((qtlresult['Mb']-start_cm-startMb)*plotXScale)/((qtlresult['Mb']-start_cm-startMb+self.GraphInterval)*plotXScale)) else: Xc = startPosX + (qtlresult['Mb']-startMb)*plotXScale diff --git a/wqflask/wqflask/marker_regression/rqtl_mapping.py b/wqflask/wqflask/marker_regression/rqtl_mapping.py index c1a56787..8c294460 100644 --- a/wqflask/wqflask/marker_regression/rqtl_mapping.py +++ b/wqflask/wqflask/marker_regression/rqtl_mapping.py @@ -11,7 +11,7 @@ from utility.tools import locate, TEMPDIR import utility.logger logger = utility.logger.getLogger(__name__ ) -def run_rqtl_geno(vals, samples, dataset, method, model, permCheck, num_perm, perm_strata_list, do_control, control_marker, manhattan_plot, pair_scan, cofactors): +def run_rqtl_geno(vals, samples, dataset, mapping_scale, method, model, permCheck, num_perm, perm_strata_list, do_control, control_marker, manhattan_plot, pair_scan, cofactors): ## 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 @@ -33,7 +33,13 @@ def run_rqtl_geno(vals, samples, dataset, method, model, permCheck, num_perm, pe # genofilelocation = locate(crossname + ".RData", "genotype/rdata") # cross_object = read_cross_from_rdata(genofilelocation) # Map the local GENOtoCSVR function #except: - generate_cross_from_geno(dataset) + + if mapping_scale == "morgan": + scale_units = "cM" + else: + scale_units = "Mb" + + generate_cross_from_geno(dataset, scale_units) GENOtoCSVR = ro.r["GENOtoCSVR"] # Map the local GENOtoCSVR function crossfilelocation = TMPDIR + crossname + ".cross" if dataset.group.genofile: @@ -47,6 +53,8 @@ def run_rqtl_geno(vals, samples, dataset, method, model, permCheck, num_perm, pe else: cross_object = calc_genoprob(cross_object, step=1, stepwidth="max") + logger.debug("VAL LEN:", len(vals)) + pheno_string = sanitize_rqtl_phenotype(vals) cross_object = add_phenotype(cross_object, pheno_string, "the_pheno") # Add the phenotype @@ -94,11 +102,9 @@ def run_rqtl_geno(vals, samples, dataset, method, model, permCheck, num_perm, pe perm_data_frame = scanone(cross_object, pheno_col = "the_pheno", n_perm = num_perm, model=model, method=method) perm_output, suggestive, significant = process_rqtl_perm_results(num_perm, perm_data_frame) # Functions that sets the thresholds for the webinterface - the_scale = check_mapping_scale(genofilelocation) - return perm_output, suggestive, significant, process_rqtl_results(result_data_frame, dataset.group.species), the_scale + return perm_output, suggestive, significant, process_rqtl_results(result_data_frame, dataset.group.species) else: - the_scale = check_mapping_scale(genofilelocation) - return process_rqtl_results(result_data_frame, dataset.group.species), the_scale + return process_rqtl_results(result_data_frame, dataset.group.species) def generate_cross_from_rdata(dataset): rdata_location = locate(dataset.group.name + ".RData", "genotype/rdata") @@ -110,7 +116,7 @@ def generate_cross_from_rdata(dataset): } """ % (rdata_location)) -def generate_cross_from_geno(dataset): # TODO: Need to figure out why some genofiles have the wrong format and don't convert properly +def generate_cross_from_geno(dataset, scale_units): # TODO: Need to figure out why some genofiles have the wrong format and don't convert properly ro.r(""" trim <- function( x ) { gsub("(^[[:space:]]+|[[:space:]]+$)", "", x) } @@ -127,21 +133,23 @@ def generate_cross_from_geno(dataset): # TODO: Need to figure out why som type <- getGenoCode(header, 'type') if(type == '4-way'){ genocodes <- c('1','2','3','4') + genodata <- read.csv(genotypes, sep='\t', skip=toskip, header=TRUE, na.strings=getGenoCode(header,'unk'), colClasses='character', comment.char = '#', crosstype="4way") } else { genocodes <- c(getGenoCode(header, 'mat'), getGenoCode(header, 'het'), getGenoCode(header, 'pat')) # Get the genotype codes + genodata <- read.csv(genotypes, sep='\t', skip=toskip, header=TRUE, na.strings=getGenoCode(header,'unk'), colClasses='character', comment.char = '#') } - genodata <- read.csv(genotypes, sep='\t', skip=toskip, header=TRUE, na.strings=getGenoCode(header,'unk'), colClasses='character', comment.char = '#') 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 c('sex', '', '', sex), # Sex phenotype for the mice - cbind(genodata[,c('Locus','Chr', 'cM')], genodata[, 5:ncol(genodata)])) # Genotypes + cbind(genodata[,c('Locus','Chr', '%s')], 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) if(type == '4-way'){ cat('Loading in as 4-WAY\n') - cross = read.cross(file=out, 'csvr', genotypes=genocodes, crosstype="4way", convertXdata=FALSE) # Load the created cross file using R/qtl read.cross + cross = read.cross(file=out, 'csvr', genotypes=genocodes) + #cross = read.cross(file=out, 'csvr', genotypes=genocodes, crosstype="4way", convertXdata=FALSE) # Load the created cross file using R/qtl read.cross }else{ cat('Loading in as normal\n') cross = read.cross(file=out, 'csvr', genotypes=genocodes) # Load the created cross file using R/qtl read.cross @@ -152,7 +160,7 @@ def generate_cross_from_geno(dataset): # TODO: Need to figure out why som } return(cross) } - """ % (dataset.group.genofile)) + """ % (dataset.group.genofile, scale_units)) def add_perm_strata(cross, perm_strata): col_string = 'c("the_strata")' @@ -300,20 +308,4 @@ def process_rqtl_results(result, species_name): # TODO: how to make this marker['lod_score'] = output[i][2] qtl_results.append(marker) - return qtl_results - -def check_mapping_scale(genofile_location): - scale = "physic" - with open(genofile_location, "r") as geno_fh: - for line in geno_fh: - if line[0] == "@" or line[0] == "#": - - if "@scale" in line: - scale = line.split(":")[1].strip() - break - else: - continue - else: - break - - return scale
\ No newline at end of file + return qtl_results
\ No newline at end of file diff --git a/wqflask/wqflask/marker_regression/run_mapping.py b/wqflask/wqflask/marker_regression/run_mapping.py index 589be702..7449d8ce 100644 --- a/wqflask/wqflask/marker_regression/run_mapping.py +++ b/wqflask/wqflask/marker_regression/run_mapping.py @@ -156,6 +156,8 @@ class RunMapping(object): self.transform = "" self.score_type = "LRS" #ZS: LRS or LOD self.mapping_scale = "physic" + if "mapping_scale" in start_vars: + self.mapping_scale = start_vars['mapping_scale'] self.num_perm = 0 self.perm_output = [] self.bootstrap_results = [] @@ -255,9 +257,9 @@ class RunMapping(object): #if start_vars['pair_scan'] == "true": # self.pair_scan = True if self.permCheck and self.num_perm > 0: - self.perm_output, self.suggestive, self.significant, results, self.mapping_scale = rqtl_mapping.run_rqtl_geno(self.vals, self.samples, self.dataset, self.method, self.model, self.permCheck, self.num_perm, perm_strata, self.do_control, self.control_marker, self.manhattan_plot, self.pair_scan, self.covariates) + self.perm_output, self.suggestive, self.significant, results= rqtl_mapping.run_rqtl_geno(self.vals, self.samples, self.dataset, self.mapping_scale, self.method, self.model, self.permCheck, self.num_perm, perm_strata, self.do_control, self.control_marker, self.manhattan_plot, self.pair_scan, self.covariates) else: - results, self.mapping_scale = rqtl_mapping.run_rqtl_geno(self.vals, self.samples, self.dataset, self.method, self.model, self.permCheck, self.num_perm, perm_strata, self.do_control, self.control_marker, self.manhattan_plot, self.pair_scan, self.covariates) + results = rqtl_mapping.run_rqtl_geno(self.vals, self.samples, self.dataset, self.mapping_scale, self.method, self.model, self.permCheck, self.num_perm, perm_strata, self.do_control, self.control_marker, self.manhattan_plot, self.pair_scan, self.covariates) elif self.mapping_method == "reaper": if "startMb" in start_vars: #ZS: Check if first time page loaded, so it can default to ON if "additiveCheck" in start_vars: @@ -429,7 +431,7 @@ class RunMapping(object): with Bench("Trimming Markers for Table"): self.trimmed_markers = trim_markers_for_table(results) - chr_lengths = get_chr_lengths(self.mapping_scale, self.dataset, self.qtl_results) + chr_lengths = get_chr_lengths(self.mapping_scale, self.mapping_method, self.dataset, self.qtl_results) #ZS: For zooming into genome browser, need to pass chromosome name instead of number if self.dataset.group.species == "mouse": @@ -643,7 +645,7 @@ def geno_db_exists(this_dataset): except: return "False" -def get_chr_lengths(mapping_scale, dataset, qtl_results): +def get_chr_lengths(mapping_scale, mapping_method, dataset, qtl_results): chr_lengths = [] if mapping_scale == "physic": for i, the_chr in enumerate(dataset.species.chromosomes.chromosomes): @@ -666,8 +668,12 @@ def get_chr_lengths(mapping_scale, dataset, qtl_results): this_chr = chr_as_num highest_pos = 0 else: - if float(result['Mb']) > highest_pos: - highest_pos = float(result['Mb']) + if mapping_method == "reaper": + if float(result['cM']) > highest_pos: + highest_pos = float(result['cM']) + else: + if float(result['Mb']) > highest_pos: + highest_pos = float(result['Mb']) return chr_lengths diff --git a/wqflask/wqflask/show_trait/show_trait.py b/wqflask/wqflask/show_trait/show_trait.py index 8883e627..40e344b8 100644 --- a/wqflask/wqflask/show_trait/show_trait.py +++ b/wqflask/wqflask/show_trait/show_trait.py @@ -22,6 +22,7 @@ from base import webqtlConfig from base import webqtlCaseData from wqflask.show_trait.SampleList import SampleList from utility import webqtlUtil, Plot, Bunch, helper_functions +from utility.tools import locate_ignore_error from base.trait import GeneralTrait from base import data_set from db import webqtlDatabaseFunction @@ -170,6 +171,17 @@ class ShowTrait(object): self.genofiles = self.dataset.group.get_genofiles() + if "QTLReaper" or "R/qtl" in dataset.group.mapping_names: #ZS: No need to grab scales from .geno file unless it's using a mapping method that reads .geno files + if self.genofiles: + self.scales_in_geno = get_genotype_scales(self.genofiles) + else: + self.scales_in_geno = get_genotype_scales(self.dataset.group + ".geno") + + if len(self.scales_in_geno) < 2: + hddn['mapping_scale'] = self.scales_in_geno[self.scales_in_geno.keys()[0]][0] + else: + self.scales_in_geno = {} + self.has_num_cases = has_num_cases(self.this_trait) self.stats_table_width, self.trait_table_width = get_table_widths(self.sample_groups, self.has_num_cases) @@ -239,6 +251,7 @@ class ShowTrait(object): #hddn['control_marker'] = self.nearest_marker1+","+self.nearest_marker2 hddn['do_control'] = False hddn['maf'] = 0.05 + hddn['mapping_scale'] = "physic" hddn['compare_traits'] = [] hddn['export_data'] = "" hddn['export_format'] = "excel" @@ -251,6 +264,7 @@ class ShowTrait(object): short_description = short_description, unit_type = trait_units, dataset_type = self.dataset.type, + scales_in_geno = self.scales_in_geno, data_scale = self.dataset.data_scale, sample_group_types = self.sample_group_types, sample_lists = sample_lists, @@ -597,4 +611,77 @@ def get_categorical_variables(this_trait, sample_list): if num_distinct < 10: categorical_var_list.append(sample_list.attributes[attribute].name) - return categorical_var_list
\ No newline at end of file + return categorical_var_list + +def get_genotype_scales(genofiles): + geno_scales = {} + if type(genofiles) is list: + for the_file in genofiles: + file_location = the_file['location'] + geno_scales[file_location] = get_scales_from_genofile(file_location) + else: + geno_scales[genofiles] = get_scales_from_genofile(genofiles) + + return geno_scales + +def get_scales_from_genofile(file_location): + geno_path = locate_ignore_error(file_location, 'genotype') + + if not geno_path: #ZS: This is just to allow the code to run when + return [["physic", "Mb"]] + cm_and_mb_cols_exist = True + cm_column = None + mb_column = None + with open(geno_path, "r") as geno_fh: + for i, line in enumerate(geno_fh): + if line[0] == "#" or line[0] == "@": + if "@scale" in line: #ZS: If the scale is made explicit in the metadata, use that + scale = line.split(":")[1].strip() + if scale == "morgan": + return [["morgan", "cM"]] + else: + return [["physic", "Mb"]] + else: + continue + if line[:3] == "Chr": + first_marker_line = i + 1 + if line.split("\t")[2].strip() == "cM": + cm_column = 2 + elif line.split("\t")[3].strip() == "cM": + cm_column = 3 + if line.split("\t")[2].strip() == "Mb": + mb_column = 2 + elif line.split("\t")[3].strip() == "Mb": + mb_column = 3 + break + + #ZS: This attempts to check whether the cM and Mb columns are 'real', since some .geno files have one column be a copy of the other column, or have one column that is all 0s + cm_all_zero = True + mb_all_zero = True + cm_mb_all_equal = True + for i, line in enumerate(geno_fh): + if first_marker_line <= i < first_marker_line + 10: #ZS: I'm assuming there won't be more than 10 markers where the position is listed as 0 + if cm_column: + cm_val = line.split("\t")[cm_column].strip() + if cm_val != "0": + cm_all_zero = False + if mb_column: + mb_val = line.split("\t")[mb_column].strip() + if mb_val != "0": + mb_all_zero = False + if cm_column and mb_column: + if cm_val != mb_val: + cm_mb_all_equal = False + else: + if i > first_marker_line + 10: + break + + #ZS: This assumes that both won't be all zero, since if that's the case mapping shouldn't be an option to begin with + if mb_all_zero: + return [["morgan", "cM"]] + elif cm_mb_all_equal: + return [["physic", "Mb"]] + elif cm_and_mb_cols_exist: + return [["physic", "Mb"], ["morgan", "cM"]] + else: + return [["physic", "Mb"]] diff --git a/wqflask/wqflask/static/new/javascript/show_trait_mapping_tools.js b/wqflask/wqflask/static/new/javascript/show_trait_mapping_tools.js index 478ed87e..7176a0da 100644 --- a/wqflask/wqflask/static/new/javascript/show_trait_mapping_tools.js +++ b/wqflask/wqflask/static/new/javascript/show_trait_mapping_tools.js @@ -165,6 +165,7 @@ url = "/loading"; $('input[name=method]').val("rqtl_geno"); $('input[name=selected_chr]').val($('#chr_rqtl_geno').val()); + $('input[name=mapping_scale]').val($('#scale_rqtl_geno').val()); $('input[name=genofile]').val($('#genofile_rqtl_geno').val()); $('input[name=num_perm]').val($('input[name=num_perm_rqtl_geno]').val()); $('input[name=categorical_vars]').val(js_data.categorical_vars) @@ -210,6 +211,7 @@ url = "/loading"; $('input[name=method]').val("reaper"); $('input[name=selected_chr]').val($('#chr_reaper').val()); + $('input[name=mapping_scale]').val($('#scale_reaper').val()); $('input[name=genofile]').val($('#genofile_reaper').val()); $('input[name=num_perm]').val($('input[name=num_perm_reaper]').val()); $('input[name=control_marker]').val($('input[name=control_reaper]').val()); @@ -289,4 +291,21 @@ return toggle_enable_disable("#suggestive_lrs"); }); + $('#genofile_rqtl_geno').change(function() { + geno_location = $(this).children("option:selected").val().split(":")[0] + $('#scale_rqtl_geno').empty() + the_scales = js_data.scales_in_geno[geno_location] + for (var i = 0; i < the_scales.length; i++){ + $('#scale_rqtl_geno').append($("<option></option>").attr("value", the_scales[i][0]).text(the_scales[i][1])); + } + }); + $('#genofile_reaper').change(function() { + geno_location = $(this).children("option:selected").val().split(":")[0] + $('#scale_reaper').empty() + the_scales = js_data.scales_in_geno[geno_location] + for (var i = 0; i < the_scales.length; i++){ + $('#scale_reaper').append($("<option></option>").attr("value", the_scales[i][0]).text(the_scales[i][1])); + } + }); + }).call(this); diff --git a/wqflask/wqflask/templates/mapping_results.html b/wqflask/wqflask/templates/mapping_results.html index 7e05be18..94ac0350 100644 --- a/wqflask/wqflask/templates/mapping_results.html +++ b/wqflask/wqflask/templates/mapping_results.html @@ -274,7 +274,11 @@ {% endif %} <td align="right">{{marker.chr}}</td> {% if plotScale != "physic" %} + {% if 'cM' in marker %} + <td align="right">{{ '%0.3f' | format(marker.cM|float) }}</td> + {% else %} <td align="right">{{ '%0.3f' | format(marker.Mb|float) }}</td> + {% endif %} {% else %} <td align="right">{{ '%0.6f' | format(marker.Mb|float) }}</td> {% endif %} diff --git a/wqflask/wqflask/templates/show_trait_mapping_tools.html b/wqflask/wqflask/templates/show_trait_mapping_tools.html index 777d4a2d..7c897409 100755 --- a/wqflask/wqflask/templates/show_trait_mapping_tools.html +++ b/wqflask/wqflask/templates/show_trait_mapping_tools.html @@ -121,6 +121,16 @@ </div> {% if genofiles and genofiles|length>0 %} <div class="mapping_method_fields form-group"> + <label for="scale_select" style="text-align: right;" class="col-xs-3 control-label">Map Scale</label> + <div style="margin-left:20px;" class="col-xs-2 controls"> + <select id="scale_reaper" class="form-control" style="width: 80px;"> + {% for item in scales_in_geno[genofiles[0]['location']] %} + <option value="{{ item[0] }}">{{ item[1] }}</option> + {% endfor %} + </select> + </div> + </div> + <div class="mapping_method_fields form-group"> <label style="text-align: right;" for="genofiles" class="col-xs-3 control-label">Genotypes</label> <div style="margin-left:20px;" class="col-xs-4 controls"> <select id="genofile_reaper" class="form-control"> @@ -225,7 +235,17 @@ </select> </div> </div> - {% if genofiles and genofiles|length>0 %} + {% if genofiles and genofiles|length > 0 %} + <div class="mapping_method_fields form-group"> + <label for="scale_select" style="text-align: right;" class="col-xs-3 control-label">Map Scale</label> + <div style="margin-left:20px;" class="col-xs-2 controls"> + <select id="scale_rqtl_geno" class="form-control" style="width: 80px;"> + {% for item in scales_in_geno[genofiles[0]['location']] %} + <option value="{{ item[0] }}">{{ item[1] }}</option> + {% endfor %} + </select> + </div> + </div> <div class="mapping_method_fields form-group"> <label style="text-align:right;" for="genofiles" class="col-xs-3 control-label">Genotypes</label> <div style="margin-left:20px;" class="col-xs-4 controls"> |