########################################################################### # # # 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 Mus musculus 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/M_musculus/mmusculus.genedb" % cisRoot def loadChromosome(db, chromID, chromPath, chromOut): seqArray = [] mmGenome = Genome("mmusculus", 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() mmGenome.addChromosomeEntry(chromID, chromOut, "file") def loadGeneEntries(db, gFile, cDict): """ FIXME - NEED TO DEAL WITH ALTERNATIVE SPLICING ENTRIES """ geneEntries = [] alreadySeen = [] mmGenome = Genome("mmusculus", dbFile=db) geneFile = open(gFile, "r") for line in geneFile: cols = line.split("\t") if cols[11].strip() != "GENE" or cols[12] != "C57BL/6J": continue chrom = cols[1].strip() if chrom not in cDict: 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] if sense == "+": sense = "F" else: sense = "R" geneID = ("mmusculus", gid) gidVersion = 1 geneEntries.append((geneID, chrom, start, stop, sense, "gene", gidVersion)) print "Adding %d gene entries" % len(geneEntries) mmGenome.addGeneEntryBatch(geneEntries) def loadGeneFeatures(db, gFile, cDict): """ Load gene features such as CDS, UTR, and PSEUDO from the gene file. """ featureEntries = [] mmGenome = Genome("mmusculus", dbFile=db) featureFile = open(gFile, "r") for line in featureFile: cols = line.split("\t") if cols[11].strip() not in ["CDS", "UTR", "PSEUDO"] or cols[12] != "C57BL/6J": continue chrom = cols[1].strip() if chrom not in cDict: 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] if sense == "+": sense = "F" else: sense = "R" geneID = ("mmusculus", gid) gidVersion = 1 featureEntries.append((geneID, gidVersion, chrom, start, stop, sense, fType)) print "Adding %d feature entries" % len(featureEntries) mmGenome.addFeatureEntryBatch(featureEntries) def loadGeneAnnotations(db): geneAnnotations = [] idb = geneinfoDB() mmGenome = Genome("mmusculus", dbFile=db) gidList = mmGenome.allGIDs() for locID in gidList: gID = ("mmusculus", locID) geneDescArray = idb.getDescription(gID) geneDesc = "" for entry in geneDescArray: geneDesc += "," geneDesc += entry.strip() if len(geneDescArray) > 0: geneAnnotations.append((gID, string.replace(geneDesc[1:], "'", "p"))) print "Adding %d annotations" % len(geneAnnotations) mmGenome.addAnnotationBatch(geneAnnotations) def loadGeneOntology(db, goPath, goDefPath): mmGenome = Genome("mmusculus", 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] != "10116": continue locID = fields[1].strip() gID = ("mmusculus", locID) if prevGID != gID: prevGID = gID gene_name = "" synonyms = idb.geneIDSynonyms(gID) if len(synonyms) >0: for entry in synonyms: gene_name += "," gene_name += entry else: gene_name = " " goArray.append((gID, fields[2], "", gene_name[1:], "", 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) mmGenome.addGoInfoBatch(goArray) def createDBFile(db): mmGenome = Genome("mmusculus", dbFile=db) mmGenome.createGeneDB(db) def createDBindices(db): mmGenome = Genome("mmusculus", dbFile=db) mmGenome.createIndices() def buildMmusculusDB(db=geneDB, downloadDir="%s/download" % cisRoot): genePath = "%s/seq_gene.md" % downloadDir # ftp://ftp.ncbi.nih.gov/genomes/M_musculus/mapview/seq_gene.md goDefPath = "%s/GO.terms_and_ids" % downloadDir # ftp://ftp.geneontology.org/pub/go/doc/GO.terms_and_ids goPath = "%s/gene2go" % downloadDir # ftp://ftp.ncbi.nih.gov/gene/DATA/gene2go.gz # chromosomes are from ftp://hgdownload.cse.ucsc.edu/goldenPath/mm9/chromosomes # but ignoring all random chromosomes chromDict = {"1": "%s/chr1.fa" % downloadDir, "2": "%s/chr2.fa" % downloadDir, "3": "%s/chr3.fa" % downloadDir, "4": "%s/chr4.fa" % downloadDir, "5": "%s/chr5.fa" % downloadDir, "6": "%s/chr6.fa" % downloadDir, "7": "%s/chr7.fa" % downloadDir, "8": "%s/chr8.fa" % downloadDir, "9": "%s/chr9.fa" % downloadDir, "10": "%s/chr10.fa" % downloadDir, "11": "%s/chr11.fa" % downloadDir, "12": "%s/chr12.fa" % downloadDir, "13": "%s/chr13.fa" % downloadDir, "14": "%s/chr14.fa" % downloadDir, "15": "%s/chr15.fa" % downloadDir, "16": "%s/chr16.fa" % downloadDir, "17": "%s/chr17.fa" % downloadDir, "18": "%s/chr18.fa" % downloadDir, "19": "%s/chr19.fa" % downloadDir, "X": "%s/chrX.fa" % downloadDir, "Y": "%s/chrY.fa" % downloadDir, "M": "%s/chrM.fa" % downloadDir } print "Creating database %s" % db createDBFile(db) print "Adding gene entries" loadGeneEntries(db, genePath, chromDict) print "Adding gene features" loadGeneFeatures(db, genePath, chromDict) print "Adding gene annotations" loadGeneAnnotations(db) print "Adding gene ontology" loadGeneOntology(db, goPath, goDefPath) for chromID in chromDict.keys(): print "Loading chromosome %s" % chromID loadChromosome(db, chromID, chromDict[chromID], "/M_musculus/chromo%s.bin" % chromID) print "Creating Indices" createDBindices(db) print "Finished creating database %s" % db