kuroの覚え書き

96の個人的覚え書き

sequence alignmentをBokehを使ってインタラクティブに表示してみる

clustalwでsequenceのアライメントをとると、テキストで.alnというファイルが生成されるが、文字の並びとアスタリスクではわかりにくいことが多い。
なのでclustalxだとかMEGAだとかで表示するとカラフルに色分けで表示できるのでぱっと直感的に分かるのだが、いちいちファイルをアプリで開き直すとか、面倒くさいし。
ということでpythonでそういうのができないかやってみた。できたらそれをいつものウェブアプリに実装してやる。
Bioinformatics and other bits - A sequence alignment viewer with Bokeh and Panel
こちらの記事を参考に

import string
import numpy as np
from Bio import AlignIO
import panel as pn
pn.extension()
from bokeh.plotting import figure, output_file, save
from bokeh.models import ColumnDataSource, Plot, Grid, Range1d
from bokeh.models.glyphs import Text, Rect
from bokeh.layouts import gridplot

#output_file("out.html")
aln = AlignIO.read('Phylo.aln','clustal')
seqs = [rec.seq for rec in (aln)]
ids = [rec.id for rec in aln]    
text = [i for s in list(seqs) for i in s]
clrs =  {'A':'red','T':'green','G':'orange','C':'blue','-':'white'}
colors = [clrs[i] for i in text]
N = len(seqs[0])
S = len(seqs)
width = .4

x = np.arange(1,N+1)
y = np.arange(0,S,1)
xx, yy = np.meshgrid(x, y)

gx = xx.ravel()
gy = yy.flatten()

recty = gy+.5
h= 1/S

source = ColumnDataSource(dict(x=gx, y=gy, recty=recty, text=text, colors=colors))
plot_height = len(seqs)*15+50
x_range = Range1d(0,N+1, bounds='auto')
plot_width=800
fontsize="9pt"
if N>100:
    viewlen=100
else:
     viewlen=N

view_range = (0,viewlen)
tools="xpan, xwheel_zoom, reset, save"

p = figure(title=None, plot_width= plot_width, plot_height=50,
        x_range=x_range, y_range=(0,S), tools=tools,
        min_border=0, toolbar_location='below')

rects = Rect(x="x", y="recty",  width=1, height=1, fill_color="colors",
        line_color=None, fill_alpha=0.6)
p.add_glyph(source, rects)
p.yaxis.visible = False
p.grid.visible = False
p1 = figure(title=None, plot_width=plot_width, plot_height=plot_height,
                x_range=view_range, y_range=ids, tools="xpan,reset",
                min_border=0, toolbar_location='below')#, lod_factor=1)          
glyph = Text(x="x", y="y", text="text", text_align='center',text_color="black",
                text_font="monospace",text_font_size=fontsize)
rects = Rect(x="x", y="recty",  width=1, height=1, fill_color="colors",
                line_color=None, fill_alpha=0.4)
p1.add_glyph(source, glyph)
p1.add_glyph(source, rects)

p1.grid.visible = False
p1.xaxis.major_label_text_font_style = "bold"
p1.yaxis.minor_tick_line_width = 0
p1.yaxis.major_tick_line_width = 0

p = gridplot([[p],[p1]])

pn.pane.Bokeh(p)
#save(p)

できた。
f:id:k-kuro:20200902171838p:plain
Bokehって色々可能性を感じさせるね。

一番最初に
output_file("out.html")
と保存先を指定して
save(p)
でhtmlファイルとして保存してやるとwebアプリからも利用しやすいな。

f:id:k-kuro:20200902172541p:plain

系統樹とアライメントを一気に解析。ああなんて快適。