########################################################################### # # # 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 Danio Rerio import string, os from cistematic.genomes import Genome from cistematic.core.geneinfo import geneinfoDB from os import environ if environ.get("CISTEMATIC_ROOT"): cisRoot = environ.get("CISTEMATIC_ROOT") else: cisRoot = "/proj/genome" geneDB = "%s/D_rerio/drerio.genedb" % cisRoot def loadChromosome(db, chromPath, chromOutPath): seqArray = [] seqLen = 0 drGenome = Genome("drerio", dbFile=db) files = os.listdir(chromPath) for filename in files: inFile = open("%s/%s" % (chromPath, filename), "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 < 250000: print "Added contig %s to database" % chromID drGenome.addSequence(("drerio", chromID), seq, "chromosome", str(seqLen)) drGenome.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 drGenome.addChromosomeEntry(chromID, outFileName, "file") header = currentLine inFile.close() def loadGeneEntries(db, gFile): geneEntries = [] seenGIDs = [] drGenome = Genome("drerio", dbFile=db) geneFile = open(gFile, "r") idb = geneinfoDB() for line in geneFile: cols = line.split("\t") gid = cols[0] try: tempID = idb.getGeneID("drerio", gid) if not len(tempID): print "could not find %s" % gid continue geneInfo = idb.getGeneInfo(tempID) gid = geneInfo[1] except: continue if gid == "": continue start = int(cols[6]) stop = int(cols[7]) sense = cols[3] chrom = cols[2] if sense == "-": sense = "R" else: sense = "F" geneID = ("drerio", gid) if geneID in seenGIDs: gidVersion = "2" else: gidVersion = "1" seenGIDs.append(geneID) geneEntries.append((geneID, chrom, start, stop, sense, "gene", gidVersion)) print "Adding %d gene entries" % len(geneEntries) drGenome.addGeneEntryBatch(geneEntries) def loadGeneFeatures(db, gfile): geneFile = open(gfile, "r") idb = geneinfoDB() seenGIDs = [] 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] sense = senseArray[geneFields[3]] gstart = int(geneFields[6]) - 1 gstop = int(geneFields[7]) - 1 gid = geneFields[0] try: tempID = idb.getGeneID("drerio", gid) if not len(tempID): print "could not find %s" % gid continue geneInfo = idb.getGeneInfo(tempID) gid = geneInfo[1] except: continue if gid == "": continue geneID = ("drerio", gid) if geneID in seenGIDs: gidVersion = "2" else: gidVersion = "1" 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() drGenome = Genome("drerio", dbFile=db) print "Adding %d features" % len(insertArray) drGenome.addFeatureEntryBatch(insertArray) def createDBFile(db): drGenome = Genome("drerio", dbFile=db) drGenome.createGeneDB(db) def createDBindices(db): drGenome = Genome("drerio", dbFile=db) drGenome.createIndices() def buildDrerioDB(db=geneDB): """ genes and annotations are from UCSC (dr3). """ #genePath = "%s/download/xenoRefFlat.txt" % cisRoot chromoPath = "%s/download/dr3" % cisRoot chromoOutPath = "/D_rerio/" 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