X-Git-Url: http://woldlab.caltech.edu/gitweb/?p=erange.git;a=blobdiff_plain;f=cistematic%2Fgenomes%2Fxtropicalis.py;fp=cistematic%2Fgenomes%2Fxtropicalis.py;h=eb37500e68a5cfb6f54cf8d054ad47397d2c312b;hp=0000000000000000000000000000000000000000;hb=bc30aca13e5ec397c92e67002fbf7a103130b828;hpb=0d3e3112fd04c2e6b44a25cacef1d591658ad181 diff --git a/cistematic/genomes/xtropicalis.py b/cistematic/genomes/xtropicalis.py new file mode 100644 index 0000000..eb37500 --- /dev/null +++ b/cistematic/genomes/xtropicalis.py @@ -0,0 +1,266 @@ +########################################################################### +# # +# 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 Xenopus tropicalis +import string +from cistematic.genomes import Genome +from os import environ + +if environ.get("CISTEMATIC_ROOT"): + cisRoot = environ.get("CISTEMATIC_ROOT") +else: + cisRoot = "/proj/genome" + +geneDB = "%s/X_tropicalis/xtropicalis.genedb" % cisRoot + + +def loadChromosome(db, chromPath, chromOutPath): + seqArray = [] + seqLen = 0 + xtGenome = Genome("xtropicalis", dbFile=db) + inFile = open(chromPath, "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 < 500000: + print "Added contig %s to database" % chromID + xtGenome.addSequence(("xtropicalis", chromID), seq, "chromosome", str(seqLen)) + xtGenome.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 + xtGenome.addChromosomeEntry(chromID, outFileName, "file") + + header = currentLine + + inFile.close() + + +def loadGeneEntries(db, gFile): + """ FIXME - NEED TO DEAL WITH ALTERNATIVE SPLICING ENTRIES + """ + geneEntries = [] + xtGenome = Genome("xtropicalis", dbFile=db) + geneFile = open(gFile, "r") + for line in geneFile: + cols = line.split("\t") + gid = cols[0] + start = int(cols[5]) + stop = int(cols[6]) + sense = cols[2] + chrom = cols[1] + if sense == "+": + sense = "F" + else: + sense = "R" + + geneID = ("xtropicalis", gid) + gidVersion = 1 + geneEntries.append((geneID, chrom, start, stop, sense, "gene", gidVersion)) + + print "Adding %d gene entries" % len(geneEntries) + xtGenome.addGeneEntryBatch(geneEntries) + + +def loadGeneAnnotations(db, annotPath): + geneAnnotations = [] + annotFile = open(annotPath, "r") + xtGenome = Genome("xtropicalis", dbFile=db) + for line in annotFile: + try: + cols = line.split("\t") + locID = cols[0] + geneDesc = cols[6] + if len(locID) > 0: + geneAnnotations.append((("xtropicalis", locID), string.replace(geneDesc.strip(), "'", "p"))) + except: + pass + + print "Adding %d annotations" % len(geneAnnotations) + xtGenome.addAnnotationBatch(geneAnnotations) + + +def loadGeneFeatures(db, gfile): + geneFile = open(gfile, "r") + senseArray = {"+": "F", + "-": "R", + ".": "F" + } + + seenArray = [] + insertArray = [] + for geneLine in geneFile: + geneFields = geneLine.split("\t") + exonNum = int(geneFields[7]) + exonStarts = geneFields[8].split(",") + exonStops = geneFields[9].split(",") + chrom = geneFields[1] + sense = senseArray[geneFields[2]] + gstop = int(geneFields[6]) - 1 + gstart = int(geneFields[5]) - 1 + geneid = geneFields[0] + try: + geneID = ("xtropicalis", geneid) + except: + continue + + gidVersion = "1" + if geneID in seenArray: + gidVersion = "2" # doesn't deal with more than 2 refseq's for the same locus, yet. + else: + seenArray.append(geneID) + + 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, 1, chrom, estart, estop, sense, fType)) + elif estart >= gstop: + if sense == "F": + fType = "3UTR" + else: + fType = "5UTR" + + insertArray.append((geneID, 1, chrom, estart, estop, sense, fType)) + elif estart <= gstop: + if sense == "F": + fType = "3UTR" + else: + fType = "5UTR" + + insertArray.append((geneID, 1, chrom, estart, gstop, sense, "CDS")) + insertArray.append((geneID, 1, chrom, gstop + 1, estop, sense, fType)) + else: + if sense == "F": + fType = "5UTR" + else: + fType = "3UTR" + + insertArray.append((geneID, 1, chrom, gstart, estop, sense, "CDS")) + insertArray.append((geneID, 1, chrom, estart, gstart - 1, sense, fType)) + + geneFile.close() + xtGenome = Genome("xtropicalis", dbFile=db) + print "Adding %d features" % len(insertArray) + xtGenome.addFeatureEntryBatch(insertArray) + + +def loadGeneOntology(db, goPath, goDefPath, annotPath): + xtGenome = Genome("xtropicalis", dbFile=db) + goDefFile = open(goDefPath, "r") + goFile = open(goPath, "r") + annotFile = open(annotPath, "r") + annotEntries = annotFile.readlines() + annotFile.close() + goDefEntries = goDefFile.readlines() + goDefs = {} + locus = {} + goArray = [] + for goDefEntry in goDefEntries: + if goDefEntry[0] != "!": + cols = goDefEntry.split("\t") + goDefs[cols[0]] = (cols[1], cols[2].strip()) + + goEntries = goFile.readlines() + for annotEntry in annotEntries: + try: + cols = annotEntry.split("\t") + locID = cols[0] + geneName = cols[1] + geneDesc = cols[6] + mimID = "" + if len(locID) > 0: + locus[locID] = (geneName, geneDesc, mimID) + except: + pass + + for entry in goEntries: + try: + fields = entry.split("\t") + locID = fields[0].strip() + (gene_name, gene_desc, mimID) = locus[locID] + goArray.append((("xtropicalis", locID), fields[1], "", gene_name, "", string.replace(goDefs[fields[1]][0], "'", "p"), goDefs[fields[1]][1], mimID)) + except: + pass + + print "adding %d go entries" % len(goArray) + xtGenome.addGoInfoBatch(goArray) + + +def createDBFile(db): + xtGenome = Genome("xtropicalis", dbFile=db) + xtGenome.createGeneDB(db) + + +def createDBindices(db): + xtGenome = Genome("xtropicalis", dbFile=db) + xtGenome.createIndices() + + +def buildXtropicalisDB(db=geneDB): + genePath = "%s/download/xt1/jgiFilteredModels.txt" % cisRoot + chromoPath = "%s/download/xt1/xenTro1.softmask2.fa" % cisRoot + chromoOutPath = "/X_tropicalis/" + + print "Creating database %s" % db + createDBFile(db) + + print "Adding gene entries" + loadGeneEntries(db, genePath) + + print "Adding gene features" + loadGeneFeatures(db, genePath) + + print "Loading sequences" + loadChromosome(db, chromoPath, chromoOutPath) + + print "Creating Indices" + createDBindices(db) + + print "Finished creating database %s" % db \ No newline at end of file