diff options
Diffstat (limited to 'wqflask')
| -rw-r--r-- | wqflask/wqflask/interval_mapping/interval_mapping.py | 370 |
1 files changed, 161 insertions, 209 deletions
diff --git a/wqflask/wqflask/interval_mapping/interval_mapping.py b/wqflask/wqflask/interval_mapping/interval_mapping.py index 3ec08f6b..48e8018e 100644 --- a/wqflask/wqflask/interval_mapping/interval_mapping.py +++ b/wqflask/wqflask/interval_mapping/interval_mapping.py @@ -52,7 +52,7 @@ class IntervalMapping(object): self.set_options(start_vars) - self.gen_data(tempdata) + self.gen_qtl_results(tempdata) #Get chromosome lengths for drawing the interval map plot chromosome_mb_lengths = {} @@ -60,217 +60,11 @@ class IntervalMapping(object): chromosome_mb_lengths[key] = self.species.chromosomes.chromosomes[key].mb_length self.js_data = dict( + lrs_lod = self.lrs_lod, chromosomes = chromosome_mb_lengths, qtl_results = self.qtl_results, ) - def get_qtl_results(self): - """Gets the LOD (or LRS) score at each marker in order do the qtl mapping""" - - - #INTERCROSS = (self.genotype.type=="intercross") - - #THESE ARE OLD COMMENTS - #draw the LRS scale - #We first determine whether or not we are using a sliding scale. - #If so, we need to compute the maximum LRS value to determine where the max y-value should be, and call this LRSMax. - #LRSTop is then defined to be above the LRSMax by enough to add one additional LRSScale increment. - #if we are using a set-scale, then we set LRSTop to be the user's value, and LRSMax doesn't matter. - - #if self.lrs_lod == 'LOD': - # lodm = self.LODFACTOR - #else: - # lodm = 1.0 - - if self.lrsMax <= 0: #sliding scale - LRSMax = max(map(max, self.qtlresults)).lrs - #genotype trait will give infinite LRS - LRSMax = min(LRSMax, webqtlConfig.MAXLRS) - if self.permChecked and not self.multipleInterval: - self.significance = min(self.significance, webqtlConfig.MAXLRS) - self.suggestive = min(self.suggestive, webqtlConfig.MAXLRS) - LRSMax = max(self.significance, LRSMax) - else: - LRSMax = self.lrsMax*lodm - - if LRSMax/lodm > 100: - LRSScale = 20.0 - elif LRSMax/lodm > 20: - LRSScale = 5.0 - elif LRSMax/lodm > 7.5: - LRSScale = 2.5 - else: - LRSScale = 1.0 - - LRSAxisList = Plot.frange(LRSScale, LRSMax/lodm, LRSScale) - #make sure the user's value appears on the y-axis - #update by NL 6-21-2011: round the LOD value to 100 when LRSMax is equal to 460 - LRSAxisList.append(round(LRSMax/lodm)) - - #draw the "LRS" or "LOD" string to the left of the axis - LRSScaleFont=pid.Font(ttf="verdana", size=14*fontZoom, bold=0) - LRSLODFont=pid.Font(ttf="verdana", size=14*zoom*1.5, bold=0) - yZero = yTopOffset + plotHeight - - - canvas.drawString(self.LRS_LOD, xLeftOffset - canvas.stringWidth("999.99", font=LRSScaleFont) - 10*zoom, \ - yZero - 150, font=LRSLODFont, color=pid.black, angle=90) - - for item in LRSAxisList: - yLRS = yZero - (item*lodm/LRSMax) * LRSHeightThresh - canvas.drawLine(xLeftOffset, yLRS, xLeftOffset - 4, yLRS, color=self.LRS_COLOR, width=1*zoom) - scaleStr = "%2.1f" % item - canvas.drawString(scaleStr, xLeftOffset-4-canvas.stringWidth(scaleStr, font=LRSScaleFont)-5, yLRS+3, font=LRSScaleFont, color=self.LRS_COLOR) - - - #"Significant" and "Suggestive" Drawing Routine - # ======= Draw the thick lines for "Significant" and "Suggestive" ===== (crowell: I tried to make the SNPs draw over these lines, but piddle wouldn't have it...) - if self.permChecked and not self.multipleInterval: - significantY = yZero - self.significance*LRSHeightThresh/LRSMax - suggestiveY = yZero - self.suggestive*LRSHeightThresh/LRSMax - startPosX = xLeftOffset - for i, _chr in enumerate(self.genotype): - rightEdge = int(startPosX + self.ChrLengthDistList[i]*plotXScale - self.SUGGESTIVE_WIDTH/1.5) - canvas.drawLine(startPosX+self.SUGGESTIVE_WIDTH/1.5, suggestiveY, rightEdge, suggestiveY, color=self.SUGGESTIVE_COLOR, - width=self.SUGGESTIVE_WIDTH*zoom, clipX=(xLeftOffset, xLeftOffset + plotWidth-2)) - canvas.drawLine(startPosX+self.SUGGESTIVE_WIDTH/1.5, significantY, rightEdge, significantY, color=self.SIGNIFICANT_COLOR, - width=self.SIGNIFICANT_WIDTH*zoom, clipX=(xLeftOffset, xLeftOffset + plotWidth-2)) - sugg_coords = "%d, %d, %d, %d" % (startPosX, suggestiveY-2, rightEdge + 2*zoom, suggestiveY+2) - sig_coords = "%d, %d, %d, %d" % (startPosX, significantY-2, rightEdge + 2*zoom, significantY+2) - if self.LRS_LOD == 'LRS': - sugg_title = "Suggestive LRS = %0.2f" % self.suggestive - sig_title = "Significant LRS = %0.2f" % self.significance - else: - sugg_title = "Suggestive LOD = %0.2f" % (self.suggestive/4.61) - sig_title = "Significant LOD = %0.2f" % (self.significance/4.61) - Areas1 = HT.Area(shape='rect',coords=sugg_coords,title=sugg_title) - Areas2 = HT.Area(shape='rect',coords=sig_coords,title=sig_title) - gifmap.areas.append(Areas1) - gifmap.areas.append(Areas2) - - startPosX += (self.ChrLengthDistList[i]+self.GraphInterval)*plotXScale - - - if self.multipleInterval: - lrsEdgeWidth = 1 - else: - additiveMax = max(map(lambda X : abs(X.additive), self.qtlresults[0])) - if INTERCROSS: - dominanceMax = max(map(lambda X : abs(X.dominance), self.qtlresults[0])) - else: - dominanceMax = -1 - lrsEdgeWidth = 2 - for i, qtlresult in enumerate(self.qtlresults): - m = 0 - startPosX = xLeftOffset - thisLRSColor = self.colorCollection[i] - for j, _chr in enumerate(self.genotype): - LRSCoordXY = [] - AdditiveCoordXY = [] - DominanceCoordXY = [] - for k, _locus in enumerate(_chr): - if self.plotScale == 'physic': - Xc = startPosX + (_locus.Mb-startMb)*plotXScale - else: - Xc = startPosX + (_locus.cM-_chr[0].cM)*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 - if qtlresult[m].lrs> 460 or qtlresult[m].lrs=='inf': - Yc = yZero - webqtlConfig.MAXLRS*LRSHeightThresh/LRSMax - else: - Yc = yZero - qtlresult[m].lrs*LRSHeightThresh/LRSMax - - LRSCoordXY.append((Xc, Yc)) - if not self.multipleInterval and self.additiveChecked: - Yc = yZero - qtlresult[m].additive*AdditiveHeightThresh/additiveMax - AdditiveCoordXY.append((Xc, Yc)) - if not self.multipleInterval and INTERCROSS and self.additiveChecked: - Yc = yZero - qtlresult[m].dominance*DominanceHeightThresh/dominanceMax - DominanceCoordXY.append((Xc, Yc)) - m += 1 - canvas.drawPolygon(LRSCoordXY,edgeColor=thisLRSColor,closed=0, edgeWidth=lrsEdgeWidth, clipX=(xLeftOffset, xLeftOffset + plotWidth)) - if not self.multipleInterval and self.additiveChecked: - plusColor = self.ADDITIVE_COLOR_POSITIVE - minusColor = self.ADDITIVE_COLOR_NEGATIVE - for k, aPoint in enumerate(AdditiveCoordXY): - if k > 0: - Xc0, Yc0 = AdditiveCoordXY[k-1] - Xc, Yc = aPoint - if (Yc0-yZero)*(Yc-yZero) < 0: - if Xc == Xc0: #genotype , locus distance is 0 - Xcm = Xc - else: - Xcm = (yZero-Yc0)/((Yc-Yc0)/(Xc-Xc0)) +Xc0 - if Yc0 < yZero: - canvas.drawLine(Xc0, Yc0, Xcm, yZero, color=plusColor, width=1, clipX=(xLeftOffset, xLeftOffset + plotWidth)) - canvas.drawLine(Xcm, yZero, Xc, yZero-(Yc-yZero), color=minusColor, width=1, clipX=(xLeftOffset, xLeftOffset + plotWidth)) - else: - canvas.drawLine(Xc0, yZero - (Yc0-yZero), Xcm, yZero, color=minusColor, width=1, clipX=(xLeftOffset, xLeftOffset + plotWidth)) - canvas.drawLine(Xcm, yZero, Xc, Yc, color=plusColor, width=1, clipX=(xLeftOffset, xLeftOffset + plotWidth)) - elif (Yc0-yZero)*(Yc-yZero) > 0: - if Yc < yZero: - canvas.drawLine(Xc0, Yc0, Xc, Yc, color=plusColor, width=1, clipX=(xLeftOffset, xLeftOffset + plotWidth)) - else: - canvas.drawLine(Xc0, yZero - (Yc0-yZero), Xc, yZero - (Yc-yZero), color=minusColor, width=1, clipX=(xLeftOffset, xLeftOffset + plotWidth)) - else: - minYc = min(Yc-yZero, Yc0-yZero) - if minYc < 0: - canvas.drawLine(Xc0, Yc0, Xc, Yc, color=plusColor, width=1, clipX=(xLeftOffset, xLeftOffset + plotWidth)) - else: - canvas.drawLine(Xc0, yZero - (Yc0-yZero), Xc, yZero - (Yc-yZero), color=minusColor, width=1, clipX=(xLeftOffset, xLeftOffset + plotWidth)) - if not self.multipleInterval and INTERCROSS and self.dominanceChecked: - plusColor = self.DOMINANCE_COLOR_POSITIVE - minusColor = self.DOMINANCE_COLOR_NEGATIVE - for k, aPoint in enumerate(DominanceCoordXY): - if k > 0: - Xc0, Yc0 = DominanceCoordXY[k-1] - Xc, Yc = aPoint - if (Yc0-yZero)*(Yc-yZero) < 0: - if Xc == Xc0: #genotype , locus distance is 0 - Xcm = Xc - else: - Xcm = (yZero-Yc0)/((Yc-Yc0)/(Xc-Xc0)) +Xc0 - if Yc0 < yZero: - canvas.drawLine(Xc0, Yc0, Xcm, yZero, color=plusColor, width=1, clipX=(xLeftOffset, xLeftOffset + plotWidth)) - canvas.drawLine(Xcm, yZero, Xc, yZero-(Yc-yZero), color=minusColor, width=1, clipX=(xLeftOffset, xLeftOffset + plotWidth)) - else: - canvas.drawLine(Xc0, yZero - (Yc0-yZero), Xcm, yZero, color=minusColor, width=1, clipX=(xLeftOffset, xLeftOffset + plotWidth)) - canvas.drawLine(Xcm, yZero, Xc, Yc, color=plusColor, width=1, clipX=(xLeftOffset, xLeftOffset + plotWidth)) - elif (Yc0-yZero)*(Yc-yZero) > 0: - if Yc < yZero: - canvas.drawLine(Xc0, Yc0, Xc, Yc, color=plusColor, width=1, clipX=(xLeftOffset, xLeftOffset + plotWidth)) - else: - canvas.drawLine(Xc0, yZero - (Yc0-yZero), Xc, yZero - (Yc-yZero), color=minusColor, width=1, clipX=(xLeftOffset, xLeftOffset + plotWidth)) - else: - minYc = min(Yc-yZero, Yc0-yZero) - if minYc < 0: - canvas.drawLine(Xc0, Yc0, Xc, Yc, color=plusColor, width=1, clipX=(xLeftOffset, xLeftOffset + plotWidth)) - else: - canvas.drawLine(Xc0, yZero - (Yc0-yZero), Xc, yZero - (Yc-yZero), color=minusColor, width=1, clipX=(xLeftOffset, xLeftOffset + plotWidth)) - startPosX += (self.ChrLengthDistList[j]+self.GraphInterval)*plotXScale - - ###draw additive scale - if not self.multipleInterval and self.additiveChecked: - additiveScaleFont=pid.Font(ttf="verdana",size=12*fontZoom,bold=0) - additiveScale = Plot.detScaleOld(0,additiveMax) - additiveStep = (additiveScale[1]-additiveScale[0])/additiveScale[2] - additiveAxisList = Plot.frange(0, additiveScale[1], additiveStep) - maxAdd = additiveScale[1] - addPlotScale = AdditiveHeightThresh/additiveMax - - additiveAxisList.append(additiveScale[1]) - for item in additiveAxisList: - additiveY = yZero - item*addPlotScale - canvas.drawLine(xLeftOffset + plotWidth,additiveY,xLeftOffset+4+ plotWidth,additiveY,color=self.ADDITIVE_COLOR_POSITIVE, width=1*zoom) - scaleStr = "%2.3f" % item - canvas.drawString(scaleStr,xLeftOffset + plotWidth +6,additiveY+5,font=additiveScaleFont,color=self.ADDITIVE_COLOR_POSITIVE) - - canvas.drawLine(xLeftOffset+plotWidth,additiveY,xLeftOffset+plotWidth,yZero,color=self.ADDITIVE_COLOR_POSITIVE, width=1*zoom) - - canvas.drawLine(xLeftOffset, yZero, xLeftOffset, yTopOffset, color=self.LRS_COLOR, width=1*zoom) #the blue line running up the y axis - - def set_options(self, start_vars): """Sets various options (physical/genetic mapping, # permutations, which chromosome""" @@ -281,6 +75,38 @@ class IntervalMapping(object): self.do_bootstrap = start_vars['do_bootstrap'] self.control_locus = start_vars['control_locus'] self.selected_chr = start_vars['selected_chr'] + self.weighted_regression = start_vars['weighted'] + self.lrs_lod = start_vars['lrs_lod'] + + + def gen_qtl_results(self, tempdata): + """Generates qtl results for plotting interval map""" + + self.dataset.group.get_markers() + + pheno_vector = np.array([val == "x" and np.nan or float(val) for val in self.vals]) + + #if self.dataset.group.species == "human": + # p_values, t_stats = self.gen_human_results(pheno_vector, tempdata) + #else: + genotype_data = [marker['genotypes'] for marker in self.dataset.group.markers.markers] + + no_val_samples = self.identify_empty_samples() + trimmed_genotype_data = self.trim_genotypes(genotype_data, no_val_samples) + + genotype_matrix = np.array(trimmed_genotype_data).T + + t_stats, p_values = lmm.run( + pheno_vector, + genotype_matrix, + restricted_max_likelihood=True, + refit=False, + temp_data=tempdata + ) + + self.dataset.group.markers.add_pvalues(p_values) + + self.qtl_results = self.dataset.group.markers.markers def identify_empty_samples(self): @@ -290,6 +116,7 @@ class IntervalMapping(object): no_val_samples.append(sample_count) return no_val_samples + def trim_genotypes(self, genotype_data, no_value_samples): trimmed_genotype_data = [] for marker in genotype_data: @@ -304,4 +131,129 @@ class IntervalMapping(object): pass new_genotypes.append(genotype) trimmed_genotype_data.append(new_genotypes) - return trimmed_genotype_data
\ No newline at end of file + 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 = [] + |