cistematic/genomes/athaliana.py

   1 ###########################################################################
   2 #                                                                         #
   3 # C O P Y R I G H T   N O T I C E                                         #
   4 #  Copyright (c) 2003-10 by:                                              #
   5 #    * California Institute of Technology                                 #
   6 #                                                                         #
   7 #    All Rights Reserved.                                                 #
   8 #                                                                         #
   9 # Permission is hereby granted, free of charge, to any person             #
  10 # obtaining a copy of this software and associated documentation files    #
  11 # (the "Software"), to deal in the Software without restriction,          #
  12 # including without limitation the rights to use, copy, modify, merge,    #
  13 # publish, distribute, sublicense, and/or sell copies of the Software,    #
  14 # and to permit persons to whom the Software is furnished to do so,       #
  15 # subject to the following conditions:                                    #
  16 #                                                                         #
  17 # The above copyright notice and this permission notice shall be          #
  18 # included in all copies or substantial portions of the Software.         #
  19 #                                                                         #
  20 # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,         #
  21 # EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF      #
  22 # MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND                   #
  23 # NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS     #
  24 # BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN      #
  25 # ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN       #
  26 # CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE        #
  27 # SOFTWARE.                                                               #
  28 ###########################################################################
  29 #
  30 # data for Arabidopsis thaliana
  31 import string
  32 from cistematic.genomes import Genome
  33 from os import environ
  34
  35 if environ.get("CISTEMATIC_ROOT"):
  36     cisRoot = environ.get("CISTEMATIC_ROOT")
  37 else:
  38     cisRoot = "/proj/genome"
  39
  40 geneDB = "%s/A_thaliana/athaliana.genedb" % cisRoot
  41
  42 chromSize = {"1": 30432563,
  43              "2": 19705359,
  44              "3": 23470805,
  45              "4": 18585042,
  46              "5": 26992738
  47 }
  48
  49 background = [0.3180185, 0.1819815, 0.1819815, 0.3180185]
  50 genomeSize = 119186497
  51
  52
  53 def decodeGFF3(cols):
  54     fType = cols[2]
  55     chrom = cols[0][3:]
  56     start = int(cols[3])
  57     stop = int(cols[4])
  58     sense = cols[6]
  59     if sense == "+":
  60         sense = "F"
  61     else:
  62         sense = "R"
  63
  64     other = cols[-1]
  65     otherList = other.split(";")
  66     otherDict = {}
  67     for otherItem in otherList:
  68         try:
  69             (name, value) = otherItem.split("=")
  70         except:
  71             continue
  72
  73         otherDict[name] = value
  74         if name == "Name":
  75             gid = value.strip()
  76
  77         if name == "Parent":
  78             if "," in value:
  79                 value = value.split(",")[0]
  80
  81             gid = value.strip()
  82
  83     return (fType, gid, chrom, start, stop, sense, otherDict)
  84
  85
  86 def loadChromosome(db, chromID, chromPath, chromOut):
  87     seqArray = []
  88     atGenome = Genome("athaliana", dbFile=db)
  89
  90     inFile = open(chromPath, "r")
  91     line = inFile.readline()
  92     for line in inFile:
  93         seqArray.append(line.strip())
  94
  95     seq = string.join(seqArray,"")
  96     seqLen = len(seq)
  97     if seqLen < 1:
  98         print "Problems reading sequence from file"
  99
 100     print "writing to file %s" % (chromOut)
 101     outFile = open(cisRoot + chromOut, "w")
 102     outFile.write(seq)
 103     outFile.close()
 104     seq = ""
 105
 106     atGenome.addChromosomeEntry(chromID, chromOut, "file")
 107     # Add alternative chromID - should be A-O and 01-09
 108     atGenome.addChromosomeEntry("chromo%s" % chromID, chromOut, "file")
 109
 110
 111 def loadGeneEntries(db, gFile):
 112     geneEntries = []
 113     atGenome = Genome("athaliana", dbFile=db)
 114     geneFile = open(gFile, "r")
 115     for line in geneFile:
 116         if line[0] == "#" or len(line) < 10:
 117             continue
 118
 119         fields = line.strip().split("\t")
 120         if fields[2] != "gene":
 121             continue
 122
 123         (fType, gid, chrom, start, stop, sense, otherDict) = decodeGFF3(fields)
 124         geneID = ("athaliana", gid)
 125         version = 1
 126         geneEntries.append((geneID, chrom, start, stop, sense, "gene", version))
 127
 128     print "inserting %d gene entries" % len(geneEntries)
 129     atGenome.addGeneEntryBatch(geneEntries)
 130
 131
 132 def loadFeatureEntries(db, gFile):
 133     featureEntries = []
 134     trackedGenes = []
 135     atGenome = Genome("athaliana", dbFile=db)
 136     featureTranslation = {"CDS": "CDS",
 137                           "three_prime_UTR": "3UTR",
 138                           "five_prime_UTR": "5UTR",
 139                           "miRNA": "5UTR",
 140                           "exon": "5UTR"
 141     }
 142     geneFile = open(gFile, "r")
 143     for line in geneFile:
 144         fields = line.split("\t")
 145         (fType, gid, chrom, start, stop, sense, otherDict) = decodeGFF3(fields)
 146         if fType in ["ncRNA"]:
 147             (locus, rev) = gid.split(".")
 148             if gid not in trackedGenes:
 149                 trackedGenes.append(locus)
 150
 151     geneFile = open(gFile, "r")
 152     for line in geneFile:
 153         if line[0] == "c" or len(line) < 10:
 154             continue
 155
 156         fields = line.split("\t")
 157         (fType, gid, chrom, start, stop, sense, otherDict) = decodeGFF3(fields)
 158         locusField = gid.split('.')
 159         try:
 160             (locus, rev) = locusField
 161             rev = rev.strip()
 162         except:
 163             locus = gid
 164             rev = 1
 165
 166         if fType not in featureTranslation:
 167             continue
 168
 169         elif fType == "exon" and locus not in trackedGenes:
 170             continue
 171
 172         geneID = ("athaliana", locus)
 173         featureEntries.append((geneID, rev, chrom, start, stop, sense, featureTranslation[fType]))
 174
 175     print "inserted %d feature entries" % len(featureEntries)
 176     atGenome.addFeatureEntryBatch(featureEntries)
 177
 178
 179 def loadGeneAnnotations(db, annotPath):
 180     geneAnnotations = []
 181     annotFile = open(annotPath, "r")
 182     annotFile.readline()
 183     lines = annotFile.readlines()
 184     annotFile.close()
 185
 186     atGenome = Genome("athaliana", dbFile=db)
 187     for line in lines:
 188         field = line.split("\t")
 189         try:
 190             orfName = field[0].strip()
 191             if "." in orfName:
 192                 (locus, rev) = orfName.split(".")
 193                 orfName = locus
 194
 195             description = field[2].strip()
 196             geneAnnotations.append((("athaliana", orfName), string.replace(description, "'", "p")))
 197         except:
 198             pass
 199
 200     print "Adding %d annotations" % len(geneAnnotations)
 201     atGenome.addAnnotationBatch(geneAnnotations)
 202
 203
 204 def loadGeneOntology(db, goPath):
 205     atGenome = Genome("athaliana", dbFile=db)
 206     goFile = open(goPath, "r")
 207     goEntries = goFile.readlines()
 208     goArray = []
 209     for line in goEntries:
 210         fields = line.split("\t")
 211         gID = fields[0]
 212         GOID = fields[4]
 213         objType = string.replace(fields[3], "'", "p")
 214         objName = string.replace(fields[2], "'", "p")
 215         isNot = ""
 216         GOterm  = fields[7]
 217         evidence = fields[8]
 218         goArray.append((("athaliana", gID), GOID[3:], objType, objName, isNot, GOterm, evidence, fields[9]))
 219
 220     atGenome.addGoInfoBatch(goArray)
 221
 222
 223 def createDBFile(db):
 224     atGenome = Genome("athaliana", dbFile=db)
 225     atGenome.createGeneDB(db)
 226
 227
 228 def createDBindices(db):
 229     atGenome = Genome("athaliana", dbFile=db)
 230     atGenome.createIndices()
 231
 232
 233 def buildArabidopsisDB(db=geneDB, downloadDir="%s/download" % cisRoot):
 234     genePath = "%s/TAIR9_GFF3_genes_transposons.gff" % downloadDir
 235     annotPath = "%s/TAIR9_functional_descriptions" % downloadDir
 236     goPath = "%s/ATH_GO_GOSLIM.txt" % downloadDir
 237
 238     chromos = {"1": "%s/chr1.fas" % downloadDir,
 239                "2": "%s/chr2.fas" % downloadDir,
 240                "3": "%s/chr3.fas" % downloadDir,
 241                "4": "%s/chr4.fas" % downloadDir,
 242                "5": "%s/chr5.fas" % downloadDir,
 243                "C": "%s/chrC.fas" % downloadDir,
 244                "M": "%s/chrM.fas" % downloadDir
 245     }
 246
 247     print "Creating database %s" % db
 248     createDBFile(db)
 249
 250     print "Adding gene entries"
 251     loadGeneEntries(db, genePath)
 252
 253     print "Adding feature entries"
 254     loadFeatureEntries(db, genePath)
 255
 256     print "Adding gene annotations"
 257     loadGeneAnnotations(db, annotPath)
 258
 259     print "Adding gene ontology"
 260     loadGeneOntology(db, goPath)
 261
 262     for chromID in chromos.keys():
 263         print "Loading chromosome %s" % chromID
 264         loadChromosome(db, chromID, chromos[chromID],  "/A_thaliana/chr%s.bin" % chromID)
 265
 266     print "Creating Indices"
 267     createDBindices(db)
 268
 269     print "Finished creating database %s" % db