diff options
Diffstat (limited to 'web/webqtl/heatmap/Heatmap.py')
| -rwxr-xr-x | web/webqtl/heatmap/Heatmap.py | 437 |
1 files changed, 437 insertions, 0 deletions
diff --git a/web/webqtl/heatmap/Heatmap.py b/web/webqtl/heatmap/Heatmap.py new file mode 100755 index 00000000..c4543cee --- /dev/null +++ b/web/webqtl/heatmap/Heatmap.py @@ -0,0 +1,437 @@ +# Copyright (C) University of Tennessee Health Science Center, Memphis, TN. +# +# This program is free software: you can redistribute it and/or modify it +# under the terms of the GNU Affero General Public License +# as published by the Free Software Foundation, either version 3 of the +# License, or (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. +# See the GNU Affero General Public License for more details. +# +# This program is available from Source Forge: at GeneNetwork Project +# (sourceforge.net/projects/genenetwork/). +# +# Contact Drs. Robert W. Williams and Xiaodong Zhou (2010) +# at rwilliams@uthsc.edu and xzhou15@uthsc.edu +# +# +# +# This module is used by GeneNetwork project (www.genenetwork.org) +# +# Created by GeneNetwork Core Team 2010/08/10 +# +# Last updated by GeneNetwork Core Team 2010/10/20 + +import os +import string +import piddle as pid +import cPickle + +from base import webqtlConfig +from base.webqtlTrait import webqtlTrait +from dbFunction import webqtlDatabaseFunction +from utility import webqtlUtil +from utility import Plot +import slink + + +# XZ, 09/09/2008: After adding several traits to collection, click "QTL Heatmap" button, +# XZ, 09/09/2008: This class will generate what you see. +######################################### +# QTL heatmap Page +######################################### +class Heatmap: + + def __init__(self, fd=None, searchResult=None, colorScheme=None, userPrivilege=None, userName=None): + cursor = webqtlDatabaseFunction.getCursor() + if (not cursor): + return + targetDescriptionChecked = fd.formdata.getvalue('targetDescriptionCheck', '') + clusterChecked = fd.formdata.getvalue('clusterCheck', '') + sessionfile = fd.formdata.getvalue("session") + genotype = fd.genotype + strainlist = fd.strainlist + ppolar = fd.ppolar + mpolar = fd.mpolar + traitList = [] + traitDataList = [] + for item in searchResult: + thisTrait = webqtlTrait(fullname=item, cursor=cursor) + thisTrait.retrieveInfo() + thisTrait.retrieveData(fd.strainlist) + traitList.append(thisTrait) + traitDataList.append(thisTrait.exportData(fd.strainlist)) + self.buildCanvas(colorScheme=colorScheme, targetDescriptionChecked=targetDescriptionChecked, clusterChecked=clusterChecked, sessionfile=sessionfile, genotype=genotype, strainlist=strainlist, ppolar=ppolar, mpolar=mpolar, traitList=traitList, traitDataList=traitDataList, userPrivilege=userPrivilege, userName=userName) + + def buildCanvas(self, colorScheme='', targetDescriptionChecked='', clusterChecked='', sessionfile='', genotype=None, strainlist=None, ppolar=None, mpolar=None, traitList=None, traitDataList=None, userPrivilege=None, userName=None): + labelFont = pid.Font(ttf="tahoma",size=14,bold=0) + topHeight = 0 + NNN = len(traitList) + #XZ: It's necessory to define canvas here + canvas = pid.PILCanvas(size=(80+NNN*20,880)) + names = map(webqtlTrait.displayName, traitList) + #XZ, 7/29/2009: create trait display and find max strWidth + strWidth = 0 + for j in range(len(names)): + thisTrait = traitList[j] + if targetDescriptionChecked: + if thisTrait.db.type == 'ProbeSet': + if thisTrait.probe_target_description: + names[j] += ' [%s at Chr %s @ %2.3fMB, %s]' % (thisTrait.symbol, thisTrait.chr, thisTrait.mb, thisTrait.probe_target_description) + else: + names[j] += ' [%s at Chr %s @ %2.3fMB]' % (thisTrait.symbol, thisTrait.chr, thisTrait.mb) + elif thisTrait.db.type == 'Geno': + names[j] += ' [Chr %s @ %2.3fMB]' % (thisTrait.chr, thisTrait.mb) + elif thisTrait.db.type == 'Publish': + if thisTrait.confidential: + if webqtlUtil.hasAccessToConfidentialPhenotypeTrait(privilege=userPrivilege, userName=userName, authorized_users=thisTrait.authorized_users): + if thisTrait.post_publication_abbreviation: + names[j] += ' [%s]' % (thisTrait.post_publication_abbreviation) + else: + if thisTrait.pre_publication_abbreviation: + names[j] += ' [%s]' % (thisTrait.pre_publication_abbreviation) + else: + if thisTrait.post_publication_abbreviation: + names[j] += ' [%s]' % (thisTrait.post_publication_abbreviation) + else: + pass + + i = canvas.stringWidth(names[j], font=labelFont) + if i > strWidth: + strWidth = i + + width = NNN*20 + xoffset = 40 + yoffset = 40 + cellHeight = 3 + nLoci = reduce(lambda x,y: x+y, map(lambda x: len(x),genotype),0) + + if nLoci > 2000: + cellHeight = 1 + elif nLoci > 1000: + cellHeight = 2 + elif nLoci < 200: + cellHeight = 10 + else: + pass + + pos = range(NNN) + neworder = [] + BWs = Plot.BWSpectrum() + colors100 = Plot.colorSpectrum() + colors = Plot.colorSpectrum(130) + finecolors = Plot.colorSpectrum(250) + colors100.reverse() + colors.reverse() + finecolors.reverse() + + scaleFont=pid.Font(ttf="tahoma",size=10,bold=0) + + if not clusterChecked: #XZ: this part is for original order + for i in range(len(names)): + neworder.append((xoffset+20*(i+1), i)) + + canvas = pid.PILCanvas(size=(80+NNN*20+240,80+ topHeight +5+5+strWidth+nLoci*cellHeight+80+20*cellHeight)) + + self.drawTraitNameBottom(canvas,names,yoffset,neworder,strWidth,topHeight,labelFont) + else: #XZ: this part is to cluster traits + topHeight = 400 + canvas = pid.PILCanvas(size=(80+NNN*20+240,80+ topHeight +5+5+strWidth+nLoci*cellHeight+80+20*cellHeight)) + + corArray = [([0] * (NNN))[:] for i in range(NNN)] + + nnCorr = len(strainlist) + + #XZ, 08/04/2009: I commented out pearsonArray, spearmanArray + for i, thisTrait in enumerate(traitList): + names1 = [thisTrait.db.name, thisTrait.name, thisTrait.cellid] + for j, thisTrait2 in enumerate(traitList): + names2 = [thisTrait2.db.name, thisTrait2.name, thisTrait2.cellid] + if j < i: + corr,nOverlap = webqtlUtil.calCorrelation(traitDataList[i], traitDataList[j],nnCorr) + if (1-corr) < 0: + distance = 0.0 + else: + distance = 1-corr + corArray[i][j] = distance + corArray[j][i] = distance + elif j == i: + corArray[i][j] = 0.0 + else: + pass + + #XZ, 7/29/2009: The parameter d has info of cluster (group member and distance). The format of d is tricky. Print it out to see it's format. + d = slink.slink(corArray) + + #XZ, 7/29/2009: Attention: The 'neworder' is changed by the 'draw' function + #XZ, 7/30/2009: Only toppos[1][0] and top[1][1] are used later. Then what toppos[0] is used for? + toppos = self.draw(canvas,names,d,xoffset,yoffset,neworder,topHeight) + self.drawTraitNameTop(canvas,names,yoffset,neworder,strWidth,topHeight,labelFont) + + #XZ, 7/29/2009: draw the top vertical line + canvas.drawLine(toppos[1][0],toppos[1][1],toppos[1][0],yoffset) + + #XZ: draw string 'distance = 1-r' + canvas.drawString('distance = 1-r',neworder[-1][0] + 50, topHeight*3/4,font=labelFont,angle=90) + + #draw Scale + scaleFont=pid.Font(ttf="tahoma",size=10,bold=0) + x = neworder[-1][0] + canvas.drawLine(x+5, topHeight+yoffset, x+5, yoffset, color=pid.black) + y = 0 + while y <=2: + canvas.drawLine(x+5, topHeight*y/2.0+yoffset, x+10, topHeight*y/2.0+yoffset) + canvas.drawString('%2.1f' % (2-y), x+12, topHeight*y/2.0+yoffset, font=scaleFont) + y += 0.5 + + + chrname = 0 + chrnameFont=pid.Font(ttf="tahoma",size=24,bold=0) + Ncol = 0 + + nearestMarkers = self.getNearestMarker(traitList, genotype) + + # import cPickle + if sessionfile: + fp = open(os.path.join(webqtlConfig.TMPDIR, sessionfile + '.session'), 'rb') + permData = cPickle.load(fp) + fp.close() + else: + permData = {} + + areas = [] + #XZ, 7/31/2009: This for loop is to generate the heatmap + #XZ: draw trait by trait instead of marker by marker + for order in neworder: + #startHeight = 40+400+5+5+strWidth + startHeight = topHeight + 40+5+5+strWidth + startWidth = order[0]-5 + if Ncol and Ncol % 5 == 0: + drawStartPixel = 8 + else: + drawStartPixel = 9 + + tempVal = traitDataList[order[1]] + _vals = [] + _strains = [] + for i in range(len(strainlist)): + if tempVal[i] != None: + _strains.append(strainlist[i]) + _vals.append(tempVal[i]) + + qtlresult = genotype.regression(strains = _strains, trait = _vals) + + if sessionfile: + LRSArray = permData[str(traitList[order[1]])] + else: + LRSArray = genotype.permutation(strains = _strains, trait = _vals, nperm = 1000) + permData[str(traitList[order[1]])] = LRSArray + + sugLRS = LRSArray[369] + sigLRS = LRSArray[949] + prechr = 0 + chrstart = 0 + nearest = nearestMarkers[order[1]] + midpoint = [] + + for item in qtlresult: + if item.lrs > webqtlConfig.MAXLRS: + adjustlrs = webqtlConfig.MAXLRS + else: + adjustlrs = item.lrs + + if item.locus.chr != prechr: + if prechr: + canvas.drawRect(startWidth-drawStartPixel, startHeight, startWidth+10, startHeight+3,edgeColor=pid.white, edgeWidth=0, fillColor=pid.white) + startHeight+= 3 + if not chrname: + canvas.drawString(prechr,xoffset-20,(chrstart+startHeight)/2,font = chrnameFont,color=pid.dimgray) + prechr = item.locus.chr + chrstart = startHeight + if colorScheme == '0': + if adjustlrs <= sugLRS: + colorIndex = int(65*adjustlrs/sugLRS) + else: + colorIndex = int(65 + 35*(adjustlrs-sugLRS)/(sigLRS-sugLRS)) + if colorIndex > 99: + colorIndex = 99 + colorIndex = colors100[colorIndex] + elif colorScheme == '1': + sugLRS = LRSArray[369]/2.0 + if adjustlrs <= sugLRS: + colorIndex = BWs[20+int(50*adjustlrs/sugLRS)] + else: + if item.additive > 0: + colorIndex = int(80 + 50*(adjustlrs-sugLRS)/(sigLRS-sugLRS)) + else: + colorIndex = int(50 - 50*(adjustlrs-sugLRS)/(sigLRS-sugLRS)) + if colorIndex > 129: + colorIndex = 129 + if colorIndex < 0: + colorIndex = 0 + colorIndex = colors[colorIndex] + elif colorScheme == '2': + if item.additive > 0: + colorIndex = int(150 + 100*(adjustlrs/sigLRS)) + else: + colorIndex = int(100 - 100*(adjustlrs/sigLRS)) + if colorIndex > 249: + colorIndex = 249 + if colorIndex < 0: + colorIndex = 0 + colorIndex = finecolors[colorIndex] + else: + colorIndex = pid.white + + if startHeight > 1: + canvas.drawRect(startWidth-drawStartPixel, startHeight, startWidth+10, startHeight+cellHeight,edgeColor=colorIndex, edgeWidth=0, fillColor=colorIndex) + else: + canvas.drawLine(startWidth-drawStartPixel, startHeight, startWidth+10, startHeight, Color=colorIndex) + + if item.locus.name == nearest: + midpoint = [startWidth,startHeight-5] + startHeight+=cellHeight + + #XZ, map link to trait name and band + COORDS = "%d,%d,%d,%d" %(startWidth-drawStartPixel,topHeight+40,startWidth+10,startHeight) + HREF = "javascript:showDatabase2('%s','%s','%s');" % (traitList[order[1]].db.name, traitList[order[1]].name, traitList[order[1]].cellid) + area = (COORDS, HREF, '%s' % names[order[1]]) + areas.append(area) + + if midpoint: + traitPixel = ((midpoint[0],midpoint[1]),(midpoint[0]-6,midpoint[1]+12),(midpoint[0]+6,midpoint[1]+12)) + canvas.drawPolygon(traitPixel,edgeColor=pid.black,fillColor=pid.orange,closed=1) + + if not chrname: + canvas.drawString(prechr,xoffset-20,(chrstart+startHeight)/2,font = chrnameFont,color=pid.dimgray) + chrname = 1 + Ncol += 1 + + + #draw Spectrum + startSpect = neworder[-1][0] + 30 + startHeight = topHeight + 40+5+5+strWidth + + if colorScheme == '0': + for i in range(100): + canvas.drawLine(startSpect+i,startHeight+20,startSpect+i,startHeight+40,color=colors100[i]) + scaleFont=pid.Font(ttf="tahoma",size=10,bold=0) + canvas.drawLine(startSpect,startHeight+45,startSpect,startHeight+39,color=pid.black) + canvas.drawString('LRS = 0',startSpect,startHeight+55,font=scaleFont) + canvas.drawLine(startSpect+64,startHeight+45,startSpect+64,startHeight+39,color=pid.black) + canvas.drawString('Suggestive LRS',startSpect+64,startHeight+55,font=scaleFont) + canvas.drawLine(startSpect+99,startHeight+45,startSpect+99,startHeight+39,color=pid.black) + canvas.drawString('Significant LRS',startSpect+105,startHeight+40,font=scaleFont) + elif colorScheme == '1': + for i in range(50): + canvas.drawLine(startSpect+i,startHeight,startSpect+i,startHeight+40,color=BWs[20+i]) + for i in range(50,100): + canvas.drawLine(startSpect+i,startHeight,startSpect+i,startHeight+20,color=colors[100-i]) + canvas.drawLine(startSpect+i,startHeight+20,startSpect+i,startHeight+40,color=colors[30+i]) + + canvas.drawLine(startSpect,startHeight+45,startSpect,startHeight+39,color=pid.black) + canvas.drawString('LRS = 0',startSpect,startHeight+60,font=scaleFont) + canvas.drawLine(startSpect+50,startHeight+45,startSpect+50,startHeight+39,color=pid.black) + canvas.drawString('0.5*Suggestive LRS',startSpect+50,startHeight+ 60,font=scaleFont) + canvas.drawLine(startSpect+99,startHeight+45,startSpect+99,startHeight+39,color=pid.black) + canvas.drawString('Significant LRS',startSpect+105,startHeight+50,font=scaleFont) + textFont=pid.Font(ttf="verdana",size=18,bold=0) + canvas.drawString('%s +' % ppolar,startSpect+120,startHeight+ 35,font=textFont,color=pid.red) + canvas.drawString('%s +' % mpolar,startSpect+120,startHeight+ 15,font=textFont,color=pid.blue) + elif colorScheme == '2': + for i in range(100): + canvas.drawLine(startSpect+i,startHeight,startSpect+i,startHeight+20,color=finecolors[100-i]) + canvas.drawLine(startSpect+i,startHeight+20,startSpect+i,startHeight+40,color=finecolors[150+i]) + + canvas.drawLine(startSpect,startHeight+45,startSpect,startHeight+39,color=pid.black) + canvas.drawString('LRS = 0',startSpect,startHeight+60,font=scaleFont) + canvas.drawLine(startSpect+99,startHeight+45,startSpect+99,startHeight+39,color=pid.black) + canvas.drawString('Significant LRS',startSpect+105,startHeight+50,font=scaleFont) + textFont=pid.Font(ttf="verdana",size=18,bold=0) + canvas.drawString('%s +' % ppolar,startSpect+120,startHeight+ 35,font=textFont,color=pid.red) + canvas.drawString('%s +' % mpolar,startSpect+120,startHeight+ 15,font=textFont,color=pid.blue) + + filename= webqtlUtil.genRandStr("Heatmap_") + canvas.save(webqtlConfig.IMGDIR+filename, format='png') + if not sessionfile: + sessionfile = webqtlUtil.generate_session() + webqtlUtil.dump_session(permData, os.path.join(webqtlConfig.TMPDIR, sessionfile +'.session')) + self.filename=filename + self.areas=areas + self.sessionfile=sessionfile + + def getResult(self): + return self.filename, self.areas, self.sessionfile + + #XZ, 7/31/2009: This function put the order of traits into parameter neworder, + #XZ: return the position of the top vertical line of the hierarchical tree, draw the hierarchical tree. + def draw(self,canvas,names,d,xoffset,yoffset,neworder,topHeight): + maxDistance = topHeight + fontoffset = 4 #XZ, 7/31/2009: used only for drawing tree + if type(d[0]) == type(1) and type(d[1]) == type(1): + neworder.append((xoffset+20,d[0])) + neworder.append((xoffset+40,d[1])) + height = d[2]*maxDistance/2 + canvas.drawLine(xoffset+20-fontoffset,maxDistance+yoffset,xoffset+20-fontoffset,maxDistance-height+yoffset) + canvas.drawLine(xoffset+40-fontoffset,maxDistance+yoffset,xoffset+40-fontoffset,maxDistance-height+yoffset) + canvas.drawLine(xoffset+40-fontoffset,maxDistance+yoffset-height,xoffset+20-fontoffset,maxDistance-height+yoffset) + return (40,(xoffset+30-fontoffset,maxDistance-height+yoffset)) + elif type(d[0]) == type(1): + neworder.append((xoffset+20,d[0])) + d2 = self.draw(canvas,names,d[1],xoffset+20,yoffset,neworder,topHeight) + height = d[2]*maxDistance/2 + canvas.drawLine(xoffset+20-fontoffset,maxDistance+yoffset,xoffset+20-fontoffset,maxDistance-height+yoffset) + canvas.drawLine(d2[1][0],d2[1][1],d2[1][0],maxDistance-height+yoffset) + canvas.drawLine(d2[1][0],maxDistance-height+yoffset,xoffset+20-fontoffset,maxDistance-height+yoffset) + return (20+d2[0],((d2[1][0]+xoffset+20-fontoffset)/2,maxDistance-height+yoffset)) + elif type(d[1]) == type(1): + d1 = self.draw(canvas,names,d[0],xoffset,yoffset,neworder,topHeight) + neworder.append((xoffset+d1[0]+20,d[1])) + height = d[2]*maxDistance/2 + canvas.drawLine(xoffset+d1[0]+20-fontoffset,maxDistance+yoffset,xoffset+d1[0]+20-fontoffset,maxDistance-height+yoffset) + canvas.drawLine(d1[1][0],d1[1][1],d1[1][0],maxDistance-height+yoffset) + canvas.drawLine(d1[1][0],maxDistance-height+yoffset,xoffset+d1[0]+20-fontoffset,maxDistance-height+yoffset) + return (d1[0]+20,((d1[1][0]+xoffset+d1[0]+20-fontoffset)/2,maxDistance-height+yoffset)) + else: + d1 = self.draw(canvas,names,d[0],xoffset,yoffset,neworder,topHeight) + d2 = self.draw(canvas,names,d[1],xoffset+d1[0],yoffset,neworder,topHeight) + height = d[2]*maxDistance/2 + canvas.drawLine(d2[1][0],d2[1][1],d2[1][0],maxDistance-height+yoffset) + canvas.drawLine(d1[1][0],d1[1][1],d1[1][0],maxDistance-height+yoffset) + canvas.drawLine(d1[1][0],maxDistance-height+yoffset,d2[1][0],maxDistance-height+yoffset) + return (d1[0]+d2[0],((d1[1][0]+d2[1][0])/2,maxDistance-height+yoffset)) + + #XZ, 7/31/2009: dras trait names + def drawTraitNameBottom (self,canvas,names,yoffset,neworder,strWidth,topHeight,labelFont): + maxDistance = topHeight + for oneOrder in neworder: + canvas.drawString(names[oneOrder[1]],oneOrder[0]-5,maxDistance+yoffset+5+strWidth-canvas.stringWidth(names[oneOrder[1]],font=labelFont),font=labelFont,color=pid.black,angle=270) + + def drawTraitNameTop (self,canvas,names,yoffset,neworder,strWidth,topHeight,labelFont): + maxDistance = topHeight + for oneOrder in neworder: + canvas.drawString(names[oneOrder[1]],oneOrder[0]-5,maxDistance+yoffset+5,font=labelFont,color=pid.black,angle=270) + + + def getNearestMarker(self, traitList, genotype): + out = [] + if not genotype.Mbmap: + return [None]* len(traitList) + for item in traitList: + try: + nearest = None + for _chr in genotype: + if _chr.name != item.chr: + continue + distance = 1e30 + for _locus in _chr: + if abs(_locus.Mb-item.mb) < distance: + distance = abs(_locus.Mb-item.mb) + nearest = _locus.name + out.append(nearest) + except: + out.append(None) + + return out |