--- /dev/null
+###########################################################################
+# #
+# C O P Y R I G H T N O T I C E #
+# Copyright (c) 2003-10 by: #
+# * California Institute of Technology #
+# #
+# All Rights Reserved. #
+# #
+# Permission is hereby granted, free of charge, to any person #
+# obtaining a copy of this software and associated documentation files #
+# (the "Software"), to deal in the Software without restriction, #
+# including without limitation the rights to use, copy, modify, merge, #
+# publish, distribute, sublicense, and/or sell copies of the Software, #
+# and to permit persons to whom the Software is furnished to do so, #
+# subject to the following conditions: #
+# #
+# The above copyright notice and this permission notice shall be #
+# included in all copies or substantial portions of the Software. #
+# #
+# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, #
+# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF #
+# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND #
+# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS #
+# BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN #
+# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN #
+# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE #
+# SOFTWARE. #
+###########################################################################
+#
+# data for Danio Rerio
+import string, os
+from cistematic.genomes import Genome
+from cistematic.core.geneinfo import geneinfoDB
+from os import environ
+
+if environ.get("CISTEMATIC_ROOT"):
+ cisRoot = environ.get("CISTEMATIC_ROOT")
+else:
+ cisRoot = "/proj/genome"
+
+geneDB = "%s/D_rerio/drerio.genedb" % cisRoot
+
+
+def loadChromosome(db, chromPath, chromOutPath):
+ seqArray = []
+ seqLen = 0
+ drGenome = Genome("drerio", dbFile=db)
+ files = os.listdir(chromPath)
+ for filename in files:
+ inFile = open("%s/%s" % (chromPath, filename), "r")
+ header = inFile.readline()
+ while header != "":
+ seqArray = []
+ seqLen = 0
+ chromID = header.strip()[1:]
+ currentLine = inFile.readline()
+ while currentLine != "" and currentLine[0] != ">":
+ lineSeq = currentLine.strip()
+ seqLen += len(lineSeq)
+ seqArray.append(lineSeq)
+ currentLine = inFile.readline()
+
+ seq = string.join(seqArray, "")
+ if seqLen < 250000:
+ print "Added contig %s to database" % chromID
+ drGenome.addSequence(("drerio", chromID), seq, "chromosome", str(seqLen))
+ drGenome.addChromosomeEntry(chromID, chromID, "db")
+ else:
+ outFileName = "%s%s.bin" % (chromOutPath, chromID)
+ outFile = open("%s%s" % (cisRoot, outFileName), "w")
+ outFile.write(seq)
+ outFile.close()
+ print "Added contig file %s to database" % outFileName
+ drGenome.addChromosomeEntry(chromID, outFileName, "file")
+
+ header = currentLine
+
+ inFile.close()
+
+
+def loadGeneEntries(db, gFile):
+ geneEntries = []
+ seenGIDs = []
+ drGenome = Genome("drerio", dbFile=db)
+ geneFile = open(gFile, "r")
+ idb = geneinfoDB()
+ for line in geneFile:
+ cols = line.split("\t")
+ gid = cols[0]
+ try:
+ tempID = idb.getGeneID("drerio", gid)
+ if not len(tempID):
+ print "could not find %s" % gid
+ continue
+
+ geneInfo = idb.getGeneInfo(tempID)
+ gid = geneInfo[1]
+ except:
+ continue
+
+ if gid == "":
+ continue
+
+ start = int(cols[6])
+ stop = int(cols[7])
+ sense = cols[3]
+ chrom = cols[2]
+ if sense == "-":
+ sense = "R"
+ else:
+ sense = "F"
+
+ geneID = ("drerio", gid)
+ if geneID in seenGIDs:
+ gidVersion = "2"
+ else:
+ gidVersion = "1"
+ seenGIDs.append(geneID)
+
+ geneEntries.append((geneID, chrom, start, stop, sense, "gene", gidVersion))
+
+ print "Adding %d gene entries" % len(geneEntries)
+ drGenome.addGeneEntryBatch(geneEntries)
+
+
+def loadGeneFeatures(db, gfile):
+ geneFile = open(gfile, "r")
+ idb = geneinfoDB()
+ seenGIDs = []
+ senseArray = {"+": "F",
+ "-": "R",
+ ".": "F"
+ }
+
+ insertArray = []
+ for geneLine in geneFile:
+ geneFields = geneLine.split("\t")
+ exonNum = int(geneFields[8])
+ exonStarts = geneFields[9].split(",")
+ exonStops = geneFields[10].split(",")
+ chrom = geneFields[2]
+ sense = senseArray[geneFields[3]]
+ gstart = int(geneFields[6]) - 1
+ gstop = int(geneFields[7]) - 1
+ gid = geneFields[0]
+ try:
+ tempID = idb.getGeneID("drerio", gid)
+ if not len(tempID):
+ print "could not find %s" % gid
+ continue
+
+ geneInfo = idb.getGeneInfo(tempID)
+ gid = geneInfo[1]
+ except:
+ continue
+
+ if gid == "":
+ continue
+
+ geneID = ("drerio", gid)
+ if geneID in seenGIDs:
+ gidVersion = "2"
+ else:
+ gidVersion = "1"
+
+ for index in range(exonNum):
+ estart = int(exonStarts[index]) - 1
+ estop = int(exonStops[index]) - 1
+ if estart >= gstart and estop <= gstop:
+ insertArray.append((geneID, gidVersion, chrom, estart, estop, sense, "CDS"))
+ elif estop <= gstart:
+ if sense == "F":
+ fType = "5UTR"
+ else:
+ fType = "3UTR"
+
+ insertArray.append((geneID, gidVersion, chrom, estart, estop, sense, fType))
+ elif estart >= gstop:
+ if sense == "F":
+ fType = "3UTR"
+ else:
+ fType = "5UTR"
+
+ insertArray.append((geneID, gidVersion, chrom, estart, estop, sense, fType))
+ elif estart <= gstop and estart > gstart:
+ if sense == "F":
+ fType = "3UTR"
+ else:
+ fType = "5UTR"
+
+ insertArray.append((geneID, gidVersion, chrom, estart, gstop, sense, "CDS"))
+ insertArray.append((geneID, gidVersion, chrom, gstop + 1, estop, sense, fType))
+ elif estart < gstart and estop <= gstop:
+ if sense == "F":
+ fType = "5UTR"
+ else:
+ fType = "3UTR"
+
+ insertArray.append((geneID, gidVersion, chrom, estart, gstart - 1, sense, fType))
+ insertArray.append((geneID, gidVersion, chrom, gstart, estop, sense, "CDS"))
+ else:
+ if sense == "F":
+ fType1 = "5UTR"
+ fType2 = "3UTR"
+ else:
+ fType1 = "3UTR"
+ fType2 = "5UTR"
+
+ insertArray.append((geneID, gidVersion, chrom, estart, gstart - 1, sense, fType1))
+ insertArray.append((geneID, gidVersion, chrom, gstart, gstop, sense, "CDS"))
+ insertArray.append((geneID, gidVersion, chrom, gstop + 1, estop - 1, sense, fType2))
+
+ geneFile.close()
+ drGenome = Genome("drerio", dbFile=db)
+ print "Adding %d features" % len(insertArray)
+ drGenome.addFeatureEntryBatch(insertArray)
+
+
+def createDBFile(db):
+ drGenome = Genome("drerio", dbFile=db)
+ drGenome.createGeneDB(db)
+
+
+def createDBindices(db):
+ drGenome = Genome("drerio", dbFile=db)
+ drGenome.createIndices()
+
+
+def buildDrerioDB(db=geneDB):
+ """ genes and annotations are from UCSC (dr3).
+ """
+ #genePath = "%s/download/xenoRefFlat.txt" % cisRoot
+ chromoPath = "%s/download/dr3" % cisRoot
+ chromoOutPath = "/D_rerio/"
+ print "Creating database %s" % db
+ createDBFile(db)
+
+ #print "Adding gene entries"
+ #loadGeneEntries(db, genePath)
+
+ #print "Adding gene features"
+ #loadGeneFeatures(db, genePath)
+
+ print "Loading chromosomes"
+ loadChromosome(db, chromoPath, chromoOutPath)
+
+ print "Creating Indices"
+ createDBindices(db)
+
+ print "Finished creating database %s" % db
\ No newline at end of file
projects / erange.git / blobdiff