1 ###########################################################################
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 #
7 # All Rights Reserved. #
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: #
17 # The above copyright notice and this permission notice shall be #
18 # included in all copies or substantial portions of the Software. #
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 #
28 ###########################################################################
30 # data for Danio Rerio
32 from cistematic.genomes import Genome
33 from cistematic.core.geneinfo import geneinfoDB
34 from os import environ
36 if environ.get("CISTEMATIC_ROOT"):
37 cisRoot = environ.get("CISTEMATIC_ROOT")
39 cisRoot = "/proj/genome"
41 geneDB = "%s/D_rerio/drerio.genedb" % cisRoot
44 def loadChromosome(db, chromPath, chromOutPath):
47 drGenome = Genome("drerio", dbFile=db)
48 files = os.listdir(chromPath)
49 for filename in files:
50 inFile = open("%s/%s" % (chromPath, filename), "r")
51 header = inFile.readline()
55 chromID = header.strip()[1:]
56 currentLine = inFile.readline()
57 while currentLine != "" and currentLine[0] != ">":
58 lineSeq = currentLine.strip()
59 seqLen += len(lineSeq)
60 seqArray.append(lineSeq)
61 currentLine = inFile.readline()
63 seq = string.join(seqArray, "")
65 print "Added contig %s to database" % chromID
66 drGenome.addSequence(("drerio", chromID), seq, "chromosome", str(seqLen))
67 drGenome.addChromosomeEntry(chromID, chromID, "db")
69 outFileName = "%s%s.bin" % (chromOutPath, chromID)
70 outFile = open("%s%s" % (cisRoot, outFileName), "w")
73 print "Added contig file %s to database" % outFileName
74 drGenome.addChromosomeEntry(chromID, outFileName, "file")
81 def loadGeneEntries(db, gFile):
84 drGenome = Genome("drerio", dbFile=db)
85 geneFile = open(gFile, "r")
88 cols = line.split("\t")
91 tempID = idb.getGeneID("drerio", gid)
93 print "could not find %s" % gid
96 geneInfo = idb.getGeneInfo(tempID)
113 geneID = ("drerio", gid)
114 if geneID in seenGIDs:
118 seenGIDs.append(geneID)
120 geneEntries.append((geneID, chrom, start, stop, sense, "gene", gidVersion))
122 print "Adding %d gene entries" % len(geneEntries)
123 drGenome.addGeneEntryBatch(geneEntries)
126 def loadGeneFeatures(db, gfile):
127 geneFile = open(gfile, "r")
130 senseArray = {"+": "F",
136 for geneLine in geneFile:
137 geneFields = geneLine.split("\t")
138 exonNum = int(geneFields[8])
139 exonStarts = geneFields[9].split(",")
140 exonStops = geneFields[10].split(",")
141 chrom = geneFields[2]
142 sense = senseArray[geneFields[3]]
143 gstart = int(geneFields[6]) - 1
144 gstop = int(geneFields[7]) - 1
147 tempID = idb.getGeneID("drerio", gid)
149 print "could not find %s" % gid
152 geneInfo = idb.getGeneInfo(tempID)
160 geneID = ("drerio", gid)
161 if geneID in seenGIDs:
166 for index in range(exonNum):
167 estart = int(exonStarts[index]) - 1
168 estop = int(exonStops[index]) - 1
169 if estart >= gstart and estop <= gstop:
170 insertArray.append((geneID, gidVersion, chrom, estart, estop, sense, "CDS"))
171 elif estop <= gstart:
177 insertArray.append((geneID, gidVersion, chrom, estart, estop, sense, fType))
178 elif estart >= gstop:
184 insertArray.append((geneID, gidVersion, chrom, estart, estop, sense, fType))
185 elif estart <= gstop and estart > gstart:
191 insertArray.append((geneID, gidVersion, chrom, estart, gstop, sense, "CDS"))
192 insertArray.append((geneID, gidVersion, chrom, gstop + 1, estop, sense, fType))
193 elif estart < gstart and estop <= gstop:
199 insertArray.append((geneID, gidVersion, chrom, estart, gstart - 1, sense, fType))
200 insertArray.append((geneID, gidVersion, chrom, gstart, estop, sense, "CDS"))
209 insertArray.append((geneID, gidVersion, chrom, estart, gstart - 1, sense, fType1))
210 insertArray.append((geneID, gidVersion, chrom, gstart, gstop, sense, "CDS"))
211 insertArray.append((geneID, gidVersion, chrom, gstop + 1, estop - 1, sense, fType2))
214 drGenome = Genome("drerio", dbFile=db)
215 print "Adding %d features" % len(insertArray)
216 drGenome.addFeatureEntryBatch(insertArray)
219 def createDBFile(db):
220 drGenome = Genome("drerio", dbFile=db)
221 drGenome.createGeneDB(db)
224 def createDBindices(db):
225 drGenome = Genome("drerio", dbFile=db)
226 drGenome.createIndices()
229 def buildDrerioDB(db=geneDB):
230 """ genes and annotations are from UCSC (dr3).
232 #genePath = "%s/download/xenoRefFlat.txt" % cisRoot
233 chromoPath = "%s/download/dr3" % cisRoot
234 chromoOutPath = "/D_rerio/"
235 print "Creating database %s" % db
238 #print "Adding gene entries"
239 #loadGeneEntries(db, genePath)
241 #print "Adding gene features"
242 #loadGeneFeatures(db, genePath)
244 print "Loading chromosomes"
245 loadChromosome(db, chromoPath, chromoOutPath)
247 print "Creating Indices"
250 print "Finished creating database %s" % db