From 6c49d8bcf0c17057b441d4985b5813a13c990c2b Mon Sep 17 00:00:00 2001 From: zsloan Date: Thu, 14 Jan 2016 23:50:37 +0000 Subject: Regular and single chr views work for qtl reaper now Single chr now works for R/qtl (didn't before) Merged qtl reaper code into marker_regression.py --- .../wqflask/interval_mapping/interval_mapping.py | 196 +-------------------- .../wqflask/marker_regression/marker_regression.py | 75 +++++++- .../marker_regression/marker_regression_gn1.py | 18 +- .../new/javascript/show_trait_mapping_tools.js | 4 +- .../wqflask/templates/marker_regression_gn1.html | 5 + .../templates/show_trait_mapping_tools.html | 18 ++ wqflask/wqflask/views.py | 44 ++--- 7 files changed, 135 insertions(+), 225 deletions(-) diff --git a/wqflask/wqflask/interval_mapping/interval_mapping.py b/wqflask/wqflask/interval_mapping/interval_mapping.py index 46aac86c..e8ac99d7 100755 --- a/wqflask/wqflask/interval_mapping/interval_mapping.py +++ b/wqflask/wqflask/interval_mapping/interval_mapping.py @@ -58,7 +58,7 @@ class IntervalMapping(object): self.json_data = {} self.json_data['lodnames'] = ['lod.hk'] - self.gen_reaper_results(tempdata) + self.gen_reaper_results() #Get chromosome lengths for drawing the interval map plot chromosome_mb_lengths = {} @@ -102,76 +102,8 @@ class IntervalMapping(object): else: self.control_locus = None - def gen_qtl_results(self, tempdata): - """Generates qtl results for plotting interval map""" - self.dataset.group.get_markers() - genotype = self.dataset.group.read_genotype_file() - - samples, values, variances = self.this_trait.export_informative() - - trimmed_samples = [] - trimmed_values = [] - for i in range(0, len(samples)): - if samples[i] in self.dataset.group.samplelist: - trimmed_samples.append(samples[i]) - trimmed_values.append(values[i]) - - #if self.weighted_regression: - # self.lrs_array = self.genotype.permutation(strains = trimmed_samples, - # trait = trimmed_values, - # variance = variances, - # nperm=self.num_permutations) - #else: - self.lrs_array = genotype.permutation(strains = trimmed_samples, - trait = trimmed_values, - nperm= self.num_permutations) - self.suggestive = self.lrs_array[int(self.num_permutations*0.37-1)] - self.significant = self.lrs_array[int(self.num_permutations*0.95-1)] - - print("samples:", trimmed_samples) - - if self.control_locus: - #if self.weighted_regression: - # self.qtl_results = self.dataset.genotype.regression(strains = samples, - # trait = values, - # variance = variances, - # control = self.control_locus) - #else: - reaper_results = genotype.regression(strains = trimmed_samples, - trait = trimmed_values, - control = self.control_locus) - else: - #if self.weighted_regression: - # self.qtl_results = self.dataset.genotype.regression(strains = samples, - # trait = values, - # variance = variances) - #else: - print("strains:", trimmed_samples) - print("trait:", trimmed_values) - reaper_results = genotype.regression(strains = trimmed_samples, - trait = trimmed_values) - - #Need to convert the QTL objects that qtl reaper returns into a json serializable dictionary - self.qtl_results = [] - self.json_data['chr'] = [] - self.json_data['pos'] = [] - self.json_data['lod.hk'] = [] - self.json_data['markernames'] = [] - for qtl in reaper_results: - reaper_locus = qtl.locus - self.json_data['pos'].append(reaper_locus.Mb) - self.json_data['lod.hk'].append(qtl.lrs) - self.json_data['markernames'].append(reaper_locus.name) - if self.additive: - self.json_data['additive'] = [] - self.json_data['additive'].append(qtl.additive) - locus = {"name":reaper_locus.name, "chr":reaper_locus.chr, "cM":reaper_locus.cM, "Mb":reaper_locus.Mb} - qtl = {"lrs": qtl.lrs, "locus": locus, "additive": qtl.additive} - self.qtl_results.append(qtl) - - - def gen_reaper_results(self, tempdata): + def gen_reaper_results(self): genotype = self.dataset.group.read_genotype_file() samples, values, variances = self.this_trait.export_informative() @@ -249,127 +181,3 @@ class IntervalMapping(object): trimmed_genotype_data.append(new_genotypes) return trimmed_genotype_data - #def get_qtl_results(self): - # """Gets the LOD (or LRS) score at each marker in order do the qtl mapping""" - # - # - # - # #self.genotype = self.genotype.addinterval() - # #resultSlice = [] - # #controlGeno = [] - # - # #if self.multipleInterval: - # # self.suggestive = 0 - # # self.significance = 0 - # # if self.selectedChr > -1: - # # self.genotype.chromosome = [self.genotype[self.selectedChr]] - # #else: - # #single interval mapping - # #try: - # # self.suggestive = float(fd.formdata.getvalue('permSuggestive')) - # # self.significance = float(fd.formdata.getvalue('permSignificance')) - # #except: - # # self.suggestive = None - # # self.significance = None - # - # #NOT DOING MULTIPLE TRAITS YET, BUT WILL NEED TO LATER - # #_strains, _vals, _vars = self.traitList[0].exportInformative(weightedRegression) - # - # #if webqtlUtil.ListNotNull(_vars): - # # pass - # #else: - # # weightedRegression = 0 - # # _strains, _vals, _vars = self.traitList[0].exportInformative() - # - # ##locate genotype of control Locus - # #if self.controlLocus: - # # controlGeno2 = [] - # # _FIND = 0 - # # for _chr in self.genotype: - # # for _locus in _chr: - # # if _locus.name == self.controlLocus: - # # controlGeno2 = _locus.genotype - # # _FIND = 1 - # # break - # # if _FIND: - # # break - # # if controlGeno2: - # # _prgy = list(self.genotype.prgy) - # # for _strain in _strains: - # # _idx = _prgy.index(_strain) - # # controlGeno.append(controlGeno2[_idx]) - # # else: - # # return "The control marker you selected is not in the genofile." - # # - # # - # # if self.significance and self.suggestive: - # # pass - # # else: - # # if self.permChecked: - # # if weightedRegression: - # # self.LRSArray = self.genotype.permutation(strains = _strains, trait = _vals, - # # variance = _vars, nperm=fd.nperm) - # # else: - # # self.LRSArray = self.genotype.permutation(strains = _strains, trait = _vals, - # # nperm=fd.nperm) - # # self.suggestive = self.LRSArray[int(fd.nperm*0.37-1)] - # # self.significance = self.LRSArray[int(fd.nperm*0.95-1)] - # # - # # else: - # # self.suggestive = 9.2 - # # self.significance = 16.1 - # # - # # #calculating bootstrap - # # #from now on, genotype could only contain a single chromosome - # # #permutation need to be performed genome wide, this is not the case for bootstrap - # # - # # #due to the design of qtlreaper, composite regression need to be performed genome wide - # # if not self.controlLocus and self.selectedChr > -1: - # # self.genotype.chromosome = [self.genotype[self.selectedChr]] - # # elif self.selectedChr > -1: #self.controlLocus and self.selectedChr > -1 - # # lociPerChr = map(len, self.genotype) - # # resultSlice = reduce(lambda X, Y: X+Y, lociPerChr[:self.selectedChr], 0) - # # resultSlice = [resultSlice,resultSlice+lociPerChr[self.selectedChr]] - # # else: - # # pass - # - # #calculate QTL for each trait - # self.qtl_results = [] - # - # #for thisTrait in self.traitList: - # _strains, _vals, _vars = thisTrait.exportInformative(weightedRegression) - # if self.controlLocus: - # if weightedRegression: - # qtlresult = self.genotype.regression(strains = _strains, trait = _vals, - # variance = _vars, control = self.controlLocus) - # else: - # qtlresult = self.genotype.regression(strains = _strains, trait = _vals, - # control = self.controlLocus) - # if resultSlice: - # qtlresult = qtlresult[resultSlice[0]:resultSlice[1]] - # else: - # if weightedRegression: - # qtlresult = self.genotype.regression(strains = _strains, trait = _vals, - # variance = _vars) - # else: - # qtlresult = self.genotype.regression(strains = _strains, trait = _vals) - # - # self.qtlresults.append(qtlresult) - # - # if not self.multipleInterval: - # if self.controlLocus and self.selectedChr > -1: - # self.genotype.chromosome = [self.genotype[self.selectedChr]] - # - # if self.bootChecked: - # if controlGeno: - # self.bootResult = self.genotype.bootstrap(strains = _strains, trait = _vals, - # control = controlGeno, nboot=fd.nboot) - # elif weightedRegression: - # self.bootResult = self.genotype.bootstrap(strains = _strains, trait = _vals, - # variance = _vars, nboot=fd.nboot) - # else: - # self.bootResult = self.genotype.bootstrap(strains = _strains, trait = _vals, - # nboot=fd.nboot) - # else: - # self.bootResult = [] - diff --git a/wqflask/wqflask/marker_regression/marker_regression.py b/wqflask/wqflask/marker_regression/marker_regression.py index c396b14b..605fde63 100755 --- a/wqflask/wqflask/marker_regression/marker_regression.py +++ b/wqflask/wqflask/marker_regression/marker_regression.py @@ -110,6 +110,15 @@ class MarkerRegression(object): results = self.run_rqtl_geno() #print("qtl_results:", results) + elif self.mapping_method == "reaper": + if start_vars['num_perm'] == "": + self.num_perm = 0 + else: + self.num_perm = int(start_vars['num_perm']) + self.additive = False + self.control = start_vars['control_marker'] + self.do_control = start_vars['do_control'] + results = self.gen_reaper_results() elif self.mapping_method == "plink": results = self.run_plink() #print("qtl_results:", pf(results)) @@ -130,7 +139,7 @@ class MarkerRegression(object): if marker['chr1'] > 0 or marker['chr1'] == "X" or marker['chr1'] == "X/Y": if marker['chr1'] > highest_chr or marker['chr1'] == "X" or marker['chr1'] == "X/Y": highest_chr = marker['chr1'] - if 'lod_score' in marker: + if 'lod_score' in marker.keys(): self.qtl_results.append(marker) for qtl in enumerate(self.qtl_results): @@ -157,7 +166,7 @@ class MarkerRegression(object): if marker['chr'] > 0 or marker['chr'] == "X" or marker['chr'] == "X/Y": if marker['chr'] > highest_chr or marker['chr'] == "X" or marker['chr'] == "X/Y": highest_chr = marker['chr'] - if 'lod_score' in marker: + if ('lod_score' in marker.keys()) or ('lrs_value' in marker.keys()): self.qtl_results.append(marker) self.json_data['chr'] = [] @@ -178,7 +187,7 @@ class MarkerRegression(object): else: self.json_data['chr'].append(str(qtl['chr'])) self.json_data['pos'].append(qtl['Mb']) - if 'lrs_value' in qtl: + if 'lrs_value' in qtl.keys(): self.json_data['lod.hk'].append(str(qtl['lrs_value'])) else: self.json_data['lod.hk'].append(str(qtl['lod_score'])) @@ -588,6 +597,66 @@ class MarkerRegression(object): return sample_list + def gen_reaper_results(self): + genotype = self.dataset.group.read_genotype_file() + + samples, values, variances = self.this_trait.export_informative() + + trimmed_samples = [] + trimmed_values = [] + for i in range(0, len(samples)): + if samples[i] in self.dataset.group.samplelist: + trimmed_samples.append(samples[i]) + trimmed_values.append(values[i]) + + self.lrs_array = genotype.permutation(strains = trimmed_samples, + trait = trimmed_values, + nperm= self.num_perm) + + self.suggestive = self.lrs_array[int(self.num_perm*0.37-1)] + self.significant = self.lrs_array[int(self.num_perm*0.95-1)] + self.json_data['suggestive'] = self.suggestive + self.json_data['significant'] = self.significant + + print("samples:", trimmed_samples) + + if self.control != "" and self.do_control == "true": + reaper_results = genotype.regression(strains = trimmed_samples, + trait = trimmed_values, + control = self.control) + else: + reaper_results = genotype.regression(strains = trimmed_samples, + trait = trimmed_values) + + self.json_data['chr'] = [] + self.json_data['pos'] = [] + self.json_data['lod.hk'] = [] + self.json_data['markernames'] = [] + if self.additive: + self.json_data['additive'] = [] + + #Need to convert the QTL objects that qtl reaper returns into a json serializable dictionary + qtl_results = [] + for qtl in reaper_results: + reaper_locus = qtl.locus + #ZS: Convert chr to int + converted_chr = reaper_locus.chr + if reaper_locus.chr != "X" and reaper_locus.chr != "X/Y": + converted_chr = int(reaper_locus.chr) + self.json_data['chr'].append(converted_chr) + self.json_data['pos'].append(reaper_locus.Mb) + self.json_data['lod.hk'].append(qtl.lrs) + self.json_data['markernames'].append(reaper_locus.name) + if self.additive: + self.json_data['additive'].append(qtl.additive) + locus = {"name":reaper_locus.name, "chr":reaper_locus.chr, "cM":reaper_locus.cM, "Mb":reaper_locus.Mb} + qtl = {"lrs_value": qtl.lrs, "chr":converted_chr, "Mb":reaper_locus.Mb, + "cM":reaper_locus.cM, "name":reaper_locus.name, "additive":qtl.additive, "dominance":qtl.dominance} + qtl_results.append(qtl) + + return qtl_results + + def parse_plink_output(self, output_filename): plink_results={} diff --git a/wqflask/wqflask/marker_regression/marker_regression_gn1.py b/wqflask/wqflask/marker_regression/marker_regression_gn1.py index 23d5a631..a7b5fad3 100644 --- a/wqflask/wqflask/marker_regression/marker_regression_gn1.py +++ b/wqflask/wqflask/marker_regression/marker_regression_gn1.py @@ -178,9 +178,13 @@ class MarkerRegression(object): self.this_trait = start_vars['this_trait'] self.species = start_vars['species'] - self.vals = start_vars['vals'] #Needed to put into form for when it is resubmitted for single chr views or remapping - - self.mapping_method = start_vars['mapping_method'] #Needing for form submission when doing single chr mapping or remapping after changing options + #Needing for form submission when doing single chr mapping or remapping after changing options + self.vals = start_vars['vals'] + self.mapping_method = start_vars['mapping_method'] + if self.mapping_method == "rqtl_geno": + self.mapmethod_rqtl_geno = start_vars['method'] + self.mapmodel_rqtl_geno = start_vars['model'] + self.pair_scan = start_vars['pair_scan'] self.js_data = start_vars['js_data'] @@ -225,10 +229,15 @@ class MarkerRegression(object): #self.permChecked = fd.formdata.getvalue('permCheck', True) self.bootChecked = False #ZS: For now setting to False, I'll add this option later once rest of figure works #self.bootChecked = fd.formdata.getvalue('bootCheck', '') + if 'do_control' in start_vars.keys(): + self.doControl = start_vars['do_control'] + else: + self.doControl = "false" if 'control' in start_vars.keys(): self.controlLocus = start_vars['control'] else: self.controlLocus = "" + #self.controlLocus = fd.formdata.getvalue('controlLocus', '') #try: @@ -1760,6 +1769,7 @@ class MarkerRegression(object): lineColor = pid.lightblue startPosX = xLeftOffset for j, ChrInfo in enumerate(ChrAInfo): + if ChrInfo == self.selectedChr: preLpos = -1 for i, item in enumerate(ChrInfo): Lname,Lpos = item @@ -1932,8 +1942,6 @@ class MarkerRegression(object): previous_chr = qtlresult['chr'] previous_chr_as_int += 1 - print("ChrLengthDistList:", self.ChrLengthDistList) - newStartPosX = (self.ChrLengthDistList[previous_chr_as_int - 1]+self.GraphInterval)*plotXScale if newStartPosX != oldStartPosX: startPosX += newStartPosX 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 16eeb8df..a48ac2b7 100755 --- a/wqflask/wqflask/static/new/javascript/show_trait_mapping_tools.js +++ b/wqflask/wqflask/static/new/javascript/show_trait_mapping_tools.js @@ -236,9 +236,11 @@ var form_data, url; console.log("In interval mapping"); //$("#progress_bar_container").modal(); - url = "/interval_mapping"; + url = "/marker_regression"; $('input[name=method]').val("reaper"); $('input[name=manhattan_plot]').val($('input[name=manhattan_plot_reaper]:checked').val()); + $('input[name=control_marker]').val($('input[name=control_reaper]').val()); + $('input[name=do_control]').val($('input[name=do_control_reaper]:checked').val()); $('input[name=mapping_display_all]').val($('input[name=display_all_reaper]')); $('input[name=suggestive]').val($('input[name=suggestive_reaper]')); form_data = $('#trait_data_form').serialize(); diff --git a/wqflask/wqflask/templates/marker_regression_gn1.html b/wqflask/wqflask/templates/marker_regression_gn1.html index 05623462..9f7dd41f 100644 --- a/wqflask/wqflask/templates/marker_regression_gn1.html +++ b/wqflask/wqflask/templates/marker_regression_gn1.html @@ -22,6 +22,11 @@ + + + + +
{{ gifmap|safe }} diff --git a/wqflask/wqflask/templates/show_trait_mapping_tools.html b/wqflask/wqflask/templates/show_trait_mapping_tools.html index de25f4ee..5b0c74c3 100755 --- a/wqflask/wqflask/templates/show_trait_mapping_tools.html +++ b/wqflask/wqflask/templates/show_trait_mapping_tools.html @@ -187,6 +187,24 @@
+
+ +
+ {% if dataset.type == 'ProbeSet' and this_trait.locus_chr != "" %} + + {% else %} + + {% endif %} + + +
+