########################################################################### # # # 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 Gallus gallus import string 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/G_gallus/ggallus.genedb" % cisRoot def loadChromosome(db, chromID, chromPath, chromOut): seqArray = [] ggGenome = Genome("ggallus", 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() ggGenome.addChromosomeEntry(chromID, chromOut, "file") def loadGeneEntries(db, gFile): """ FIXME - NEED TO DEAL WITH ALTERNATIVE SPLICING ENTRIES """ geneEntries = [] alreadySeen = [] ggGenome = Genome("ggallus", dbFile=db) geneFile = open(gFile, "r") geneFile.readline() for line in geneFile: if "|" in line: continue 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 = ("ggallus", gid) gidVersion = 1 geneEntries.append((geneID, chrom, start, stop, sense, "gene", gidVersion)) print "Adding %d gene entries" % len(geneEntries) ggGenome.addGeneEntryBatch(geneEntries) def loadGeneFeatures(db, gFile): """ Load gene features such as CDS, UTR, and PSEUDO from the gene file. """ featureEntries = [] ggGenome = Genome("ggallus", dbFile=db) featureFile = open(gFile, "r") featureFile.readline() for line in featureFile: if "|" in line: continue 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 = ("ggallus", gid) gidVersion = 1 featureEntries.append((geneID, gidVersion, chrom, start, stop, sense, fType)) print "Adding %d feature entries" % len(featureEntries) ggGenome.addFeatureEntryBatch(featureEntries) def loadGeneOntology(db, goPath, goDefPath): ggGenome = Genome("ggallus", dbFile=db) goDefFile = open(goDefPath, "r") goFile = open(goPath, "r") idb = geneinfoDB() goDefs = {} goArray = [] for goDefEntry in goDefFile: if goDefEntry[0] != "!": cols = goDefEntry.split("\t") goDefs[cols[0]] = (cols[1], cols[2].strip()) goEntries = goFile.readlines() prevGID = "" for entry in goEntries: try: fields = entry.split("\t") if fields[0] != "9031": continue locID = fields[1].strip() gID = ("ggallus", locID) if prevGID != gID: prevGID = gID gene_name = "" synonyms = idb.geneIDSynonyms(gID) if len(synonyms) >0: gene_name = string.join(synonyms, ",") goArray.append((gID, fields[2], "", gene_name, "", string.replace(goDefs[fields[2]][0], "'", "p"), goDefs[fields[2]][1], "")) except: print "locus ID %s could not be added" % locID pass print "adding %d go entries" % len(goArray) ggGenome.addGoInfoBatch(goArray) def createDBFile(db): ggGenome = Genome("ggallus", dbFile=db) ggGenome.createGeneDB(db) def createDBindices(db): ggGenome = Genome("ggallus", dbFile=db) ggGenome.createIndices() def buildChickenDB(db=geneDB): genePath = "%s/download/seq_gene.md" % cisRoot goDefPath = "%s/download/GO.terms_and_ids" % cisRoot # ftp://ftp.geneontology.org/pub/go/doc/GO.terms_and_ids goPath = "%s/download/gene2go" % cisRoot # ftp://ftp.ncbi.nih.gov/gene/DATA/gene2go.gz 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, "32": "%s/download/chr32.fa" % cisRoot, "W": "%s/download/chrW.fa" % cisRoot, "Z": "%s/download/chrZ.fa" % cisRoot, "M": "%s/download/chrM.fa" % cisRoot, "E22C19W28_E50C23": "%s/download/chrE22C19W28_E50C23.fa" % cisRoot, "E64": "%s/download/chrE64.fa" % cisRoot, "1_random": "%s/download/chr1_random.fa" % cisRoot, "2_random": "%s/download/chr2_random.fa" % cisRoot, "4_random": "%s/download/chr4_random.fa" % cisRoot, "5_random": "%s/download/chr5_random.fa" % cisRoot, "6_random": "%s/download/chr6_random.fa" % cisRoot, "7_random": "%s/download/chr7_random.fa" % cisRoot, "8_random": "%s/download/chr8_random.fa" % cisRoot, "10_random": "%s/download/chr10_random.fa" % cisRoot, "11_random": "%s/download/chr11_random.fa" % cisRoot, "12_random": "%s/download/chr12_random.fa" % cisRoot, "13_random": "%s/download/chr13_random.fa" % cisRoot, "16_random": "%s/download/chr16_random.fa" % cisRoot, "17_random": "%s/download/chr17_random.fa" % cisRoot, "18_random": "%s/download/chr18_random.fa" % cisRoot, "20_random": "%s/download/chr20_random.fa" % cisRoot, "22_random": "%s/download/chr22_random.fa" % cisRoot, "25_random": "%s/download/chr25_random.fa" % cisRoot, "28_random": "%s/download/chr28_random.fa" % cisRoot, "Un_random": "%s/download/chrUn_random.fa" % cisRoot, "W_random": "%s/download/chrW_random.fa" % cisRoot, "E64_random": "%s/download/chrE64_random.fa" % cisRoot, "Z_random": "%s/download/chrZ_random.fa" % cisRoot, "E22C19W28_E50C23_random": "%s/download/chrE22C19W28_E50C23_random.fa" % cisRoot } print "Creating database %s" % db createDBFile(db) print "Adding gene entries" loadGeneEntries(db, genePath) #print "Adding gene annotations" #loadGeneAnnotations(db, annotPath) print "Adding gene features" loadGeneFeatures(db, genePath) for chromID in chromos.keys(): print "Loading chromosome %s" % chromID loadChromosome(db, chromID, chromos[chromID], "/G_gallus/chromo%s.bin" % chromID) print "Adding gene ontology" loadGeneOntology(db, goPath, goDefPath) print "Creating Indices" createDBindices(db) print "Finished creating database %s" % db