########################################################################### # # # 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 Xenopus tropicalis 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/X_tropicalis/xtropicalis.genedb" % cisRoot def loadChromosome(db, chromPath, chromOutPath): seqArray = [] seqLen = 0 xtGenome = Genome("xtropicalis", 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 xtGenome.addSequence(("xtropicalis", chromID), seq, "chromosome", str(seqLen)) xtGenome.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 xtGenome.addChromosomeEntry(chromID, outFileName, "file") header = currentLine inFile.close() def loadGeneEntries(db, gFile): """ FIXME - NEED TO DEAL WITH ALTERNATIVE SPLICING ENTRIES """ geneEntries = [] xtGenome = Genome("xtropicalis", 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 = ("xtropicalis", gid) gidVersion = 1 geneEntries.append((geneID, chrom, start, stop, sense, "gene", gidVersion)) print "Adding %d gene entries" % len(geneEntries) xtGenome.addGeneEntryBatch(geneEntries) def loadGeneAnnotations(db, annotPath): geneAnnotations = [] annotFile = open(annotPath, "r") xtGenome = Genome("xtropicalis", dbFile=db) for line in annotFile: try: cols = line.split("\t") locID = cols[0] geneDesc = cols[6] if len(locID) > 0: geneAnnotations.append((("xtropicalis", locID), string.replace(geneDesc.strip(), "'", "p"))) except: pass print "Adding %d annotations" % len(geneAnnotations) xtGenome.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 = ("xtropicalis", 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() xtGenome = Genome("xtropicalis", dbFile=db) print "Adding %d features" % len(insertArray) xtGenome.addFeatureEntryBatch(insertArray) def loadGeneOntology(db, goPath, goDefPath, annotPath): xtGenome = Genome("xtropicalis", 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") locID = cols[0] geneName = cols[1] geneDesc = cols[6] mimID = "" if len(locID) > 0: locus[locID] = (geneName, geneDesc, mimID) except: pass for entry in goEntries: try: fields = entry.split("\t") locID = fields[0].strip() (gene_name, gene_desc, mimID) = locus[locID] goArray.append((("xtropicalis", locID), fields[1], "", gene_name, "", string.replace(goDefs[fields[1]][0], "'", "p"), goDefs[fields[1]][1], mimID)) except: pass print "adding %d go entries" % len(goArray) xtGenome.addGoInfoBatch(goArray) def createDBFile(db): xtGenome = Genome("xtropicalis", dbFile=db) xtGenome.createGeneDB(db) def createDBindices(db): xtGenome = Genome("xtropicalis", dbFile=db) xtGenome.createIndices() def buildXtropicalisDB(db=geneDB): genePath = "%s/download/xt1/jgiFilteredModels.txt" % cisRoot chromoPath = "%s/download/xt1/xenTro1.softmask2.fa" % cisRoot chromoOutPath = "/X_tropicalis/" print "Creating database %s" % db createDBFile(db) print "Adding gene entries" loadGeneEntries(db, genePath) print "Adding gene features" loadGeneFeatures(db, genePath) print "Loading sequences" loadChromosome(db, chromoPath, chromoOutPath) print "Creating Indices" createDBindices(db) print "Finished creating database %s" % db