8 print "psyco not running"
10 from cistematic.genomes import Genome
11 from commoncode import getConfigParser, getConfigIntOption
18 verstring = "getsplicefa: version 1.1"
22 usage = "usage: python %prog 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
26 parser = makeParser(usage)
27 parser.set_defaults(doVerbose=False, spacer=2)
28 (options, args) = parser.parse_args(argv[1:])
32 datafilename = args[1]
39 getSpliceFasta(genome, datafilename, outfilename, maxBorder, options.doVerbose, options.spacer, delimiter)
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")
49 configParser = getConfigParser()
50 section = "getsplicefa"
51 doVerbose = getConfigIntOption(configParser, section, "doVerbose", False)
52 spacer = getConfigIntOption(configParser, section, "spacer", 2)
54 parser.set_defaults(doVerbose=doVerbose, spacer=spacer)
59 def getSpliceFasta(genome, datafilename, outfilename, maxBorder, doVerbose=False, spacer=2, delimiter="|"):
60 spacerseq = "N" * spacer
62 datafile = open(datafilename)
70 nameToComplementDict = {}
77 spliceCount = int(fields[7]) - 1
81 counter += spliceCount
82 spliceCountDict[name] = spliceCount
87 nameToChromDict[name] = chrom
88 if chrom not in alreadySeen:
89 alreadySeen[chrom] = []
91 nameToComplementDict[name] = fields[2]
94 for val in fields[8].split(",")[:-1]:
95 exonStarts.append(int(val))
97 for val in fields[9].split(",")[:-1]:
98 exonStops.append(int(val))
100 exonStartDict[name] = exonStarts
101 exonStopDict[name] = exonStops
103 for index in range(spliceCount + 1):
104 exonLengths.append(exonStops[index] - exonStarts[index] + 1)
106 exonLengthDict[name] = exonLengths
108 print len(spliceCountDict)
115 outfile = open(outfilename, "w")
116 for name in nameToChromDict:
118 spliceCount = spliceCountDict[name]
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]:
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:
138 if (beforeLen + afterLen) < 2 * maxBorder:
141 if beforeLen > maxBorder:
142 beforeLen = maxBorder
144 if afterLen > maxBorder:
148 beforeSplice = hg.sequence(chrom, exonStops[index] - maxBorder, maxBorder)
149 afterSplice = hg.sequence(chrom, exonStarts[index + 1], maxBorder)
152 print "could not get chr%s:%d-%d" % (chrom, exonStops[index], exonStarts[index + 1])
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)
162 if spliceCounter > 10000:
163 print "%d genes" % splicefileindex
168 print "%d splices too short to be seen" % missedCount
169 print "%d splices will be under-reported" % depressedCount
171 if __name__ == "__main__":