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