seq = SeqIO.read('/home/tboland/SAD7743_P01_E09_LC_E09.ab1', 'abi')
quality = seq.letter_annotations['phred_quality']
peak_calls = seq.annotations['abif_raw']['PLOC2']
order = seq.annotations['abif_raw']['FWO_1']
sequence = str(seq.seq)
# Extract color data
peaks = {}
for i, (base, r) in enumerate(zip(order, [9, 10, 11, 12])):
peaks[base] = seq.annotations['abif_raw']['DATA{}'.format(r)]
# Plot
p = figure(
plot_width=900,
plot_height=400,
output_backend="webgl",
tools=['xpan', 'reset','xwheel_zoom', BoxZoomTool(dimensions='width')],
active_drag='xpan',
active_scroll='xwheel_zoom',
title=seq.id
)
# Plot quality
p.extra_y_ranges = {"q": Range1d(start=0, end=100)}
p.add_layout(LinearAxis(y_range_name="q"), 'right')
p.rect(x=peak_calls, y=np.array(quality)/2, width=14, height=quality, y_range_name='q', color='#e7ecf6')
# Plot colors
for color, base in zip(['red', 'green', 'blue', 'black'], ['T', 'A', 'C', 'G']):
peak = peaks[base]
p.line(range(len(peak)), peak, line_color=color, legend=base)
# Adjust tickers/ranges
p.xaxis.ticker = FixedTicker(ticks=peak_calls)
p.xaxis.major_label_overrides = dict(zip(peak_calls, sequence))
p.y_range = Range1d(0, max(max(x) for x in peaks.values()))
show(p)