X-Git-Url: http://woldlab.caltech.edu/gitweb/?p=erange.git;a=blobdiff_plain;f=cistematic%2Fgenomes%2Fdmelanogaster.py;fp=cistematic%2Fgenomes%2Fdmelanogaster.py;h=5c9f53a58bbc4466689f9ee2dc40378db3367c10;hp=0000000000000000000000000000000000000000;hb=bc30aca13e5ec397c92e67002fbf7a103130b828;hpb=0d3e3112fd04c2e6b44a25cacef1d591658ad181 diff --git a/cistematic/genomes/dmelanogaster.py b/cistematic/genomes/dmelanogaster.py new file mode 100644 index 0000000..5c9f53a --- /dev/null +++ b/cistematic/genomes/dmelanogaster.py @@ -0,0 +1,330 @@ +########################################################################### +# # +# 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 Drosophila melanogaster +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/D_melanogaster/dmelanogaster.genedb" % cisRoot + +version = {"A": "1", + "B": "2", + "C": "3", + "D": "4", + "E": "5", + "F": "6", + "G": "7", + "H": "8", + "I": "9", + "J": "10", + "K": "11", + "L": "12", + "M": "13", + "N": "14", + "O": "15", + "P": "16", + "Q": "17", + "R": "18", + "S": "19", + "T": "20", + "U": "21", + "V": "22", + "W": "23", + "X": "24", + "Y": "25", + "Z": "26" +} + +def loadChromosome(db, chromID, chromPath, chromOut): + seqArray = [] + dmGenome = Genome("dmelanogaster", 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() + dmGenome.addChromosomeEntry(chromID, chromOut, "file") + + +def loadGeneEntries(db, gFile): + geneEntries = [] + dmGenome = Genome("dmelanogaster", dbFile=db) + geneFile = open(gFile, "r") + + for line in geneFile: + cols = line.split("\t") + name = cols[1].split("-R") + gid = name[0] + start = int(cols[4]) + stop = int(cols[5]) + sense = cols[3] + chrom = cols[2][3:] + if sense == "-": + sense = "R" + else: + sense = "F" + + geneID = ("dmelanogaster", gid) + try: + gidVersion = version[name[1]] + except: + continue + + geneEntries.append((geneID, chrom, start, stop, sense, "gene", gidVersion)) + + print "Adding %d gene entries" % len(geneEntries) + dmGenome.addGeneEntryBatch(geneEntries) + + +def loadGeneFeatures(db, gfile): + geneFile = open(gfile, "r") + 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][3:] + sense = senseArray[geneFields[3]] + gstop = int(geneFields[7]) - 1 + gstart = int(geneFields[6]) - 1 + name = geneFields[1].split("-R") + geneID = ("dmelanogaster", name[0]) + try: + gidVersion = version[name[1]] + except: + continue + + 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() + dmGenome = Genome("dmelanogaster", dbFile=db) + print "Adding %d features" % len(insertArray) + dmGenome.addFeatureEntryBatch(insertArray) + + +def loadGeneAnnotations(db, annotPath): + geneAnnotations = [] + annotFile = open(annotPath, "r") + dmGenome = Genome("dmelanogaster", dbFile=db) + for line in annotFile: + try: + cols = line.split("\t") + if cols[0] != "7227": + continue + + locID = cols[3] + if "Dmel_" in locID: + locID = locID[5:] + + geneDesc = cols[4] + if len(locID) > 0: + geneAnnotations.append((("dmelanogaster", locID), string.replace(geneDesc.strip(), "'", "p"))) + except: + pass + + print "Adding %d annotations" % len(geneAnnotations) + dmGenome.addAnnotationBatch(geneAnnotations) + + +def loadGeneOntology(db, goPath, goDefPath, annotPath): + dmGenome = Genome("dmelanogaster", 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") + if cols[0] != "7227": + continue + + locID = cols[3].strip() + geneName = cols[1].strip() + if len(locID) > 0: + locus[geneName] = locID + except: + pass + + for entry in goEntries: + if entry[0] == "!": + continue + + if entry[:4] != "7227": + continue + + try: + fields = entry.split("\t") + geneName = fields[1].strip() + locID = locus[geneName] + if "Dmel_" in locID: + locID = locID[5:] + + GOID = fields[2] + goArray.append((("dmelanogaster", locID), GOID, "", geneName, "", string.replace(goDefs[GOID][0], "'", "p"), goDefs[GOID][1], "")) + except: + pass + + print "adding %d go entries" % len(goArray) + dmGenome.addGoInfoBatch(goArray) + + +def createDBFile(db): + dmGenome = Genome("dmelanogaster", dbFile=db) + dmGenome.createGeneDB(db) + + +def createDBindices(db): + dmGenome = Genome("dmelanogaster", dbFile=db) + dmGenome.createIndices() + + +def buildDmelanogasterDB(db=geneDB): + """ genes and annotations are from UCSC. GO association file is from geneontology.org. + """ + genePath = "%s/download/flyBaseGene.txt" % cisRoot + annotPath = "%s/download/gene_info" % cisRoot + goDefPath = "%s/download/GO.terms_and_ids" % cisRoot + goPath = "%s/download/gene2go" % cisRoot + chromos = {"2L": "%s/download/chr2L.fa" % cisRoot, + "2R": "%s/download/chr2R.fa" % cisRoot, + "2LHet": "%s/download/chr2LHet.fa" % cisRoot, + "2RHet": "%s/download/chr2RHet.fa" % cisRoot, + "3L": "%s/download/chr3L.fa" % cisRoot, + "3LHet": "%s/download/chr3LHet.fa" % cisRoot, + "3R": "%s/download/chr3R.fa" % cisRoot, + "3RHet": "%s/download/chr3RHet.fa" % cisRoot, + "4": "%s/download/chr4.fa" % cisRoot, + "X": "%s/download/chrX.fa" % cisRoot, + "XHet": "%s/download/chrXHet.fa" % cisRoot, + "YHet": "%s/download/chrYHet.fa" % cisRoot, + "U": "%s/download/chrU.fa" % cisRoot, + "Uextra": "%s/download/chrUextra.fa" % cisRoot, + "M": "%s/download/chrM.fa" % cisRoot + } + + print "Creating database %s" % db + createDBFile(db) + + print "Adding gene entries" + loadGeneEntries(db, genePath) + + print "Adding gene features" + loadGeneFeatures(db, genePath) + + print "Adding gene annotations" + loadGeneAnnotations(db, annotPath) + + print "Adding gene ontology" + loadGeneOntology(db, goPath, goDefPath, annotPath) + + for chromID in chromos.keys(): + print "Loading chromosome %s" % chromID + loadChromosome(db, chromID, chromos[chromID], "/D_melanogaster/chromo%s.bin" % chromID) + + print "Creating Indices" + createDBindices(db) + + print "Finished creating database %s" % db \ No newline at end of file