Skip to content

Commit ceaa185

Browse files
author
Mark Fiers
committed
ongoing work
1 parent f68df07 commit ceaa185

File tree

4 files changed

+149
-163
lines changed

4 files changed

+149
-163
lines changed

‎hagfishUtils.py

+49-2
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
11
import os
2+
import sys
23
import math
34
import copy
45

@@ -25,6 +26,28 @@
2526
COLOR6 = '#EE1904' #lighter red
2627
COLOR7 = '#32AE1C' # lighter green
2728

29+
COLDARKGREEN = '#007D1C'
30+
COLDARKRED = '#A61000'
31+
COLDARKBLUE = '#015666'
32+
COLDARKYELLOW = '#A65100'
33+
34+
COLLIGHTGREEN = '#38E05D'
35+
COLLIGHTRED = '#FF5240'
36+
COLLIGHTBLUE = '#37B6CE'
37+
COLLIGHTYELLOW = '#FF9D40'
38+
39+
COLDARKGREEN = '#007D1C'
40+
COLDARKRED = '#9F0025'
41+
COLDARKPURPLE = '#42036F'
42+
COLDARKBLUE = '#06266F'
43+
COLDARKYELLOW = '#A66F00'
44+
45+
COLLIGHTGREEN = '#38E05D'
46+
COLLIGHTRED = '#B12D4C'
47+
COLLIGHTBLUE = '#4671D5'
48+
COLLIGHTYELLOW = '#FFBF40'
49+
COLLIGHTPURPLE = '#963FD5'
50+
2851
COLMAP1 = mpl.colors.LinearSegmentedColormap.from_list(
2952
'COLMAP1',
3053
[COLOR6, COLOR3, COLOR7], N=200)
@@ -144,10 +167,22 @@ def __init__(self, options, args, seqId = None, inputFile=None):
144167
self.data['ok'] = self.ok
145168
self.data['low'] = self.low
146169
self.data['high'] = self.high
170+
147171
self.data['ok_ends'] = self.ok_ends
148172
self.data['low_ends'] = self.low_ends
149173
self.data['high_ends'] = self.high_ends
150174

175+
#see if there is gapdata to load
176+
gapbase = os.path.join('gaps', self.seqId)
177+
if np_exists(gapbase, 'nns'):
178+
self.l.info("loading gap data")
179+
self.nns = np_load(gapbase, 'nns')
180+
self.vectors.append('nns')
181+
else:
182+
self.l.info("Cannot find gap data")
183+
self.nns = None
184+
self.data['nns'] = self.nns
185+
151186
self.all = self.low + self.ok + self.high
152187

153188
self.seqLen = len(self.ok)
@@ -184,15 +219,24 @@ def __init__(self, options, data, title=None, data2=None, ymax=None):
184219
self.start = 0
185220
if options.start:
186221
self.start = int(float(options.start))
222+
if self.start > data.seqLen:
223+
self.l.critical("start plot is greater than sequence length")
224+
sys.exit(-1)
187225

188226
self.stop = self.data.seqLen
189227
if options.stop:
190228
self.stop = int(float(options.stop))
191229
if self.stop > data.seqLen:
192230
self.stop = data.seqLen
231+
self.l.info("plot stops at %d" % self.stop)
232+
233+
if self.start > self.stop:
234+
self.l.critical("plot start is greater than plot stop")
235+
sys.exit(-1)
193236

194237
self.plotLen = self.stop - self.start
195238
self.ntPerBand = int(float(options.ntPerBand))
239+
196240
if self.ntPerBand == -1:
197241
self.l.debug("nt per band is not specified")
198242
self.ntPerBand = int(1e6)
@@ -303,16 +347,18 @@ def plotBands(self, bandPlotter):
303347

304348
def save(self, tag=""):
305349
self.ax.set_xlim(0, self.ntPerBand)
350+
306351
if self.options.outfile:
307352
outFileName = self.options.outfile
308353
else:
309354
outFileName = self.data.seqId
310355

356+
if self.options.start or self.options.stop:
357+
outFileName += "_%d_%d" % (self.start, self.stop)
358+
311359
if tag:
312360
outFileName += '_%s' % tag
313361

314-
if self.options.start or self.options.stop:
315-
outFileName += "_%d_%d" % (self.start, self.stop)
316362

317363
for f in self.options.format:
318364
self.l.info("writing to %s.%s" % (outFileName, f))
@@ -352,6 +398,7 @@ def __init__(self, plot, band, options):
352398
#special form of x
353399
self.locx = self.x - self.start
354400
self.zero = np.zeros_like(self.locx)
401+
355402
def plot(self):
356403
"""
357404
Should be overridden

‎hagfish_cplot

+52-156
Original file line numberDiff line numberDiff line change
@@ -23,6 +23,11 @@ parser = getHagfishOptparser()
2323
parser.add_option('-S', dest='score', action='store_true',
2424
help='plot the hagfish score' )
2525

26+
#parser.add_option('--pe', dest='pairedends', action='store_true',
27+
# help='plot only based on paired ends')
28+
#parser.add_option('--pe', dest='pairedends', action='store_true',
29+
# help='plot only based on paired ends')
30+
2631
addPlotParameters(parser)
2732
options, args = parser.parse_args()
2833

@@ -46,165 +51,56 @@ class plotter(hagfishPlotBand):
4651
self.l.debug("Start plotting")
4752
yc = self.yCorrection
4853

49-
## print lightly colored high/low/good
50-
## to server a background
51-
#self.l.debug("Plot background")
52-
53-
# if PLOTHIGH:
54-
# ## bg high
55-
# self.ax.fill_between(
56-
# self.locx,
57-
# self.okh + yc,
58-
# self.okh + self.high + yc,
59-
# alpha = 1,
60-
# zorder = 0,
61-
# color=COLOR2)
62-
63-
# if PLOTLOW:
64-
# ### bg low
65-
# self.ax.fill_between(
66-
# self.locx,
67-
# -self.okh + yc,
68-
# (-self.okh - self.low) + yc,
69-
# alpha=1,
70-
# color=COLOR3,
71-
# zorder=0)
72-
73-
# if PLOTOK:
74-
# ### bg ok
75-
# self.ax.fill_between(
76-
# self.locx,
77-
# -self.okh + yc,
78-
# self.okh + yc,
79-
# alpha=1,
80-
# color=COLOR1,
81-
# zorder=0)
82-
83-
# if PLOTHIGH:
84-
# ## outlines
85-
# self.ax.plot(
86-
# self.locx,
87-
# self.okh + self.high + yc,
88-
# color = 'black', linewidth=0.5,
89-
# zorder=0)
90-
# if PLOTLOW:
91-
# self.ax.plot(
92-
# self.locx,
93-
# -(self.okh + self.low) + yc,
94-
# color = 'black',
95-
# linewidth=0.5,
96-
# zorder=0)
97-
98-
## Reprint High/ok/low - now with a cutof to stay
99-
## in the band
100-
101-
if PLOTHIGH:
102-
#high
103-
self.ax.fill_between(
104-
self.locx,
105-
np.minimum(self.bandTop, self.okh + yc),
106-
np.minimum(self.bandTop, self.okh + self.high + yc),
107-
alpha = 0.8,\
108-
zorder = 1,
109-
color=COLOR2)
110-
111-
if PLOTLOW:
112-
#low
113-
self.ax.fill_between(
114-
self.locx,
115-
np.maximum(self.bandBottom, -self.okh + yc),
116-
np.maximum(self.bandBottom, (-self.okh - self.low) + yc),
117-
alpha=0.8,
118-
color=COLOR3,
119-
zorder=1)
120-
#low
121-
self.ax.fill_between(
122-
self.locx,
123-
np.maximum(self.bandBottom, yc),
124-
np.maximum(self.bandBottom, -self.low_ends + yc),
125-
alpha=0.8,
126-
color="black",
127-
zorder=12)
128-
129-
if PLOTOK:
130-
#ok
131-
self.ax.fill_between(
132-
self.locx,
133-
np.minimum(self.bandTop, self.okh + yc),
134-
np.maximum(self.bandBottom, -self.okh + yc),
135-
alpha=0.8,
136-
color=COLOR1,
137-
zorder=1)
138-
54+
#print 'high'
55+
self.ax.fill_between(
56+
self.locx,
57+
np.minimum(self.bandTop, self.okh + yc),
58+
np.minimum(self.bandTop, self.okh + self.high + yc),
59+
zorder = 1,
60+
color=COLDARKRED)
61+
#print 'ok'
62+
self.ax.fill_between(
63+
self.locx,
64+
np.minimum(self.bandTop, self.okh + yc),
65+
np.maximum(self.bandBottom, -self.okh + yc),
66+
color=COLDARKGREEN,
67+
zorder=1)
68+
69+
#low
70+
self.ax.fill_between(
71+
self.locx,
72+
np.maximum(self.bandBottom, -self.okh + yc),
73+
np.maximum(self.bandBottom, -(self.okh + self.low) + yc),
74+
color=COLDARKBLUE,
75+
zorder=1)
76+
77+
78+
#ok_ends
13979
# self.ax.fill_between(
14080
# self.locx,
141-
# self.okh_ends + yc,
142-
# -self.okh_ends + yc,
143-
# alpha=0.3,
144-
# color=COLOR4,
145-
# zorder=2)
146-
147-
# self.ax.fill_between(self.locx,
148-
# self.okh + self.high + yc,
149-
# self.okh + self.high + self.high_ends + yc,
150-
# color=COLOR4,
151-
# linewidth=1,
152-
# zorder=4)
153-
154-
# self.ax.fill_between(self.locx,
155-
# -self.okh - self.low + yc,
156-
# -self.okh - self.low - self.low_ends + yc,
157-
# color=COLOR4,
158-
# linewidth=1,
159-
# zorder=4)
160-
#self.ax.plot(self.locx, -self.low_ends + yc, color=COLOR3,
161-
# linewidth=0.5, zorder=12, alpha=0.5)
81+
# np.minimum(self.bandTop, self.okh_ends + yc),
82+
# np.maximum(self.bandBottom, -self.okh_ends + yc),
83+
# color=COLDARKGREEN,
84+
# zorder=1)
85+
86+
# #print 'high'
16287
# self.ax.fill_between(
163-
# self.locx,
164-
# self.okh_ends + yc,
165-
# self.okh_ends + self.high_ends + yc,
166-
# alpha=0.3,
167-
# color='black',
168-
# zorder=2)
169-
170-
#plot new outlines
171-
#filter them, so that we only have outlines where the
172-
#plot remains inside the band
173-
174-
# if PLOTHIGH:
175-
# self.ax.plot(
176-
# self.locx,
177-
# self.okh + self.high + yc,
178-
# color = 'black', linewidth = 0.5,
179-
# zorder =5,
180-
# clip_path=self.patch
181-
# )
182-
183-
# if PLOTOK:
184-
# self.ax.plot(
185-
# self.locx,
186-
# self.okh + yc,
187-
# color = 'black', linewidth = 0.5,
188-
# zorder =5,
189-
# clip_path=self.patch
190-
# )
191-
192-
# self.ax.plot(
193-
# self.locx,
194-
# -self.okh + yc,
195-
# color = 'black',
196-
# linewidth=0.5,
197-
# clip_path=self.patch,
198-
# zorder=5)
199-
200-
# if PLOTLOW:
201-
# self.ax.plot(
202-
# self.locx,
203-
# -(self.okh + self.low) + yc,
204-
# color = 'black',
205-
# linewidth=0.5,
206-
# clip_path=self.patch,
207-
# zorder=5)
88+
# self.locx,
89+
# np.minimum(self.bandTop, self.okh_ends + yc),
90+
# np.minimum(self.bandTop, self.okh_ends + self.high_ends + yc),
91+
# zorder = 1,
92+
# color=COLDARKRED)
93+
94+
95+
96+
if self.__dict__.has_key('nns'):
97+
self.l.debug("plotting N band")
98+
self.ax.fill_between(
99+
self.locx,
100+
self.zero + yc + (self.plot.maxY * self.nns),
101+
self.zero + yc - (self.plot.maxY * self.nns),
102+
color=COLLIGHTYELLOW, alpha=0.6)
103+
208104

209105

210106
if options.score:

‎hagfish_file_util.py

+39
Original file line numberDiff line numberDiff line change
@@ -1,3 +1,5 @@
1+
import os
2+
import bz2
13
import gzip
24
import cPickle
35
import logging
@@ -21,6 +23,10 @@ def np_savez(file_base, **kwargs):
2123
cPickle.dump(kwargs[k], F)
2224
F.close()
2325

26+
def np_exists(file_base, what):
27+
fn = '%s.%s.gz' % (file_base, what)
28+
return os.path.exists(fn)
29+
2430
def np_load(file_base, what):
2531
"""
2632
reimplement the crappy np.savez function
@@ -32,3 +38,36 @@ def np_load(file_base, what):
3238
F.close()
3339
return data
3440

41+
42+
43+
def fastareader(f):
44+
if type(f) == type('hi'):
45+
#it is probably a filename
46+
if f[-4:] == '.bz2':
47+
F = bz2.BZ2File(f)
48+
else:
49+
F = open(f, 'r')
50+
else:
51+
#probably a file handle, or stdin
52+
F = f
53+
54+
name, seq = "", []
55+
while True:
56+
l = F.readline()
57+
if not l: break
58+
59+
l = l.strip()
60+
if not l: continue
61+
62+
if l[0] == '>':
63+
if name and seq:
64+
yield name, "".join(seq)
65+
seq = []
66+
name = l[1:].split()[0]
67+
else:
68+
seq.append("".join(l.split()).lower())
69+
70+
if name and seq:
71+
yield name, "".join(seq)
72+
73+
F.close()

0 commit comments

Comments
 (0)