X-Git-Url: http://woldlab.caltech.edu/gitweb/?p=erange.git;a=blobdiff_plain;f=cistematic%2Fgenomes%2Fecaballus.py;fp=cistematic%2Fgenomes%2Fecaballus.py;h=e845c426500550ef1169df4a3d4bf33fc99690d8;hp=0000000000000000000000000000000000000000;hb=bc30aca13e5ec397c92e67002fbf7a103130b828;hpb=0d3e3112fd04c2e6b44a25cacef1d591658ad181 diff --git a/cistematic/genomes/ecaballus.py b/cistematic/genomes/ecaballus.py new file mode 100644 index 0000000..e845c42 --- /dev/null +++ b/cistematic/genomes/ecaballus.py @@ -0,0 +1,198 @@ +########################################################################### +# # +# 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 Equus Caballus +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/E_caballus/ecaballus.genedb" % cisRoot + + +def loadChromosome(db, chromID, chromPath, chromOut): + seqArray = [] + ecGenome = Genome("ecaballus", 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() + ecGenome.addChromosomeEntry(chromID, chromOut, "file") + + +def loadGeneEntries(db, gFile): + """ FIXME - NEED TO DEAL WITH ALTERNATIVE SPLICING ENTRIES + """ + geneEntries = [] + alreadySeen = [] + ecGenome = Genome("ecaballus", dbFile=db) + geneFile = open(gFile, "r") + geneFile.readline() + for line in geneFile: + cols = line.split("\t") + if cols[11].strip() != "GENE": + continue + + name = cols[10].split(":") + gid = name[1] + if gid == "" or gid in alreadySeen: + continue + + alreadySeen.append(gid) + start = int(cols[2]) - 1 + stop = int(cols[3]) - 1 + sense = cols[4] + chrom = cols[1].strip() + if sense == "+": + sense = "F" + else: + sense = "R" + + geneID = ("ecaballus", gid) + gidVersion = 1 + geneEntries.append((geneID, chrom, start, stop, sense, "gene", gidVersion)) + + print "Adding %d gene entries" % len(geneEntries) + ecGenome.addGeneEntryBatch(geneEntries) + + +def loadGeneFeatures(db, gFile): + """ Load gene features such as CDS, UTR, and PSEUDO from the gene file. + """ + featureEntries = [] + ecGenome = Genome("ecaballus", dbFile=db) + featureFile = open(gFile, "r") + featureFile.readline() + for line in featureFile: + cols = line.split("\t") + if cols[11].strip() not in ["CDS", "UTR", "PSEUDO"]: + continue + + fType = cols[11] + name = cols[10].split(":") + gid = name[1] + if gid == "": + continue + + start = int(cols[2]) - 1 + stop = int(cols[3]) - 1 + sense = cols[4] + chrom = cols[1].strip() + if sense == "+": + sense = "F" + else: + sense = "R" + + geneID = ("ecaballus", gid) + gidVersion = 1 + featureEntries.append((geneID, gidVersion, chrom, start, stop, sense, fType)) + + print "Adding %d feature entries" % len(featureEntries) + ecGenome.addFeatureEntryBatch(featureEntries) + + +def createDBFile(db): + ecGenome = Genome("ecaballus", dbFile=db) + ecGenome.createGeneDB(db) + + +def createDBindices(db): + ecGenome = Genome("ecaballus", dbFile=db) + ecGenome.createIndices() + + +def buildHorseDB(db=geneDB): + genePath = "%s/download/seq_gene.md" % cisRoot + chromos = {"1": "%s/download/chr1.fa" % cisRoot, + "2": "%s/download/chr2.fa" % cisRoot, + "3": "%s/download/chr3.fa" % cisRoot, + "4": "%s/download/chr4.fa" % cisRoot, + "5": "%s/download/chr5.fa" % cisRoot, + "6": "%s/download/chr6.fa" % cisRoot, + "7": "%s/download/chr7.fa" % cisRoot, + "8": "%s/download/chr8.fa" % cisRoot, + "9": "%s/download/chr9.fa" % cisRoot, + "10": "%s/download/chr10.fa" % cisRoot, + "11": "%s/download/chr11.fa" % cisRoot, + "12": "%s/download/chr12.fa" % cisRoot, + "13": "%s/download/chr13.fa" % cisRoot, + "14": "%s/download/chr14.fa" % cisRoot, + "15": "%s/download/chr15.fa" % cisRoot, + "16": "%s/download/chr16.fa" % cisRoot, + "17": "%s/download/chr17.fa" % cisRoot, + "18": "%s/download/chr18.fa" % cisRoot, + "19": "%s/download/chr19.fa" % cisRoot, + "20": "%s/download/chr20.fa" % cisRoot, + "21": "%s/download/chr21.fa" % cisRoot, + "22": "%s/download/chr22.fa" % cisRoot, + "23": "%s/download/chr23.fa" % cisRoot, + "24": "%s/download/chr24.fa" % cisRoot, + "25": "%s/download/chr25.fa" % cisRoot, + "26": "%s/download/chr26.fa" % cisRoot, + "27": "%s/download/chr27.fa" % cisRoot, + "28": "%s/download/chr28.fa" % cisRoot, + "29": "%s/download/chr29.fa" % cisRoot, + "30": "%s/download/chr30.fa" % cisRoot, + "31": "%s/download/chr31.fa" % cisRoot, + "M": "%s/download/chrM.fa" % cisRoot, + "X": "%s/download/chrX.fa" % cisRoot, + "Un": "%s/download/chrUn.fa" % cisRoot + } + + print "Creating database %s" % db + createDBFile(db) + + print "Adding gene entries" + loadGeneEntries(db, genePath) + + print "Adding gene features" + loadGeneFeatures(db, genePath) + + for chromID in chromos.keys(): + print "Loading chromosome %s" % chromID + loadChromosome(db, chromID, chromos[chromID], "/E_caballus/chromo%s.bin" % chromID) + + print "Creating Indices" + createDBindices(db) + + print "Finished creating database %s" % db \ No newline at end of file