first pass cleanup of cistematic/genomes; change bamPreprocessing
[erange.git] / cistematic / genomes / mmusculus.py
index 4e90cfe7e271f5b704d6cb152bb83da2b0264ad9..72f1350070ff7027118e5a9c569e1da23d5e8e41 100644 (file)
@@ -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,36 +34,58 @@ 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/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())
+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
+    chromList = ["1", "2", "3", "4", "5", "6", "7", "8", "9", "10",
+                 "11", "12", "13", "14", "15", "16", "17", "18", "19",
+                 "X", "Y", "M"
+    ]
 
-    seq = string.join(seqArray, "")
-    seqLen = len(seq)
-    if seqLen < 1:
-        print "Problems reading sequence from file"
+    print "Creating database %s" % db
+    createDBFile(db)
 
-    print "writing to file %s" % chromOut
-    outFile = open("%s%s" % (cisRoot, chromOut), "w")
-    outFile.write(seq)
-    outFile.close()
-    mmGenome.addChromosomeEntry(chromID, chromOut, "file")
+    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, "/M_musculus/chromo%s.bin" % chromID)
+
+    print "Creating Indices"
+    createDBindices(db)
+
+    print "Finished creating database %s" % db
+
+
+def createDBFile(db):
+    mmGenome = Genome("mmusculus",  dbFile=db)
+    mmGenome.createGeneDB(db)
 
 
 def loadGeneEntries(db, gFile, cDict):
-    """ FIXME - NEED TO DEAL WITH ALTERNATIVE SPLICING ENTRIES
-    """
+    #TODO: - NEED TO DEAL WITH ALTERNATIVE SPLICING ENTRIES
+
     geneEntries = []
     alreadySeen = []
     mmGenome = Genome("mmusculus", dbFile=db)
@@ -184,7 +206,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,66 +220,26 @@ def loadGeneOntology(db, goPath, goDefPath):
     mmGenome.addGoInfoBatch(goArray)
 
 
-def createDBFile(db):
-    mmGenome = Genome("mmusculus",  dbFile=db)
-    mmGenome.createGeneDB(db)
-
-
-def createDBindices(db):
+def loadChromosome(db, chromID, chromPath, chromOut):
+    seqArray = []
     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)
+    inFile = open(chromPath, "r")
+    line = inFile.readline()
+    for line in inFile:
+        seqArray.append(line.strip())
 
-    print "Adding gene ontology"
-    loadGeneOntology(db, goPath, goDefPath)
+    seq = string.join(seqArray, "")
+    seqLen = len(seq)
+    if seqLen < 1:
+        print "Problems reading sequence from file"
 
-    for chromID in chromDict.keys():
-        print "Loading chromosome %s" % chromID
-        loadChromosome(db, chromID, chromDict[chromID], "/M_musculus/chromo%s.bin" % chromID)
+    print "writing to file %s" % chromOut
+    outFile = open("%s%s" % (cisRoot, chromOut), "w")
+    outFile.write(seq)
+    outFile.close()
+    mmGenome.addChromosomeEntry(chromID, chromOut, "file")
 
-    print "Creating Indices"
-    createDBindices(db)
 
-    print "Finished creating database %s" % db
\ No newline at end of file
+def createDBindices(db):
+    mmGenome = Genome("mmusculus", dbFile=db)
+    mmGenome.createIndices()