X-Git-Url: http://woldlab.caltech.edu/gitweb/?p=erange.git;a=blobdiff_plain;f=cistematic%2Fgenomes%2Fdrerio.py;fp=cistematic%2Fgenomes%2Fdrerio.py;h=7541db8d9043e2f4a082c73f6f938eb4ff639d9c;hp=0000000000000000000000000000000000000000;hb=bc30aca13e5ec397c92e67002fbf7a103130b828;hpb=0d3e3112fd04c2e6b44a25cacef1d591658ad181 diff --git a/cistematic/genomes/drerio.py b/cistematic/genomes/drerio.py new file mode 100644 index 0000000..7541db8 --- /dev/null +++ b/cistematic/genomes/drerio.py @@ -0,0 +1,250 @@ +########################################################################### +# # +# 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