From 897f3291f9d943b4deaef7bc924197c268bc8b3f Mon Sep 17 00:00:00 2001 From: zsloan Date: Fri, 3 Apr 2020 16:57:59 -0500 Subject: Added option to select cofactors from collection for R/qtl Fixed some issues with displaying the R/qtl results from the RData HET3-ITP dataa --- wqflask/utility/gen_geno_ob.py | 4 +- .../marker_regression/display_mapping_results.py | 36 ++++--- wqflask/wqflask/marker_regression/rqtl_mapping.py | 116 ++++++++++++++++----- wqflask/wqflask/marker_regression/run_mapping.py | 4 +- .../wqflask/static/new/javascript/show_trait.js | 7 ++ .../templates/show_trait_mapping_tools.html | 18 +++- 6 files changed, 140 insertions(+), 45 deletions(-) diff --git a/wqflask/utility/gen_geno_ob.py b/wqflask/utility/gen_geno_ob.py index aa5b27c4..db40f6ea 100644 --- a/wqflask/utility/gen_geno_ob.py +++ b/wqflask/utility/gen_geno_ob.py @@ -45,10 +45,12 @@ class genotype(object): chr_ob = None for marker in qtl_results: locus = Locus(self) - if str(marker['chr']) != this_chr: + if (str(marker['chr']) != this_chr) and this_chr != "X": #ZS: This is really awkward but works as a temporary fix if this_chr != "": self.chromosomes.append(chr_ob) this_chr = str(marker['chr']) + if this_chr == "20": + this_chr = "X" chr_ob = Chr(this_chr, self) if 'chr' in marker: locus.chr = str(marker['chr']) diff --git a/wqflask/wqflask/marker_regression/display_mapping_results.py b/wqflask/wqflask/marker_regression/display_mapping_results.py index 302a4ff9..b8f84721 100644 --- a/wqflask/wqflask/marker_regression/display_mapping_results.py +++ b/wqflask/wqflask/marker_regression/display_mapping_results.py @@ -334,6 +334,10 @@ class DisplayMappingResults(object): ################################################################ self.ChrList = [("All", -1)] for i, indChr in enumerate(self.genotype): + if self.dataset.group.species == "mouse" and indChr.name == "20": + self.ChrList.append(("X", i)) + elif self.dataset.group.species == "rat" and indChr.name == "21": + self.ChrList.append(("X", i)) self.ChrList.append((indChr.name, i)) self.ChrLengthMbList = g.db.execute(""" @@ -356,10 +360,7 @@ class DisplayMappingResults(object): self.ChrLengthCMList = [] for i, _chr in enumerate(self.genotype): - if self.mapping_method == "rqtl_geno" and self.genotype.filler == True: - self.ChrLengthCMList.append(_chr[-1].cM) - else: - self.ChrLengthCMList.append(_chr[-1].cM - _chr[0].cM) + self.ChrLengthCMList.append(_chr[-1].cM - _chr[0].cM) self.ChrLengthCMSum = reduce(lambda x, y:x+y, self.ChrLengthCMList, 0.0) @@ -1568,10 +1569,7 @@ class DisplayMappingResults(object): if self.selectedChr == -1: #ZS: If viewing full genome/all chromosomes for i, _chr in enumerate(self.genotype): thisChr = [] - if self.mapping_method == "rqtl_geno" and self.genotype.filler == True: - Locus0CM = 0 - else: - Locus0CM = _chr[0].cM + Locus0CM = _chr[0].cM nLoci = len(_chr) if nLoci <= 8: for _locus in _chr: @@ -1593,10 +1591,7 @@ class DisplayMappingResults(object): for i, _chr in enumerate(self.genotype): if _chr.name == self.ChrList[self.selectedChr][0]: thisChr = [] - if self.mapping_method == "rqtl_geno" and self.genotype.filler == True: - Locus0CM = 0 - else: - Locus0CM = _chr[0].cM + Locus0CM = _chr[0].cM for _locus in _chr: if _locus.name != ' - ': if _locus.cM != preLpos: @@ -1892,12 +1887,21 @@ class DisplayMappingResults(object): oldStartPosX = newStartPosX #ZS: This is beause the chromosome value stored in qtlresult['chr'] can be (for example) either X or 20 depending upon the mapping method/scale used - if self.plotScale == "physic": - this_chr = str(self.ChrList[self.selectedChr][0]) - else: + this_chr = str(self.ChrList[self.selectedChr][0]) + if self.plotScale != "physic": this_chr = str(self.ChrList[self.selectedChr][1]+1) + if self.selectedChr == -1 or str(qtlresult['chr']) == this_chr: - Xc = startPosX + (qtlresult['Mb']-startMb)*plotXScale + 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)) + else: + Xc = startPosX + (qtlresult['Mb']-startMb)*plotXScale + # updated by NL 06-18-2011: # fix the over limit LRS graph issue since genotype trait may give infinite LRS; # for any lrs is over than 460(LRS max in this system), it will be reset to 460 diff --git a/wqflask/wqflask/marker_regression/rqtl_mapping.py b/wqflask/wqflask/marker_regression/rqtl_mapping.py index 2e3ea406..d76f3812 100644 --- a/wqflask/wqflask/marker_regression/rqtl_mapping.py +++ b/wqflask/wqflask/marker_regression/rqtl_mapping.py @@ -1,14 +1,17 @@ import rpy2.robjects as ro +import rpy2.robjects.numpy2ri as np2r import numpy as np from base.webqtlConfig import TMPDIR +from base.trait import GeneralTrait +from base.data_set import create_dataset from utility import webqtlUtil from utility.tools import locate, TEMPDIR import utility.logger logger = utility.logger.getLogger(__name__ ) -def run_rqtl_geno(vals, dataset, method, model, permCheck, num_perm, do_control, control_marker, manhattan_plot, pair_scan): +def run_rqtl_geno(vals, dataset, method, model, permCheck, num_perm, do_control, control_marker, manhattan_plot, pair_scan, samples, 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 @@ -43,14 +46,22 @@ def run_rqtl_geno(vals, dataset, method, model, permCheck, num_perm, do_control, else: cross_object = calc_genoprob(cross_object, step=1, stepwidth="max") - cross_object = add_phenotype(cross_object, sanitize_rqtl_phenotype(vals)) # Add the phenotype + cross_object = add_phenotype(cross_object, sanitize_rqtl_phenotype(vals), "the_pheno") # Add the phenotype # Scan for QTLs - covar = create_covariates(control_marker, cross_object) # Create the additive covariate matrix + marker_covars = create_marker_covariates(control_marker, cross_object) # Create the additive covariate markers + + if cofactors != "": + cross_object, trait_covars = add_cofactors(cross_object, dataset, cofactors, samples) # Create the covariates from selected traits + ro.r('all_covars <- cbind(marker_covars, trait_covars)') + else: + ro.r('all_covars <- marker_covars') + + covars = ro.r['all_covars'] if pair_scan: if do_control == "true": - logger.info("Using covariate"); result_data_frame = scantwo(cross_object, pheno = "the_pheno", addcovar = covar, model=model, method=method, n_cluster = 16) + logger.info("Using covariate"); result_data_frame = scantwo(cross_object, pheno = "the_pheno", addcovar = covars, model=model, method=method, n_cluster = 16) else: logger.info("No covariates"); result_data_frame = scantwo(cross_object, pheno = "the_pheno", model=model, method=method, n_cluster = 16) @@ -61,14 +72,14 @@ def run_rqtl_geno(vals, dataset, method, model, permCheck, num_perm, do_control, return process_pair_scan_results(result_data_frame) else: - if do_control == "true": - logger.info("Using covariate"); result_data_frame = scanone(cross_object, pheno = "the_pheno", addcovar = covar, model=model, method=method) + if do_control == "true" or cofactors != "": + logger.info("Using covariate"); result_data_frame = scanone(cross_object, pheno = "the_pheno", addcovar = covars, model=model, method=method) else: logger.info("No covariates"); result_data_frame = scanone(cross_object, pheno = "the_pheno", model=model, method=method) if num_perm > 0 and permCheck == "ON": # Do permutation (if requested by user) - if do_control == "true": - perm_data_frame = scanone(cross_object, pheno_col = "the_pheno", addcovar = covar, n_perm = num_perm, model=model, method=method) + if do_control == "true" or cofactors != "": + perm_data_frame = scanone(cross_object, pheno_col = "the_pheno", addcovar = covars, n_perm = num_perm, model=model, method=method) else: perm_data_frame = scanone(cross_object, pheno_col = "the_pheno", n_perm = num_perm, model=model, method=method) @@ -118,23 +129,6 @@ def generate_cross_from_geno(dataset): # TODO: Need to figure out why som } """ % (dataset.group.name + ".geno")) -def add_phenotype(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(control_marker, cross): - ro.globalenv["the_cross"] = cross - ro.r('genotypes <- pull.geno(the_cross)') # Get the genotype matrix - userinputS = control_marker.replace(" ", "").split(",") # TODO: sanitize user input, Never Ever trust a user - covariate_names = ', '.join('"{0}"'.format(w) for w in userinputS) - 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 - return ro.r["covariates"] - def sanitize_rqtl_phenotype(vals): pheno_as_string = "c(" for i, val in enumerate(vals): @@ -151,6 +145,78 @@ def sanitize_rqtl_phenotype(vals): pheno_as_string += ")" return pheno_as_string +def add_phenotype(cross, pheno_as_string, col_name): + ro.globalenv["the_cross"] = cross + ro.r('the_cross$pheno <- cbind(pull.pheno(the_cross), ' + col_name + ' = '+ pheno_as_string +')') + return ro.r["the_cross"] + +def pull_covar(cross, covar_name_string): + ro.globalenv["the_cross"] = cross + ro.r('trait_covars <- pull.pheno(the_cross, ' + covar_name_string + ')') + + return ro.r["trait_covars"] + +def add_cofactors(cross, this_dataset, covariates, samples): + ro.numpy2ri.activate() + + covariate_list = covariates.split(",") + covar_name_string = "c(" + for i, covariate in enumerate(covariate_list): + this_covar_data = [] + covar_as_string = "c(" + trait_name = covariate.split(":")[0] + dataset_ob = create_dataset(covariate.split(":")[1]) + trait_ob = GeneralTrait(dataset=dataset_ob, + name=trait_name, + cellid=None) + + this_dataset.group.get_samplelist() + trait_samples = this_dataset.group.samplelist + trait_sample_data = trait_ob.data + for index, sample in enumerate(trait_samples): + if sample in samples: + if sample in trait_sample_data: + sample_value = trait_sample_data[sample].value + this_covar_data.append(sample_value) + else: + this_covar_data.append("NA") + + for j, item in enumerate(this_covar_data): + if j < (len(this_covar_data) - 1): + covar_as_string += str(item) + "," + else: + covar_as_string += str(item) + + covar_as_string += ")" + + col_name = "covar_" + str(i) + + if i < (len(covariate_list) - 1): + covar_name_string += '"' + col_name + '", ' + else: + covar_name_string += '"' + col_name + '"' + + cross = add_phenotype(cross, covar_as_string, col_name) + + covar_name_string += ")" + + covars = pull_covar(cross, covar_name_string) + + return cross, covars + +def create_marker_covariates(control_marker, cross): + ro.globalenv["the_cross"] = cross + ro.r('genotypes <- pull.geno(the_cross)') # Get the genotype matrix + userinputS = control_marker.replace(" ", "").split(",") # TODO: sanitize user input, Never Ever trust a user + covariate_names = ', '.join('"{0}"'.format(w) for w in userinputS) + 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('marker_covars <- genotypes[,covnames]') # Get the covariate matrix by using the marker name as index to the genotype file + + return ro.r["marker_covars"] + def process_pair_scan_results(result): pair_scan_results = [] diff --git a/wqflask/wqflask/marker_regression/run_mapping.py b/wqflask/wqflask/marker_regression/run_mapping.py index f03b046e..3006c4ff 100644 --- a/wqflask/wqflask/marker_regression/run_mapping.py +++ b/wqflask/wqflask/marker_regression/run_mapping.py @@ -216,9 +216,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 = rqtl_mapping.run_rqtl_geno(self.vals, self.dataset, self.method, self.model, self.permCheck, self.num_perm, self.do_control, self.control_marker, self.manhattan_plot, self.pair_scan) + self.perm_output, self.suggestive, self.significant, results = rqtl_mapping.run_rqtl_geno(self.vals, self.dataset, self.method, self.model, self.permCheck, self.num_perm, self.do_control, self.control_marker, self.manhattan_plot, self.pair_scan, self.samples, self.covariates) else: - results = rqtl_mapping.run_rqtl_geno(self.vals, self.dataset, self.method, self.model, self.permCheck, self.num_perm, self.do_control, self.control_marker, self.manhattan_plot, self.pair_scan) + results = rqtl_mapping.run_rqtl_geno(self.vals, self.dataset, self.method, self.model, self.permCheck, self.num_perm, self.do_control, self.control_marker, self.manhattan_plot, self.pair_scan, self.samples, 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: diff --git a/wqflask/wqflask/static/new/javascript/show_trait.js b/wqflask/wqflask/static/new/javascript/show_trait.js index 0162f858..a2fa37e0 100644 --- a/wqflask/wqflask/static/new/javascript/show_trait.js +++ b/wqflask/wqflask/static/new/javascript/show_trait.js @@ -107,6 +107,13 @@ $("#remove_covariates").click(function () { $("input[name=covariates]").val("") $(".selected_covariates").val("") }); +$(".select_covariates").click(function () { + open_covariate_selection(); +}); +$(".remove_covariates").click(function () { + $("input[name=covariates]").val("") + $(".selected_covariates").val("") +}); d3.select("#clear_compare_trait").on("click", (function(_this) { return function() { return $('.scatter-matrix-container').remove(); diff --git a/wqflask/wqflask/templates/show_trait_mapping_tools.html b/wqflask/wqflask/templates/show_trait_mapping_tools.html index 6da25a5d..66903c31 100644 --- a/wqflask/wqflask/templates/show_trait_mapping_tools.html +++ b/wqflask/wqflask/templates/show_trait_mapping_tools.html @@ -325,7 +325,23 @@ - +
+
+ + Select covariate(s) from a collection +
+
+ {% if g.user_session.num_collections < 1 %} + No collections available. Please add traits to a collection to use them as covariates. + {% else %} +
+ + +
+ + {% endif %} +
+
-- cgit v1.2.3