X-Git-Url: http://woldlab.caltech.edu/gitweb/?p=erange.git;a=blobdiff_plain;f=cistematic%2Fgenomes%2Fcelegans.py;fp=cistematic%2Fgenomes%2Fcelegans.py;h=e3df297269ddba790fc7380c6010cce2cf876f3b;hp=0000000000000000000000000000000000000000;hb=bc30aca13e5ec397c92e67002fbf7a103130b828;hpb=0d3e3112fd04c2e6b44a25cacef1d591658ad181 diff --git a/cistematic/genomes/celegans.py b/cistematic/genomes/celegans.py new file mode 100644 index 0000000..e3df297 --- /dev/null +++ b/cistematic/genomes/celegans.py @@ -0,0 +1,313 @@ +########################################################################### +# # +# 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 Caenorhaditis elegans +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/C_elegans/celegans.genedb" % cisRoot + + +def loadChromosome(db, chromID, chromPath, chromOut): + seqArray = [] + ceGenome = Genome("celegans", 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() + ceGenome.addChromosomeEntry(chromID, chromOut, "file") + + +def loadGeneEntries(db, gffFile): + ceGenome = Genome("celegans", dbFile=db) + geneFile = open(gffFile, "r") + geneEntries = [] + for line in geneFile: + if line[0] == "#": + continue + + field = line.split("\t") + if field[1] != "Coding_transcript" and field[1] != "miRNA": + continue + + if field[2] != "Transcript" and field[2] != "miRNA_primary_transcript": + continue + + gidrev = field[8].split('"') + giddots = gidrev[1].split(".") + # we are ignoring gene models!!!! + if giddots[1][-1] in string.letters: + gidGene = giddots[1][:-1] + gidLetter = giddots[1][-1] + else: + gidGene = giddots[1] + gidLetter = "a" + gid = "%s.%s" % (giddots[0], gidGene) + geneID = ("celegans", gid) + gidVersion = 1 + if gidLetter != "a": + try: + gidVersion = ord(gidLetter.lower()) - 96 + except: + print "problem processing %s - skipping" % gidrev[1] + continue + + start = int(field[3]) - 1 + stop = int(field[4]) - 1 + sense = field[6] + chrom = field[0].strip() + if sense == "+": + sense = "F" + else: + sense = "R" + + geneEntries.append((geneID, chrom, start, stop, sense, "Transcript", gidVersion)) + + print "Adding %d gene entries" % len(geneEntries) + ceGenome.addGeneEntryBatch(geneEntries) + + +def loadFeatureEntries(db, gffFile): + ceGenome = Genome("celegans", dbFile=db) + featureFile = open(gffFile, "r") + featureEntries = [] + seenFeatures = {} + featureTranslation = {"coding_exon": "CDS", + "three_prime_UTR": "3UTR", + "five_prime_UTR": "5UTR" + } + + for line in featureFile: + if line[0] == "#": + continue + + field = line.split("\t") + if field[1] not in ["Coding_transcript", "miRNA", "tRNAscan-SE-1.23"]: + continue + + if field[1] == "Coding_transcript" and field[2].strip() not in featureTranslation: + continue + + if field[1] in ["miRNA", "tRNAscan-SE-1.23"]: + featureType = "CDS" # identifying these as CDS will force their masking later on + else: + featureType = featureTranslation[field[2].strip()] + + gidrev = field[8].split('"') + giddots = gidrev[1].split(".") + # we are ignoring gene models!!!! + if giddots[1][-1] in string.letters: + gidGene = giddots[1][:-1] + gidLetter = giddots[1][-1] + else: + gidGene = giddots[1] + gidLetter = "a" + + gid = "%s.%s" % (giddots[0], gidGene) + geneID = ("celegans", gid) + gidVersion = 1 + if gidLetter != "a": + try: + gidVersion = ord(gidLetter.lower()) - 96 + except: + print "problem processing %s - skipping" % gidrev[1] + continue + + start = int(field[3]) - 1 + stop = int(field[4]) - 1 + sense = field[6] + chrom = field[0].strip() + if sense == "+": + sense = "F" + else: + sense = "R" + + if geneID not in seenFeatures: + seenFeatures[geneID] = [] + + if (gidVersion, start, stop, featureType) not in seenFeatures[geneID]: + featureEntries.append((geneID, gidVersion, chrom, start, stop, sense, featureType)) + seenFeatures[geneID].append((gidVersion, start, stop, featureType)) + + print "Adding %d feature entries" % len(featureEntries) + ceGenome.addFeatureEntryBatch(featureEntries) + + +def loadGeneAnnotations(db, geneIDPath): + geneAnnotations = [] + geneIDFile = open(geneIDPath, "r") + lines = geneIDFile.readlines() + geneIDFile.close() + ceGenome = Genome("celegans", dbFile=db) + for line in lines: + field = line.split(",") + try: + gid = field[2].strip() + geneID = "%s\t%s" % (field[0], field[1]) + geneAnnotations.append((("celegans", gid), geneID)) + except: + pass + + print "Adding %d annotations" % len(geneAnnotations) + ceGenome.addAnnotationBatch(geneAnnotations) + + +def loadGeneOntology(db, goPath, goDefPath, geneIDPath): + ceGenome = Genome("celegans", version="WS200", dbFile=db) + geneIDFile = open(geneIDPath, "r") + goDefFile = open(goDefPath, "r") + goFile = open(goPath, "r") + lines = geneIDFile.readlines() + geneIDFile.close() + geneIDmap = {} + seenGO = {} + for line in lines: + field = line.split(",") + # ugly C elegans hack - map both fields to gid, since either might be + # used by GO ! + if len(field[2].strip()) > 1: + geneIDmap[field[1]] = field[2].strip() + + 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") + name = fields[2] + if name in geneIDmap: + name = geneIDmap[name] + + if name[-1] == "a": + name = name[:-1] + + GOIDarray = fields[4].split(" ") + GOID = GOIDarray[0] + objType = fields[8] + objName = fields[10].split("|") + gID = name + isNot = fields[3] + if len(objName) > 1: + name = "%s|%s" % (name.strip(), fields[10]) + + try: + GOterm = string.replace(goDefs[GOID][0], "'", "p") + except: + print "could no map %s - using GOID only" % GOID + GOterm = "" + + evidence = fields[9] + if gID not in seenGO: + seenGO[gID] = [] + + if GOID not in seenGO[gID]: + seenGO[gID].append(GOID) + goArray.append((("celegans", gID), GOID, objType, name, isNot, GOterm, evidence, fields[1])) + + print "Adding %d GO entries" % len(goArray) + ceGenome.addGoInfoBatch(goArray) + + +def createDBFile(db): + ceGenome = Genome("celegans", version="WS200", dbFile=db) + ceGenome.createGeneDB(db) + + +def createDBindices(db): + ceGenome = Genome("celegans", version="WS200", dbFile=db) + ceGenome.createIndices() + + +def buildCelegansDB(db=geneDB, downloadRoot=""): + if downloadRoot == "": + downloadRoot = "%s/download/" % cisRoot + + geneIDPath = "%sgeneIDs.WS200" % downloadRoot + goDefPath = "%sGO.terms_and_ids" % downloadRoot + goPath = "%sgene_association.wb" % downloadRoot + + # can be found at ftp://caltech.wormbase.org/pub/schwarz/cisreg/softmasks + chromos = {"I": "%sCHROMOSOME_I_softmasked.dna" % downloadRoot, + "II": "%sCHROMOSOME_II_softmasked.dna" % downloadRoot, + "III": "%sCHROMOSOME_III_softmasked.dna" % downloadRoot, + "IV": "%sCHROMOSOME_IV_softmasked.dna" % downloadRoot, + "V": "%sCHROMOSOME_V_softmasked.dna" % downloadRoot, + "X": "%sCHROMOSOME_X_softmasked.dna" % downloadRoot + } + + # can be found at ftp://ftp.wormbase.org/pub/wormbase/genomes/elegans/genome_feature_tables/GFF2/elegansWS160.gff.gz + gffPath = "%selegansWS200.gff" % downloadRoot + + print "Creating database %s" % db + createDBFile(db) + + print "Adding gene entries" + loadGeneEntries(db, gffPath) + + print "Adding feature entries" + loadFeatureEntries(db, gffPath) + + print "Adding gene annotations" + loadGeneAnnotations(db, geneIDPath) + + print "Adding gene ontology" + loadGeneOntology(db, goPath, goDefPath, geneIDPath) + + for chromID in chromos: + print "Loading chromosome %s" % chromID + loadChromosome(db, chromID, chromos[chromID], "/C_elegans/chr%s.bin" % chromID) + + print "Creating Indices" + createDBindices(db) + + print "Finished creating database %s" % db \ No newline at end of file