<jupyter_start><jupyter_text>Interactive plotting with *bokeh*[`bokeh`](https://www.bokeh.pydata.org) is a plotting library that has some incredible example plots, usually written in about 10 lines of python... It has been around for a little while, but it's recent release spurred me to finally get around to trying it. Working with genomic data, one thing I find myself doing a lot, is plotting some metric across the genome, then wondering what genomic features are underlying the interesting bits. Currently, my workflow relies on static matplotlib plots. This means I make a plot, then replot iteratively with different axes- which can end up wasting  time. What `bokeh` does, is to provide interactivity. Allowing panning and zooming into the interesting bits, as well as helpful annotations.The aim of this post then, is to give an example of how to use the bokeh plotting library to create some interactive genome scale plots of two different organisms- *An. gambiae*, and *Pl. falciparum*. Two examples are shown directly below, but described in detail further down.TL;DRHow to make some neat generalisable genome plots<jupyter_code>%run boilerplate.ipynb<jupyter_output><empty_output><jupyter_text>Genome-wide Tajima's D of *Anopheles gambiae* from Burkina FasoTry hovering over a data point, or clicking on one of the panels on the left to wheel zoom, or select zoom!The hover label gives the `x` and `y` values, as well as any genes in the region.<jupyter_code>gf.apply(stat_plot, stat="tajimaD", label="Tajima'sD")<jupyter_output><empty_output><jupyter_text>GC% of Plasmodium falciparium<jupyter_code>gf2.apply(plot_gc_data)<jupyter_output><empty_output><jupyter_text>GenomePlot classBlogs are supposed to be narrative driven, so I will try not to spend much time discussing the below `GenomePlot` class.The main gist is that we first initialise the class with information about the genome we are plotting:1. a fasta file, which provides information about contig names and sizes 2. a list of contigs contained in the above that we wish to display (default all)3. a layout of the above, describing how many contigs on each row etc.(and optionally)4. a function that is applied to each subplot to customise layout that's not data specificThe class contains an `apply()` method, that creates one subplot per contig, and places them in a grid. The function that actually makes the subplots is passed to `apply()` as an argument, and is called for each contig.<jupyter_code>class GenomePlot:
    
    genome = None
    contigs = None
    tools = "pan,wheel_zoom,box_zoom,save,reset"
    layout = None
    pfunc = None
    # we need these min border values to ensure that different contig sizes do not squeeze axes
    min_border_left = 50
    min_border_right = 10
    major_tick_dist = 1e7
    plot_width_per_mb = 6
    chrom_label_func = lambda self, y: y
    
    def __init__(self, fasta, contigs=None, layout=None, share_y=True, pfunc=None):
        
        self.genome = pyfasta.Fasta(fasta)
        
        if contigs is None:
            self.contigs = list(self.genome.keys())
        else:
            self.contigs = contigs
            
        # handle layout
        if layout is not None:
            self.layout = self.parse_layout(layout)
        else:
            self.layout = self.auto_layout()
        
        if pfunc is not None:
            self.pfunc = pfunc
                
    def parse_layout(self, lstring):
        layout = list()
        assert type(lstring) == str, "layout must be a string"
        q = list(self.contigs)
        q.reverse()
            
        rows = lstring.split("|")
        for row in rows:
            r = list()
            for x in row:
                assert x in "ox.", "Only chars o, x, . are allowed"
                if x == ".":
                    r.append(None)
                else:
                    r.append(q.pop())
            layout.append(r)
        return layout
    
    def auto_layout():
        pass

    def apply(self, func, **kwargs):
        
        # create a figure with specified layout        
        d = [[] for i in range(len(self.layout))]
        
        for i, row in enumerate(self.layout):
            
            is_bottom = ((i + 1) == len(self.layout))
            
            for j, chrom in enumerate(row):
                
                is_left = (j == 0)
                
                if chrom is not None:
                    
                    csize = len(self.genome[chrom])
                    px = int(csize * 1e-6 * self.plot_width_per_mb) 
                    px += self.min_border_left 
                    px += self.min_border_right
                    
                    try:
                        yrange = s1.y_range
                    except NameError:
                        yrange = None

                    s1 = figure(width=px, 
                                min_border_left=self.min_border_left,
                                min_border_right=self.min_border_right,
                                plot_height=250, 
                                tools=self.tools,
                                title=self.chrom_label_func(chrom),
                                y_range=yrange,
                                x_range=(1, csize))

                    s1.xaxis.ticker = FixedTicker(
                        ticks=np.arange(0, csize, self.major_tick_dist))
                    s1.xaxis[0].formatter = NumeralTickFormatter(format="0a.0")
                    
                    # handle general plot things specific to genome not data
                    if self.pfunc is not None:
                        self.pfunc(chrom, s1)

                    # function lives here
                    func(chrom, s1, is_left=is_left, is_bottom=is_bottom, **kwargs)
                    d[i].append(s1)

                else:
                    d[i].append(None)
            
            
        # put the subplots in a gridplot
        p = gridplot(d, toolbar_location="left", sizing_mode='fixed', plot_width=None)
        
        show(p)<jupyter_output><empty_output><jupyter_text>Example 1: Anopheles gambiaeHere, I want to plot some statistics across the [*Anopheles gambiae* genome](https://www.nature.com/articles/nature24995): Wattersons's $\Theta$, Tajima's D, and Nucleotide diversity ($\Pi$). These data are precomputed in another notebook - [here](https://github.com/hardingnj/hardingnj.github.io/blob/master/_posts/generate_stats.ipynb), using [scikit-allel](http://scikit-allel.readthedocs.io/en/latest/)<jupyter_code>%run generate_stats.ipynb<jupyter_output><empty_output><jupyter_text>All these data are contained in a `dict` containing `pandas` dataframes that look like this:<jupyter_code># where X is a chromosome label
annotated_data["X"].head()<jupyter_output><empty_output><jupyter_text>Define plotting functionsAs mentioned above, the `GenomePlot` class has an `.apply()` method, that takes a plotting function, and applies it to each contig.Plotting functions must take at least 4 arguments: the contig label (chrom), the subplot object, and two booleans describing whether the subplot is to be placed on the left and/or at the bottom of the grid, to allow customisation depending on position.This particular function has 5th and 6th arguments describing which statistic is to be plotted and the label. You'll notice that we also only make a Y-axis label on plots on the left of the grid, and a X-axis label if we are positioned on the bottom.We use the HoverTool class from bokeh to annotate points with chromosomal position, contig ID, and gene names.<jupyter_code>def stat_plot(chrom, subplot, is_left=None, is_bottom=None, stat="diversity", label=None):
    
    source = ColumnDataSource(annotated_data[chrom])
    
    if label is None:
        label = stat
    
    hover = HoverTool(
        tooltips=[
            ("Position", "@start{0a.000}-@stop{0a.000}"),
            (label, "$y"),
            ("contig", chrom), 
            ("gene(s)", "@gene")],
        mode="mouse")
    
    subplot.add_tools(hover)
    
    subplot.circle("midpoint", 
                   stat,
                   source=source,
                   size=3, 
                   color="navy", 
                   alpha=0.5)
    
    if is_left:
        subplot.yaxis[0].axis_label = label
    if is_bottom:
        subplot.xaxis[0].axis_label = "genome position (bp)"<jupyter_output><empty_output><jupyter_text>Custom layout functionAnopheles have 3 chromosomes (2, 3 and X). The autosomes are split into "left" and "right" arms (i.e. 2R, 2L, 3R, 3L). We use the custom layout function to move the Y-axis of the left arms to the right hand side.<jupyter_code># Custom changes to axes that pertain to genome not data
# eg in this case move left chromosome arm axes to RHS
def anopheles_plot(chrom, subplot):
    if chrom.endswith("L"):
        subplot.yaxis.visible = False
        subplot.add_layout(LinearAxis(), 'right')<jupyter_output><empty_output><jupyter_text>LayoutThe layout string below means that we place two plots on the top row, and 3 on the bottom. 'o' represents a plot, and '|' a line break. Gaps can be left using '.'A little bit of trial and error may be required to get this right!<jupyter_code>gf = GenomePlot(fasta=phase1_ar3.genome_fn,              # filepath to fasta
                contigs=("2R", "2L", "3R", "3L", "X"),   # contigs to display in order
                layout="oo|ooo",                         # layout string
                pfunc=anopheles_plot)                    # custom layout function<jupyter_output><empty_output><jupyter_text>Any further arguments to the plotting function can be passed as `**kwargs` to the `apply()` method. Tajima's D<jupyter_code>gf.apply(stat_plot, stat="tajimaD", label="Tajima'sD")<jupyter_output><empty_output><jupyter_text>Wattersons $\Theta$<jupyter_code>gf.apply(stat_plot, stat="theta", label="\u0398W")<jupyter_output><empty_output><jupyter_text>Changing the plotting function Here, instead of the circles shown previously, we present nucleotide diversity as a line, by changing the plotting function. Note we also change the `HoverTool` to use a vline, instead of a mouseover.<jupyter_code>def stat_plot_line(chrom, subplot, is_left=None, is_bottom=None, stat="diversity", label=None):
    
    source = ColumnDataSource(annotated_data[chrom])
    
    if label is None:
        label = stat
    
    hover = HoverTool(
        tooltips=[
            ("Position", "@start{0a.000}-@stop{0a.000}"),
            (label, "$y"),
            ("contig", chrom), 
            ("gene(s)", "@gene")],
        mode="vline")
    
    subplot.add_tools(hover)
    
    subplot.line("midpoint",
                 stat,
                 source=source,
                 color="navy",
                 alpha=0.5)

    if is_left:
        subplot.yaxis[0].axis_label = label
    if is_bottom:
        subplot.xaxis[0].axis_label = "position (bp)"<jupyter_output><empty_output><jupyter_text>The nucleotide diversity ($\Pi$) plot:<jupyter_code>gf.apply(stat_plot_line, stat="diversity", label="\u03A0")<jupyter_output><empty_output><jupyter_text>Example 2:  Plasmodium falciparumPl. falciparum is known to have a highly AT rich genome- here I am going to plot the % GC content in windows across the genome using the reference genome, downloaded from [PlasmoDB](http://plasmodb.org/plasmo/).<jupyter_code>plasmo_fn = "../assets/PlasmoDB-35_Pfalciparum3D7_Genome.fasta"
fa = pyfasta.Fasta(plasmo_fn, key_fn=lambda y: y.split(" | ")[0])
list(fa.keys())<jupyter_output><empty_output><jupyter_text>This function calculates the proportion of GC bases in windows.<jupyter_code>def calc_gc(genome, window_size=10000):
    
    clen = len(genome)
    starts = np.arange(0, clen, window_size)
    stops = starts + window_size
    
    propgc = np.zeros(starts.shape[0])
    
    for i, ix in enumerate(starts):
        r = genome[ix:(ix + window_size)]
        propgc[i] = (r.count("C") + r.count("G")) / (len(r) - r.count("N"))
        
    return pd.DataFrame.from_dict({"start": starts, "stop": stops, "gc": propgc * 100})
gc_data = {c: calc_gc(fa[c]) for c in fa.keys()}<jupyter_output><empty_output><jupyter_text>This is our plotting function, no extra arguments this time.Also we don't add a custom plot layout function in this case.<jupyter_code>def plot_gc_data(chrom, subplot, is_left=False, is_bottom=False):
    
    source = ColumnDataSource(gc_data[chrom])
    
    hover = HoverTool(
        tooltips=[
            ("Position", "@start{0a.00}-@stop{0a.00}"),
            ("%GC", "$y"),
            ("contig", chrom)],
        mode="vline")
    
    subplot.add_tools(hover)
    
    subplot.line("start",
                 "gc",
                 source=source,
                 color="navy", 
                 alpha=0.5)

    if is_bottom:
        subplot.xaxis[0].axis_label = "position (bp)"
    
    if is_left:
        subplot.yaxis[0].axis_label = "%GC"
clen = {c: len(fa[c]) for c in fa.keys()}
clen<jupyter_output><empty_output><jupyter_text>The GenomePlot class stuggles to present data where some chromosomes are vastly different sizes, so here we exclude the API and the mitochondrial chromosomes.<jupyter_code>plas_contigs = sorted(fa.keys())[:14]
gf2 = GenomePlot(plasmo_fn, contigs=plas_contigs, layout="ooooo|oooo|ooo|oo")

# set some parameters default work better with anopheles
gf2.plot_width_per_mb = 100
gf2.major_tick_dist = 5e5
gf2.chrom_label_func = lambda y: " ".join(y.split("_")[:2])
gf2.apply(plot_gc_data)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09575016572144067, 'fraction_numerical': 0.009354054651248434, 'mean_word_length': 2.592063492063492, 'pattern_counts': {'":': 4, '<': 69, '<?xml version=': 0, '>': 69, 'https://': 3, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '612775', 'license': 'MIT', 'n_tokens_mistral': 4286, 'n_tokens_neox': 3955, 'n_words': 1251, 'path': '_posts/2017-12-04-interactive-plots-with-bokeh.ipynb', 'repo_name': 'hardingnj/hardingnj.github.io', 'size': '13577'}