getsplicefa.py

   1 import sys
   2 import optparse
   3 import string
   4 try:
   5     import psyco
   6     psyco.full()
   7 except:
   8     print "psyco not running"
   9
  10 from cistematic.genomes import Genome
  11 from commoncode import getConfigParser, getConfigIntOption
  12
  13
  14 def main(argv=None):
  15     if not argv:
  16         argv = sys.argv
  17
  18     verstring = "getsplicefa: version 1.1"
  19     print verstring
  20     delimiter = "|"
  21
  22     usage = "usage: python getsplicefa.py genome ucscModels outfilename maxBorder [--verbose] [--spacer num]\
  23             \n\twhere ucscModels is the gene model file in same format as the UCSC Table Browser\
  24             \n\twhere spacer is by default 2, and maxBorder should be readlen - (2 * spacer)\
  25             \n\tdelimiter is set to %s - edit the code to change it, if necessary\n" % delimiter
  26
  27     parser = makeParser(usage)
  28     parser.set_defaults(doVerbose=False, spacer=2)
  29     (options, args) = parser.parse_args(argv[1:])
  30
  31     try:
  32         genome = args[0]
  33         datafilename = args[1]
  34         outfilename = args[2]
  35         maxBorder = int(args[3])
  36     except IndexError:
  37         print usage
  38         sys.exit(1)
  39
  40     getSpliceFasta(genome, datafilename, outfilename, maxBorder, options.doVerbose, options.spacer, delimiter)
  41
  42
  43 def makeParser(usage=""):
  44     parser = optparse.OptionParser(usage=usage)
  45     parser.add_option("--verbose", action="store_true", dest="doVerbose",
  46                       help="show verbose messages [default: False]")
  47     parser.add_option("--spacer", type="int", dest="spacer",
  48                       help="number of spacer NTs to use [default: 2")
  49
  50     configParser = getConfigParser()
  51     section = "getsplicefa"
  52     doVerbose = getConfigIntOption(configParser, section, "doVerbose", False)
  53     spacer = getConfigIntOption(configParser, section, "spacer", 2)
  54
  55     parser.set_defaults(doVerbose=doVerbose, spacer=spacer)
  56
  57     return parser
  58
  59
  60 def getSpliceFasta(genome, datafilename, outfilename, maxBorder, doVerbose=False, spacer=2, delimiter="|"):
  61     spacerseq = "N" * spacer
  62
  63     datafile = open(datafilename)
  64     hg = Genome(genome)
  65
  66     spliceCountDict = {}
  67     exonStartDict = {}
  68     exonStopDict = {}
  69     exonLengthDict = {}
  70     nameToChromDict = {}
  71     nameToComplementDict = {}
  72     alreadySeen = {}
  73     counter = 0
  74
  75     for line in datafile:
  76         fields = line.split()
  77         name = fields[0]
  78         spliceCount = int(fields[7]) - 1
  79         if spliceCount < 1:
  80             continue
  81
  82         counter += spliceCount
  83         spliceCountDict[name] = spliceCount
  84         chrom = fields[1][3:]
  85         if chrom == "chrM":
  86             continue
  87
  88         nameToChromDict[name] = chrom
  89         if chrom not in alreadySeen:
  90             alreadySeen[chrom] = []
  91
  92         nameToComplementDict[name] = fields[2]
  93         exonStarts = []
  94         exonStops = []
  95         for val in fields[8].split(",")[:-1]:
  96             exonStarts.append(int(val))
  97
  98         for val in fields[9].split(",")[:-1]:
  99             exonStops.append(int(val))
 100
 101         exonStartDict[name] = exonStarts
 102         exonStopDict[name] = exonStops
 103         exonLengths = []
 104         for index in range(spliceCount + 1):
 105             exonLengths.append(exonStops[index] - exonStarts[index] + 1)
 106
 107         exonLengthDict[name] = exonLengths
 108
 109     print len(spliceCountDict)
 110     print counter
 111
 112     missedCount = 0
 113     depressedCount = 0
 114     splicefileindex = 1
 115     spliceCounter = 0
 116     outfile = open(outfilename, "w")
 117     for name in nameToChromDict:
 118         try:
 119             spliceCount = spliceCountDict[name]
 120         except:
 121             continue
 122
 123         exonStarts = exonStartDict[name]
 124         exonStops = exonStopDict[name]
 125         exonLengths = exonLengthDict[name]
 126         chrom = nameToChromDict[name]
 127         for index in range(spliceCount):
 128             if (exonStops[index], exonStarts[index + 1]) in alreadySeen[chrom]:
 129                 continue
 130
 131             regionstart = exonStops[index] - maxBorder
 132             alreadySeen[chrom].append((exonStops[index], exonStarts[index + 1]))
 133             beforeLen = exonLengths[index]
 134             afterLen = exonLengths[index + 1]
 135             if (beforeLen + afterLen) < maxBorder + spacer:
 136                 missedCount += 1
 137                 continue
 138
 139             if (beforeLen + afterLen) < 2 * maxBorder:
 140                 depressedCount += 1
 141
 142             if beforeLen > maxBorder:
 143                 beforeLen = maxBorder
 144
 145             if afterLen > maxBorder:
 146                 afterLen = maxBorder
 147
 148             try:
 149                 beforeSplice = hg.sequence(chrom, exonStops[index] - maxBorder, maxBorder)
 150                 afterSplice = hg.sequence(chrom, exonStarts[index + 1], maxBorder)
 151             except:
 152                 if doVerbose:
 153                     print "could not get chr%s:%d-%d" % (chrom, exonStops[index], exonStarts[index + 1])
 154                 continue
 155
 156             sequenceHeader = string.join([name, delimiter, str(index), delimiter, str(regionstart)], "")
 157             spliceJunctionSequence = string.join([spacerseq, beforeSplice.upper(), afterSplice.upper(), spacerseq], "")
 158             outstring = ">%s\n%s\n" % (sequenceHeader, spliceJunctionSequence)
 159             outfile.write(outstring)
 160
 161         splicefileindex += 1
 162         spliceCounter += 1
 163         if spliceCounter > 10000:
 164             print "%d genes" % splicefileindex
 165             spliceCounter = 0
 166
 167     outfile.close()
 168
 169     print "%d splices too short to be seen" % missedCount
 170     print "%d splices will be under-reported" % depressedCount
 171
 172 if __name__ == "__main__":
 173     main(sys.argv)