cistematic/genomes/btaurus.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 # data for Bos Taurus
  31 import string
  32 from cistematic.genomes import Genome
  33 from os import environ
  34
  35 if environ.get("CISTEMATIC_ROOT"):
  36     cisRoot = environ.get("CISTEMATIC_ROOT")
  37 else:
  38     cisRoot = "/proj/genome"
  39
  40 geneDB = "%s/B_taurus/btaurus.genedb" % cisRoot
  41
  42
  43 def loadChromosome(db, chromPath, chromOutPath):
  44     seqArray = []
  45     seqLen = 0
  46     btGenome = Genome("btaurus", dbFile=db)
  47     inFile = open(chromPath, "r")
  48     header = inFile.readline()
  49     while header != "":
  50         seqArray = []
  51         seqLen = 0
  52         chromID = header.strip()[1:]
  53         currentLine = inFile.readline()
  54
  55         while currentLine != "" and currentLine[0] != ">":
  56             lineSeq = currentLine.strip()
  57             seqLen += len(lineSeq)
  58             seqArray.append(lineSeq)
  59             currentLine = inFile.readline()
  60
  61         seq = string.join(seqArray, "")
  62         if seqLen < 500000:
  63             print "Added contig %s to database" % chromID
  64             btGenome.addSequence(("btaurus", chromID), seq, "chromosome", str(seqLen))
  65             btGenome.addChromosomeEntry(chromID, chromID, "db")
  66         else:
  67             outFileName = "%s%s.bin" % (chromOutPath, chromID)
  68             outFile = open("%s%s" % (cisRoot, outFileName), "w")
  69             outFile.write(seq)
  70             outFile.close()
  71             print "Added contig file %s to database" % outFileName
  72             btGenome.addChromosomeEntry(chromID, outFileName, "file")
  73
  74         header = currentLine
  75
  76     inFile.close()
  77
  78
  79 def loadGeneEntries(db, gFile):
  80     """ FIXME - NEED TO DEAL WITH ALTERNATIVE SPLICING ENTRIES
  81     """
  82     geneEntries = []
  83     btGenome = Genome("btaurus", dbFile=db)
  84     geneFile = open(gFile, "r")
  85
  86     for line in geneFile:
  87         cols = line.split("\t")
  88         gid = cols[0]
  89         start = int(cols[5])
  90         stop = int(cols[6])
  91         sense = cols[2]
  92         chrom = cols[1]
  93         if sense == "+":
  94             sense = "F"
  95         else:
  96             sense = "R"
  97
  98         geneID = ("btaurus", gid)
  99         gidVersion = 1
 100         geneEntries.append((geneID, chrom, start, stop, sense, "gene", gidVersion))
 101
 102     print "Adding %d gene entries" % len(geneEntries)
 103     btGenome.addGeneEntryBatch(geneEntries)
 104
 105
 106 def loadGeneAnnotations(db, annotPath):
 107     geneAnnotations = []
 108     annotFile = open(annotPath, "r")
 109     btGenome = Genome("btaurus", dbFile=db)
 110     for line in annotFile:
 111         try:
 112             cols = line.split("\t")
 113             locID = cols[0]
 114             geneDesc = cols[6]
 115             if len(locID) > 0:
 116                 geneAnnotations.append((("btaurus", locID), string.replace(geneDesc.strip(), "'", "p")))
 117         except:
 118             pass
 119
 120     print "Adding %d annotations" % len(geneAnnotations)
 121     btGenome.addAnnotationBatch(geneAnnotations)
 122
 123
 124 def loadGeneFeatures(db, gfile):
 125     geneFile = open(gfile, "r")
 126     senseArray = {"+": "F",
 127                   "-": "R",
 128                   ".": "F"
 129     }
 130
 131     seenArray = []
 132     insertArray = []
 133     for geneLine in geneFile:
 134         geneFields = geneLine.split("\t")
 135         exonNum = int(geneFields[7])
 136         exonStarts = geneFields[8].split(",")
 137         exonStops = geneFields[9].split(",")
 138         chrom = geneFields[1]
 139         sense = senseArray[geneFields[2]]
 140         gstop = int(geneFields[6]) - 1
 141         gstart = int(geneFields[5]) - 1
 142         geneid = geneFields[0]
 143         try:
 144             geneID = ("btaurus", geneid)
 145         except:
 146             continue
 147
 148         gidVersion = "1"
 149         if geneID in seenArray:
 150             gidVersion = "2" # doesn't deal with more than 2 refseq's for the same locus, yet.
 151         else:
 152             seenArray.append(geneID)
 153
 154         for index in range(exonNum):
 155             estart = int(exonStarts[index]) - 1
 156             estop = int(exonStops[index]) - 1
 157             if estart >= gstart and estop <= gstop:
 158                 insertArray.append((geneID, gidVersion, chrom, estart, estop, sense, "CDS"))
 159             elif estop <= gstart:
 160                 if sense == "F":
 161                     fType = "5UTR"
 162                 else:
 163                     fType = "3UTR"
 164
 165                 insertArray.append((geneID, 1, chrom, estart, estop, sense, fType))
 166             elif estart >= gstop:
 167                 if sense == "F":
 168                     fType = "3UTR"
 169                 else:
 170                     fType = "5UTR"
 171
 172                 insertArray.append((geneID, 1, chrom, estart, estop, sense, fType))
 173             elif estart <= gstop:
 174                 if sense == "F":
 175                     fType = "3UTR"
 176                 else:
 177                     fType = "5UTR"
 178
 179                 insertArray.append((geneID, 1, chrom, estart, gstop, sense, "CDS"))
 180                 insertArray.append((geneID, 1, chrom, gstop + 1, estop, sense, fType))
 181             else:
 182                 if sense == "F":
 183                     fType = "5UTR"
 184                 else:
 185                     fType = "3UTR"
 186
 187                 insertArray.append((geneID, 1, chrom, gstart, estop, sense, "CDS"))
 188                 insertArray.append((geneID, 1, chrom, estart, gstart - 1, sense, fType))
 189
 190     geneFile.close()
 191
 192     btGenome = Genome("btaurus", dbFile=db)
 193     print 'Adding %d features' % len(insertArray)
 194     btGenome.addFeatureEntryBatch(insertArray)
 195
 196
 197 def loadGeneOntology(db, goPath, goDefPath, annotPath):
 198     btGenome = Genome("btaurus", dbFile=db)
 199     goDefFile = open(goDefPath, "r")
 200     goFile = open(goPath, "r")
 201     annotFile = open(annotPath, "r")
 202     annotEntries = annotFile.readlines()
 203     annotFile.close()
 204     goDefEntries = goDefFile.readlines()
 205     goDefs = {}
 206     locus = {}
 207     goArray = []
 208
 209     for goDefEntry in goDefEntries:
 210         if goDefEntry[0] != "!":
 211             cols = goDefEntry.split("\t")
 212             goDefs[cols[0]] = (cols[1], cols[2].strip())
 213
 214     goEntries = goFile.readlines()
 215     for annotEntry in annotEntries:
 216         try:
 217             cols = annotEntry.split("\t")
 218             locID = cols[0]
 219             geneName = cols[1]
 220             geneDesc = cols[6]
 221             mimID = ""
 222             if len(locID) > 0:
 223                 locus[locID] = (geneName, geneDesc, mimID)
 224         except:
 225             pass
 226
 227     for entry in goEntries:
 228         try:
 229             fields = entry.split("\t")
 230             locID = fields[0].strip()
 231             (gene_name, gene_desc, mimID) = locus[locID]
 232             goArray.append((("btaurus", locID), fields[1], "", gene_name, "", string.replace(goDefs[fields[1]][0], "'", "p"), goDefs[fields[1]][1], mimID))
 233         except:
 234             #print "locus ID %s could not be added" % locID
 235             pass
 236
 237     print "adding %d go entries" % len(goArray)
 238     btGenome.addGoInfoBatch(goArray)
 239
 240
 241 def createDBFile(db):
 242     btGenome = Genome("btaurus",  dbFile=db)
 243     btGenome.createGeneDB(db)
 244
 245
 246 def createDBindices(db):
 247     btGenome = Genome("btaurus", dbFile=db)
 248     btGenome.createIndices()
 249
 250
 251 def buildBtaurusDB(db=geneDB):
 252     genePath = "%s/download/bt2/genscan.txt" % cisRoot
 253     chromoPath = "%s/download/bt2/bosTau2.softmask2.fa" % cisRoot
 254     chromoOutPath = "/B_taurus/"
 255
 256     print "Creating database %s" % db
 257     createDBFile(db)
 258
 259     print "Adding gene entries"
 260     loadGeneEntries(db, genePath)
 261
 262     print "Adding gene features"
 263     loadGeneFeatures(db, genePath)
 264
 265     print "Loading sequences"
 266     loadChromosome(db, chromoPath, chromoOutPath)
 267
 268     print "Creating Indices"
 269     createDBindices(db)
 270
 271     print "Finished creating database %s" % db