--- /dev/null
+###########################################################################
+# #
+# 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 Monodelphis domestica
+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/M_domestica/mdomestica.genedb" % cisRoot
+
+
+def loadChromosome(db, chromPath, chromOutPath):
+ seqArray = []
+ seqLen = 0
+ mdGenome = Genome("mdomestica", 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
+ mdGenome.addSequence(("mdomestica", chromID), seq, "chromosome", str(seqLen))
+ mdGenome.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
+ mdGenome.addChromosomeEntry(chromID, outFileName, "file")
+
+ header = currentLine
+
+ inFile.close()
+
+
+def loadGeneEntries(db, gFile):
+ """ FIXME - NEED TO DEAL WITH ALTERNATIVE SPLICING ENTRIES
+ """
+ geneEntries = []
+ mdGenome = Genome("mdomestica", 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 = ("mdomestica", gid)
+ gidVersion = 1
+ geneEntries.append((geneID, chrom, start, stop, sense, "gene", gidVersion))
+
+ print "Adding %d gene entries" % len(geneEntries)
+ mdGenome.addGeneEntryBatch(geneEntries)
+
+
+def loadGeneAnnotations(db, annotPath):
+ geneAnnotations = []
+ annotFile = open(annotPath, "r")
+ mdGenome = Genome("mdomestica", dbFile=db)
+ for line in annotFile:
+ try:
+ cols = line.split("\t")
+ locID = cols[0]
+ geneDesc = cols[6]
+ if len(locID) > 0:
+ geneAnnotations.append((("mdomestica", locID), string.replace(geneDesc.strip(), "'", "p")))
+ except:
+ pass
+
+ print "Adding %d annotations" % len(geneAnnotations)
+ mdGenome.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 = ("mdomestica", 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()
+ mdGenome = Genome("mdomestica", dbFile=db)
+ print "Adding %d features" % len(insertArray)
+ mdGenome.addFeatureEntryBatch(insertArray)
+
+
+def createDBFile(db):
+ mdGenome = Genome("mdomestica", dbFile=db)
+ mdGenome.createGeneDB(db)
+
+
+def createDBindices(db):
+ mdGenome = Genome("mdomestica", dbFile=db)
+ mdGenome.createIndices()
+
+
+def buildMdomesticaDB(db=geneDB):
+ genePath = "%s/download/mondom/genscan.txt" % cisRoot
+ chromoPath = "%s/download/mondom/softMask.fa" % cisRoot
+ chromoOutPath = "/M_domestica/"
+
+ print "Creating database %s" % db
+ createDBFile(db)
+
+ print "Adding gene entries"
+ loadGeneEntries(db, genePath)
+
+ print "Adding gene features"
+ loadGeneFeatures(db, genePath)
+
+ print "Loading chromosomes"
+ loadChromosome(db, chromoPath, chromoOutPath)
+
+ print "Creating Indices"
+ createDBindices(db)
+
+ print "Finished creating database %s" % db
\ No newline at end of file
projects / erange.git / blobdiff