X-Git-Url: http://woldlab.caltech.edu/gitweb/?p=erange.git;a=blobdiff_plain;f=cistematic%2Fgenomes%2Fscerevisiae.py;fp=cistematic%2Fgenomes%2Fscerevisiae.py;h=5866f80963fac6914c81048fb1460548a4ef8a3b;hp=0000000000000000000000000000000000000000;hb=bc30aca13e5ec397c92e67002fbf7a103130b828;hpb=0d3e3112fd04c2e6b44a25cacef1d591658ad181 diff --git a/cistematic/genomes/scerevisiae.py b/cistematic/genomes/scerevisiae.py new file mode 100644 index 0000000..5866f80 --- /dev/null +++ b/cistematic/genomes/scerevisiae.py @@ -0,0 +1,216 @@ +########################################################################### +# # +# 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 Saccharomyces cerevisiae +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/S_cerevisiae/scerevisiae.genedb" % cisRoot + + +def loadChromosome(db, chromID, chromPath, chromOut): + seqArray = [] + scGenome = Genome("scerevisiae", dbFile=db) + inFile = open(chromPath, "r") + line = inFile.readline() + for line in inFile: + seqArray.append(line.strip()) + + seq = string.join(seqArray, "") + seqLen = len(seq) + if seqLen < 1: + print "Problems reading sequence from file" + + print "writing to file %s" % chromOut + outFile = open("%s%s" % (cisRoot, chromOut), "w") + outFile.write(seq) + outFile.close() + seq = "" + print "calling scGenome()" + scGenome.addChromosomeEntry(chromID, chromOut, "file") + + +def loadGeneEntries(db, gFile): + geneEntries = [] + geneFeatures = [] + scGenome = Genome("scerevisiae", dbFile=db) + geneFile = open(gFile, "r") + for line in geneFile: + field = line.split("\t") + if field[1] != "ORF": + continue + + orfName = field[3].strip() + sense = field[11] + chrom = field[8].strip() + if sense == "W": + sense = "F" + try: + start = int(field[9].strip()) - 1 + stop = int(field[10].strip()) - 1 + except: + start = 0 + stop = 0 + else: + sense = "R" + try: + start = int(field[10].strip()) - 1 + stop = int(field[9].strip()) - 1 + except: + start = 0 + stop = 0 + + geneID = ("scerevisiae", orfName) + gidVersion = 1 + geneEntries.append((geneID, chrom, start, stop, sense, "chromosomal_feature", gidVersion)) + geneFeatures.append((geneID, gidVersion, chrom, start, stop, sense, "CDS")) + + print "loading %d gene entries" % len(geneEntries) + scGenome.addGeneEntryBatch(geneEntries) + print "loading %d gene features" % len(geneFeatures) + scGenome.addFeatureEntryBatch(geneFeatures) + + +def loadGeneAnnotations(db, annotPath): + geneAnnotations = [] + annotFile = open(annotPath, "r") + lines = annotFile.readlines() + annotFile.close() + scGenome = Genome("scerevisiae", dbFile=db) + for line in lines: + field = line.split("\t") + if field[1] != "ORF": + continue + + try: + orfName = field[6].strip() + description = field[15].strip() + geneAnnotations.append((("scerevisiae", orfName), string.replace(description, "'", "p"))) + except: + pass + + print "Adding %d annotations" % len(geneAnnotations) + scGenome.addAnnotationBatch(geneAnnotations) + + +def loadGeneOntology(db, goPath, goDefPath): + scGenome = Genome("scerevisiae", version="SGD1", dbFile=db) + goDefFile = open(goDefPath, "r") + goFile = open(goPath, "r") + goDefEntries = goDefFile.readlines() + goDefs = {} + for goDefEntry in goDefEntries: + if goDefEntry[0] != "!": + cols = goDefEntry.split("\t") + goDefs[cols[0]] = (cols[1], cols[2].strip()) + + goEntries = goFile.readlines() + goArray = [] + for line in goEntries: + if line[0] == "!": + continue + + fields = line.split("\t") + genes = fields[10].split("|") + gID = genes[0] + GOID = fields[4] + objType = fields[11] + objNameArray = fields[10].split("|") + objName = objNameArray[0] + isNot = fields[3] + try: + GOterm = string.replace(goDefs[GOID][0], "'", "p") + except: + print "Could not translate %s" % (GOID) + GOterm = "" + + evidence = fields[6] + goArray.append((("scerevisiae", gID), GOID[3:], objType, objName, isNot, GOterm, evidence, fields[1])) + + scGenome.addGoInfoBatch(goArray) + + +def createDBFile(db): + scGenome = Genome("scerevisiae", version="SGD1", dbFile=db) + scGenome.createGeneDB(db) + + +def createDBindices(db): + scGenome = Genome("scerevisiae", version="SGD1", dbFile=db) + scGenome.createIndices() + + +def buildScerevisiaeDB(db=geneDB): + genePath = "%s/download/SGD_features.tab" % cisRoot + goDefPath = "%s/download/GO.terms_and_ids" % cisRoot + goPath = "%s/download/gene_association.sgd" % cisRoot + chromos = {"1": "%s/download/chr01.fsa" % cisRoot, + "2": "%s/download/chr02.fsa" % cisRoot, + "3": "%s/download/chr03.fsa" % cisRoot, + "4": "%s/download/chr04.fsa" % cisRoot, + "5": "%s/download/chr05.fsa" % cisRoot, + "6": "%s/download/chr06.fsa" % cisRoot, + "7": "%s/download/chr07.fsa" % cisRoot, + "8": "%s/download/chr08.fsa" % cisRoot, + "9": "%s/download/chr09.fsa" % cisRoot, + "10": "%s/download/chr10.fsa" % cisRoot, + "11": "%s/download/chr11.fsa" % cisRoot, + "12": "%s/download/chr12.fsa" % cisRoot, + "13": "%s/download/chr13.fsa" % cisRoot, + "14": "%s/download/chr14.fsa" % cisRoot, + "15": "%s/download/chr15.fsa" % cisRoot, + "16": "%s/download/chr16.fsa" % cisRoot + } + + print "Creating database %s" % db + createDBFile(db) + + print "Adding gene entries" + loadGeneEntries(db, genePath) + + print "Adding gene annotations" + loadGeneAnnotations(db, genePath) + + print "Adding gene ontology" + loadGeneOntology(db, goPath, goDefPath) + + for chromID in ["1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16"]: + print "Loading chromosome %s" % chromID + loadChromosome(db, chromID, chromos[chromID], "/S_cerevisiae/chr%s.bin" % chromID) + + print "Creating Indices" + createDBindices(db) + + print "Finished creating database %s" % db \ No newline at end of file