X-Git-Url: http://woldlab.caltech.edu/gitweb/?p=erange.git;a=blobdiff_plain;f=cistematic%2Fgenomes%2Fhsapiens.py;h=b0948f1327ba90b315a2b0bc245caeb68ed4c77d;hp=e9c677cace13c0e556d0c2234029d3ce42e8c98d;hb=HEAD;hpb=4ad5495359e4322da39868020a7398676261679e diff --git a/cistematic/genomes/hsapiens.py b/cistematic/genomes/hsapiens.py index e9c677c..b0948f1 100644 --- a/cistematic/genomes/hsapiens.py +++ b/cistematic/genomes/hsapiens.py @@ -1,7 +1,7 @@ ########################################################################### # # # C O P Y R I G H T N O T I C E # -# Copyright (c) 2003-10 by: # +# Copyright (c) 2003-13 by: # # * California Institute of Technology # # # # All Rights Reserved. # @@ -34,34 +34,56 @@ from cistematic.core.geneinfo import geneinfoDB from os import environ if environ.get("CISTEMATIC_ROOT"): - cisRoot = environ.get("CISTEMATIC_ROOT") + cisRoot = environ.get("CISTEMATIC_ROOT") else: cisRoot = "/proj/genome" geneDB = "%s/H_sapiens/hsapiens.genedb" % cisRoot -def loadChromosome(db, chromID, chromPath, chromOut): - seqArray = [] - hsGenome = Genome("hsapiens", 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() - hsGenome.addChromosomeEntry(chromID, chromOut, "file") +def buildHsapiensDB(db=geneDB, downloadDir="%s/download" % cisRoot): + genePath = "%s/seq_gene.md" % downloadDir # ftp://ftp.ncbi.nih.gov/genomes/H_sapiens/mapview/seq_gene.md.gz + goDefPath = "%s/GO.terms_and_ids" % downloadDir # ftp://ftp.geneontology.org/go/doc/GO.terms_and_ids + goPath = "%s/gene2go" % downloadDir # ftp://ftp.ncbi.nih.gov/gene/gene2go.gz + # chromosomes are from UCSC - will ignore all the alternative haplotypes, chrUn, and random chromosomes + chromList = ["1", "2", "3", "4", "5", "6", "7", "8", "9", "10", + "11", "12", "13", "14", "15", "16", "17", "18", "19", "20", + "21", "22", "X", "Y" + ] + + print "Creating database %s" % db + createDBFile(db) + + print "Adding gene entries" + loadGeneEntries(db, genePath, chromList) + + print "Adding gene features" + loadGeneFeatures(db, genePath, chromList) + + print "Adding gene annotations" + loadGeneAnnotations(db) + + print "Adding gene ontology" + loadGeneOntology(db, goPath, goDefPath) + + for chromID in chromList: + print "Loading chromosome %s" % chromID + chromPath = "%s/chr%s.fa" % (downloadDir, chromID) + loadChromosome(db, chromID, chromPath, "/H_sapiens/chromo%s.bin" % chromID) + + print "Creating Indices" + createDBindices(db) + + print "Finished creating database %s" % db + + +def createDBFile(db): + hsGenome = Genome("hsapiens", dbFile=db) + hsGenome.createGeneDB(db) def loadGeneEntries(db, gFile, cDict): - """ FIXME - NEED TO DEAL WITH ALTERNATIVE SPLICING ENTRIES - """ + #TODO: - NEED TO DEAL WITH ALTERNATIVE SPLICING ENTRIES geneEntries = [] hsGenome = Genome("hsapiens", dbFile=db) geneFile = open(gFile, "r") @@ -76,7 +98,7 @@ def loadGeneEntries(db, gFile, cDict): chrom = cols[1].strip() if chrom not in cDict: continue - + name = cols[10].split(":") gid = name[1] start = int(cols[2]) @@ -105,8 +127,10 @@ def loadGeneFeatures(db, gFile, cDict): cols = line.split("\t") if cols[11] not in ["CDS", "UTR", "PSEUDO"]: continue + if cols[12] == "Celera": continue + chrom = cols[1].strip() if chrom not in cDict: continue @@ -184,7 +208,7 @@ def loadGeneOntology(db, goPath, goDefPath): synonyms = idb.geneIDSynonyms(gID) if len(synonyms) >0: for entry in synonyms: - gene_name += "," + gene_name += "," gene_name += entry else: gene_name = " " @@ -198,67 +222,25 @@ def loadGeneOntology(db, goPath, goDefPath): hsGenome.addGoInfoBatch(goArray) -def createDBFile(db): - hsGenome = Genome("hsapiens", dbFile=db) - hsGenome.createGeneDB(db) +def loadChromosome(db, chromID, chromPath, chromOut): + seqArray = [] + hsGenome = Genome("hsapiens", 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() + hsGenome.addChromosomeEntry(chromID, chromOut, "file") def createDBindices(db): hsGenome = Genome("hsapiens", dbFile=db) hsGenome.createIndices() - - -def buildHsapiensDB(db=geneDB, downloadDir="%s/download" % cisRoot): - genePath = "%s/seq_gene.md" % downloadDir # ftp://ftp.ncbi.nih.gov/genomes/H_sapiens/mapview/seq_gene.md.gz - goDefPath = "%s/GO.terms_and_ids" % downloadDir # ftp://ftp.geneontology.org/go/doc/GO.terms_and_ids - goPath = "%s/gene2go" % downloadDir # ftp://ftp.ncbi.nih.gov/gene/gene2go.gz - # chromosomes are from UCSC - will ignore all the alternative haplotypes, chrUn, and 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, - "20": "%s/chr20.fa" % downloadDir, - "21": "%s/chr21.fa" % downloadDir, - "22": "%s/chr22.fa" % downloadDir, - "X": "%s/chrX.fa" % downloadDir, - "Y": "%s/chrY.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], "/H_sapiens/chromo%s.bin" % chromID) - - print "Creating Indices" - createDBindices(db) - - print "Finished creating database %s" % db