1 ###########################################################################
3 # C O P Y R I G H T N O T I C E #
4 # Copyright (c) 2003-13 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 buildDrerioDB(db=geneDB):
45 """ genes and annotations are from UCSC (dr3).
47 #genePath = "%s/download/xenoRefFlat.txt" % cisRoot
48 chromoPath = "%s/download/dr3" % cisRoot
49 chromoOutPath = "/D_rerio/"
50 print "Creating database %s" % db
53 #print "Adding gene entries"
54 #loadGeneEntries(db, genePath)
56 #print "Adding gene features"
57 #loadGeneFeatures(db, genePath)
59 print "Loading chromosomes"
60 loadChromosome(db, chromoPath, chromoOutPath)
62 print "Creating Indices"
65 print "Finished creating database %s" % db
69 drGenome = Genome("drerio", dbFile=db)
70 drGenome.createGeneDB(db)
73 def loadGeneEntries(db, gFile):
76 drGenome = Genome("drerio", dbFile=db)
77 geneFile = open(gFile, "r")
80 cols = line.split("\t")
83 tempID = idb.getGeneID("drerio", gid)
85 print "could not find %s" % gid
88 geneInfo = idb.getGeneInfo(tempID)
105 geneID = ("drerio", gid)
106 if geneID in seenGIDs:
110 seenGIDs.append(geneID)
112 geneEntries.append((geneID, chrom, start, stop, sense, "gene", gidVersion))
114 print "Adding %d gene entries" % len(geneEntries)
115 drGenome.addGeneEntryBatch(geneEntries)
118 def loadGeneFeatures(db, gfile):
119 geneFile = open(gfile, "r")
122 senseArray = {"+": "F",
128 for geneLine in geneFile:
129 geneFields = geneLine.split("\t")
130 exonNum = int(geneFields[8])
131 exonStarts = geneFields[9].split(",")
132 exonStops = geneFields[10].split(",")
133 chrom = geneFields[2]
134 sense = senseArray[geneFields[3]]
135 gstart = int(geneFields[6]) - 1
136 gstop = int(geneFields[7]) - 1
139 tempID = idb.getGeneID("drerio", gid)
141 print "could not find %s" % gid
144 geneInfo = idb.getGeneInfo(tempID)
152 geneID = ("drerio", gid)
153 if geneID in seenGIDs:
158 for index in range(exonNum):
159 estart = int(exonStarts[index]) - 1
160 estop = int(exonStops[index]) - 1
161 if estart >= gstart and estop <= gstop:
162 insertArray.append((geneID, gidVersion, chrom, estart, estop, sense, "CDS"))
163 elif estop <= gstart:
169 insertArray.append((geneID, gidVersion, chrom, estart, estop, sense, fType))
170 elif estart >= gstop:
176 insertArray.append((geneID, gidVersion, chrom, estart, estop, sense, fType))
177 elif estart <= gstop and estart > gstart:
183 insertArray.append((geneID, gidVersion, chrom, estart, gstop, sense, "CDS"))
184 insertArray.append((geneID, gidVersion, chrom, gstop + 1, estop, sense, fType))
185 elif estart < gstart and estop <= gstop:
191 insertArray.append((geneID, gidVersion, chrom, estart, gstart - 1, sense, fType))
192 insertArray.append((geneID, gidVersion, chrom, gstart, estop, sense, "CDS"))
201 insertArray.append((geneID, gidVersion, chrom, estart, gstart - 1, sense, fType1))
202 insertArray.append((geneID, gidVersion, chrom, gstart, gstop, sense, "CDS"))
203 insertArray.append((geneID, gidVersion, chrom, gstop + 1, estop - 1, sense, fType2))
206 drGenome = Genome("drerio", dbFile=db)
207 print "Adding %d features" % len(insertArray)
208 drGenome.addFeatureEntryBatch(insertArray)
211 def loadChromosome(db, chromPath, chromOutPath):
214 drGenome = Genome("drerio", dbFile=db)
215 files = os.listdir(chromPath)
216 for filename in files:
217 inFile = open("%s/%s" % (chromPath, filename), "r")
218 header = inFile.readline()
222 chromID = header.strip()[1:]
223 currentLine = inFile.readline()
224 while currentLine != "" and currentLine[0] != ">":
225 lineSeq = currentLine.strip()
226 seqLen += len(lineSeq)
227 seqArray.append(lineSeq)
228 currentLine = inFile.readline()
230 seq = string.join(seqArray, "")
232 print "Added contig %s to database" % chromID
233 drGenome.addSequence(("drerio", chromID), seq, "chromosome", str(seqLen))
234 drGenome.addChromosomeEntry(chromID, chromID, "db")
236 outFileName = "%s%s.bin" % (chromOutPath, chromID)
237 outFile = open("%s%s" % (cisRoot, outFileName), "w")
240 print "Added contig file %s to database" % outFileName
241 drGenome.addChromosomeEntry(chromID, outFileName, "file")
248 def createDBindices(db):
249 drGenome = Genome("drerio", dbFile=db)
250 drGenome.createIndices()