cistematic/experiments/draw.py

   1 ###########################################################################
   2 #                                                                         #
   3 # C O P Y R I G H T   N O T I C E                                         #
   4 #  Copyright (c) 2003-10 by:                                              #
   5 #    * California Institute of Technology                                 #
   6 #                                                                         #
   7 #    All Rights Reserved.                                                 #
   8 #                                                                         #
   9 # Permission is hereby granted, free of charge, to any person             #
  10 # obtaining a copy of this software and associated documentation files    #
  11 # (the "Software"), to deal in the Software without restriction,          #
  12 # including without limitation the rights to use, copy, modify, merge,    #
  13 # publish, distribute, sublicense, and/or sell copies of the Software,    #
  14 # and to permit persons to whom the Software is furnished to do so,       #
  15 # subject to the following conditions:                                    #
  16 #                                                                         #
  17 # The above copyright notice and this permission notice shall be          #
  18 # included in all copies or substantial portions of the Software.         #
  19 #                                                                         #
  20 # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,         #
  21 # EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF      #
  22 # MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND                   #
  23 # NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS     #
  24 # BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN      #
  25 # ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN       #
  26 # CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE        #
  27 # SOFTWARE.                                                               #
  28 ###########################################################################
  29 #
  30 # This class contains the core code for visualizing sequences from several of the experiments.
  31
  32 from cistematic.core.motif import Motif
  33 from cistematic.core.geneinfo import geneinfoDB
  34 import math, random
  35
  36 pilPresent = False
  37
  38 try:
  39     import Image, ImageDraw, ImageFont
  40     pilPresent = True
  41 except:
  42     pass
  43
  44
  45 class Draw:
  46     """ The Draw class contains the code used to visualize the location of motifs in their
  47         genomic neighborhood. It is meant to be used as the parent of other classes,
  48         such as the orthology classes. It relies on the python imaging library, and saves
  49         PNG images.
  50     """
  51     drawable = False
  52     maxWidth = 1200
  53     leftMargin = 100
  54     rightMargin = 50
  55     lineHeight = 110
  56
  57     if pilPresent:
  58         theFont = ImageFont.load_default()
  59         drawable = True
  60     else:
  61         print "Draw: python image library missing - will not be able to draw on this system"
  62
  63
  64     def draw(self, picName, geneList=[], motifList=[], excludeGeneList=[], excludeMotifList=[],
  65             showHeader=True, showFooter=True, maxOccurences=100, skipSanity=False):
  66         """ Draws an image of motifs on the sequences. Can specifically list or exclude genes and/or motifs.
  67             Options showHeader will add experiment information while showFooter will add a motif key.
  68             This function will not show motifs that occur more than maxOccurences in the dataset.
  69         """
  70         if not self.drawable:
  71             return
  72
  73         limitGenes = False
  74         limitMotifs = False
  75         if len(geneList) > 0:
  76             limitGenes = True
  77
  78         if len(motifList) > 0:
  79             limitMotifs = True
  80
  81         idb = geneinfoDB()
  82         bound = ""
  83         if self.boundToNextGene:
  84             bound = "up to "
  85
  86         expResults = self.getResults()
  87         (up, cds, down) = self.getSeqParameters()
  88         hasORF = "NO"
  89         if cds > 0:
  90             hasORF = "YES"
  91
  92         if cds > 1:
  93             hasORF = "MASKED"
  94
  95         motColor = {}
  96         motConsensus = {}
  97         motNumber = {}
  98         maxLength = 1
  99         geneLength = {}
 100         adjustedMotLength = {}
 101         datasetIDs = self.getDatasetNames()
 102         orthologyList = []
 103         if len(datasetIDs) > 0:
 104             for datasetName in datasetIDs:
 105                 dataset = self.getSetting(datasetName)
 106                 theList = eval(dataset[0])
 107                 for gene in theList:
 108                     if limitGenes and gene not in geneList:
 109                         continue
 110
 111                     if gene not in excludeGeneList and gene not in orthologyList:
 112                         orthologyList.append(gene)
 113         else:
 114             for gene in self.genepool:
 115                 if limitGenes and gene not in geneList:
 116                     continue
 117
 118                 if gene not in excludeGeneList and gene not in orthologyList:
 119                     orthologyList.append(gene)
 120
 121         for geneID in orthologyList:
 122             geneLength[geneID] = len(self.genepool[geneID])
 123             if maxLength < geneLength[geneID]:
 124                 maxLength = geneLength[geneID]
 125
 126         seqScaler = float(self.maxWidth) / float(maxLength)
 127         for mot in expResults:
 128             if limitMotifs and mot.tagID not in motifList:
 129                 continue
 130
 131             if mot.tagID in excludeMotifList:
 132                 continue
 133
 134             if len(self.motifToGene(mot.tagID)) >= maxOccurences:
 135                 continue
 136
 137             if not skipSanity and not mot.isSane():
 138                 continue
 139
 140             motConsensus[mot.tagID] = mot.buildConsensus()
 141             currentRed = random.randint(5, 240)
 142             currentGreen = random.randint(5, 240)
 143             currentBlue = random.randint(5, 240)
 144             motColor[mot.tagID] = (currentRed, currentGreen, currentBlue)
 145
 146         motKeys = motConsensus.keys()
 147         motKeys.sort()
 148         motKeysLen = len(motKeys)
 149         for tagID in motKeys:
 150             motNumber[tagID] = 0
 151             adjustedMotLength[tagID] = int(math.ceil(len(self.findMotif(tagID)) * seqScaler))
 152
 153         numLines = len(orthologyList)
 154         if showHeader:
 155             numLines += 1
 156
 157         footerLines = motKeysLen / 3
 158         if motKeysLen % 3:
 159             footerLines += 1
 160
 161         imsize = (self.maxWidth + self.leftMargin + self.rightMargin, int(round((numLines + footerLines/2.) * self.lineHeight)))
 162         image = Image.new("RGB", imsize, color="#ffffff")
 163         draw = ImageDraw.Draw(image)
 164         currentHeight = 0
 165         if showHeader:
 166             line1 = "Experiment: %s in %s Type: %s Analysis: %s" % (self.experimentID, self.expFile, self.experimentType, self.analysisID)
 167             draw.text([10, 10], line1, font=self.theFont, fill=0)
 168             draw.line((10, 30, self.maxWidth + self.leftMargin + self.rightMargin - 10, 30), fill=0)
 169             line2 = "Upstream: %s%s ORF: %s Downstream: %s%s" % (bound, up, hasORF, bound, down)
 170             draw.text([10, 40], line2, font=self.theFont, fill=0)
 171             draw.line((10, 60, self.maxWidth + self.leftMargin + self.rightMargin - 10, 60), fill=0)
 172             currentHeight = self.lineHeight
 173
 174         for geneID in orthologyList:
 175             geneNames = ""
 176             seqLength = geneLength[geneID]
 177             try:
 178                 res = idb.geneIDSynonyms(geneID)
 179                 for entry in res[1:]:
 180                     geneNames += "%s " % str(entry)
 181             except:
 182                 pass
 183
 184             motList = self.geneToMotif(geneID)
 185             draw.text([5, currentHeight+45], str(geneID[0]), font=self.theFont, fill=0)
 186             draw.text([5, currentHeight+55], str(geneID[1]), font=self.theFont, fill=0)
 187             draw.text([5, currentHeight+65], geneNames, font=self.theFont, fill=0)
 188             adjustedSeqLength = seqLength * seqScaler
 189             seqStart = self.leftMargin + int(self.maxWidth) - adjustedSeqLength
 190             features = self.getFeatures(geneID)
 191             for (ftype, fstart, fstop, forientation) in features:
 192                 if ftype != "CDS":
 193                     continue
 194
 195                 if float(fstop) < float(fstart):
 196                     fstart = fstop
 197                     fstop = fstart
 198
 199                 start = int(math.floor(float(fstart) * seqScaler))
 200                 consLength = int(math.ceil((fstop - fstart) * seqScaler))
 201                 if start + consLength > adjustedSeqLength:
 202                     consLength = adjustedSeqLength - start
 203
 204                 start += seqStart
 205                 draw.rectangle([start, currentHeight + 51, start + consLength, currentHeight + 64], fill="#aaaaaa")
 206
 207             draw.rectangle([seqStart, currentHeight + 35, seqStart + adjustedSeqLength, currentHeight + 65], outline=0)
 208             tagIndex = 0
 209             tagPosList = []
 210             for (tagID, (pos, sense)) in motList:
 211                 if tagID not in motKeys:
 212                     continue
 213
 214                 tagIndex += 1
 215                 motNumber[tagID] += 1
 216                 start = 0
 217                 top = 0
 218                 bottom = 0
 219                 start = int(math.floor(float(pos) * seqScaler)) + seqStart
 220                 if sense == "F":
 221                     top = 12
 222                     if tagIndex % 2:
 223                         textHeight = 0
 224                     else:
 225                         textHeight = 70
 226                 else:
 227                     top = 35
 228                     bottom = 18
 229                     if tagIndex % 2:
 230                         textHeight = 15
 231                     else:
 232                         textHeight = 85
 233
 234                 for (prevStart, prevHeight) in tagPosList:
 235                     if abs(prevStart - start) <= 7 and abs(prevHeight - textHeight) <= 7:
 236                         if textHeight < 50:
 237                             textHeight -= 7
 238                         else:
 239                             textHeight += 7
 240
 241                 tagPosList.append((start, textHeight))
 242                 draw.rectangle([start, currentHeight + top, start + adjustedMotLength[tagID], currentHeight + 65 + bottom], fill=motColor[tagID])
 243                 if tagID.count("-") == 2:
 244                     tagIDlist = tagID.split("-")
 245                     tempID = tagIDlist[0] + tagIDlist[1][0] + tagIDlist[2]
 246                 else:
 247                     tempID = tagID
 248
 249                 draw.text([start - 5, currentHeight + textHeight], tempID, font=self.theFont, fill=motColor[tagID])
 250
 251             conservedWindows = []
 252             try:
 253                 conservedWindows = self.getConservedSequenceWindows(geneID)
 254             except:
 255                 pass
 256
 257             for (location, cLength, criteria) in conservedWindows:
 258                 start = int(math.floor(float(location) * seqScaler)) + seqStart
 259                 consLength = int(math.ceil(float(cLength) * seqScaler))
 260                 draw.rectangle([start, currentHeight + 38, start + consLength, currentHeight + 50], fill ='#ff0000')
 261
 262             for location in range(0, seqLength, 1000):
 263                 start = int(math.floor(float(location) * seqScaler)) + seqStart
 264                 draw.rectangle([start, currentHeight + 60, start + 1, currentHeight + 65], fill = '#000000')
 265
 266             draw.text([self.leftMargin + self.maxWidth + 5, currentHeight + 45], str(seqLength), font=self.theFont, fill=0)
 267             currentHeight += self.lineHeight
 268
 269         motNum = 0
 270         if showFooter:
 271             for motID in motKeys:
 272                 x = self.leftMargin + (self.maxWidth / 3) * (motNum % 3)
 273                 x = self.leftMargin + (self.maxWidth / 3) * (motNum % 3)
 274                 y1 = currentHeight + ((self.lineHeight/2) * (motNum/3))+ 5
 275                 y2 = currentHeight + ((self.lineHeight/2) * (motNum/3))+ 15
 276                 y3 = currentHeight + ((self.lineHeight/2) * (motNum/3))+ 25
 277                 draw.text([x, y1], motID, font=self.theFont, fill=motColor[motID])
 278                 draw.text([x, y2], motConsensus[motID], font=self.theFont, fill=0)
 279                 draw.text([x, y3], str(motNumber[motID]) + " matches", font=self.theFont, fill=0)
 280                 motNum +=1
 281
 282         del draw
 283         image.save(picName, "PNG")
 284
 285
 286     def drawMotifs(self, picName, motifList, geneList=[], excludeGeneList=[], showHeader=True,
 287                    showFooter=True, maxOccurences=100, genesWithMotifOnly=True, skipSanity=False):
 288         """ Draws an image of one or more motifs on the sequences. Can specifically list or exclude genes.
 289             Options showHeader will add experiment information while showFooter will add a motif key.
 290             This function will not show motifs that occur more than maxOccurences in the dataset, and will
 291             only show sequences with the motif by default.
 292         """
 293         restrictedGeneList = []
 294         if genesWithMotifOnly:
 295             for motID in motifList:
 296                 matches = self.motifToGene(motID)
 297                 for (loc, pos) in matches:
 298                     if loc not in restrictedGeneList and loc not in excludeGeneList:
 299                         restrictedGeneList.append(loc)
 300         else:
 301             restrictedGeneList = geneList
 302
 303         self.draw(picName, restrictedGeneList, motifList, excludeGeneList, [], showHeader, showFooter, maxOccurences, skipSanity)
 304
 305
 306     def drawGenes(self, picName, geneList, motifList=[], excludeMotifList=[], showHeader=True,
 307                   showFooter=True, maxOccurences=100, includeHomologs=False, motifsOnGeneOnly=True,
 308                   skipSanity=False):
 309         """ Draws an image of one or more motifs on the sequences. Can specifically list or exclude motif.
 310             Options showHeader will add experiment information while showFooter will add a motif key.
 311             This function will not show motifs that occur more than maxOccurences in the dataset, and will
 312             only show motifs on the sequence by default in the footer.
 313         """
 314         restrictedMotifList = []
 315         theGeneList = []
 316         if includeHomologs:
 317             for geneID in geneList:
 318                 if geneID not in theGeneList:
 319                     theGeneList.append(geneID)
 320                 try:
 321                     hgenes = self.returnHomologs(geneID)
 322                     for gID in hgenes:
 323                         if gID in self.genepool and gID not in theGeneList:
 324                             theGeneList.append(gID)
 325                 except:
 326                     pass
 327         else:
 328             theGeneList = geneList
 329
 330         if motifsOnGeneOnly:
 331             for geneID in theGeneList:
 332                 matches = self.geneToMotif(geneID)
 333                 for (motID, pos) in matches:
 334                     if motID not in restrictedMotifList and motID not in excludeMotifList:
 335                         restrictedMotifList.append(motID)
 336         else:
 337             restrictedMotifList = motifList
 338
 339         self.draw(picName, theGeneList, restrictedMotifList, [], excludeMotifList, showHeader, showFooter, maxOccurences, skipSanity)