first implementation of annotation track when a gene coordinate csv i…

…s provided

tsqc/__main__.py

@@ -14,6 +14,7 @@
 
 from . import model  # noqa
 from . import pages  # noqa
+from . import config  # noqa
 
 logger = daiquiri.getLogger("tsqc")
 
@@ -116,6 +117,7 @@ def setup_logging(log_level, no_log_filter):
 
 @click.command()
 @click.argument("path", type=click.Path(exists=True, dir_okay=False))
+@click.option("--annotate-genes", type=click.Path(exists=True, dir_okay=False))
 @click.option("--port", default=8080, help="Port to serve on")
 @click.option(
     "--show/--no-show",
@@ -129,12 +131,15 @@ def setup_logging(log_level, no_log_filter):
     is_flag=True,
     help="Do not filter the output log (advanced debugging only)",
 )
-def main(path, port, show, log_level, no_log_filter):
+def main(path, port, show, log_level, no_log_filter, annotate_genes):
     """
     Run the tsqc server.
     """
     setup_logging(log_level, no_log_filter)
+
     tsm = load_data(pathlib.Path(path))
+    if annotate_genes:
+        config.GENES_FILE = annotate_genes
 
     # Note: functools.partial doesn't work here
     def app():

tsqc/config.py

@@ -2,3 +2,4 @@
 PLOT_WIDTH = 1000
 PLOT_HEIGHT = 600
 THRESHOLD = 1000  # max number of points to overlay on a plot
+GENES_FILE = None  # set through command line option

tsqc/model.py

@@ -501,6 +501,17 @@ def trees_df(self):
             }
         )
 
+    def genes_df(self, genes_file):
+        genes_df = pd.read_csv(genes_file, sep=";")
+        # TODO file checks!
+        genes_df.columns = ["chr", "position", "end", "strand", "id", "name"]
+        genes_df = genes_df[
+            (genes_df["position"] >= self.ts.first().interval.left)
+            & (genes_df["end"] <= self.ts.last().interval.right)
+        ]
+        logger.info("Computed genes dataframe")
+        return genes_df
+
     def calc_polytomy_fractions(self):
         """
         Calculates the fraction of polytomies for each tree in the

tsqc/pages/mutations.py

@@ -1,6 +1,5 @@
 import holoviews as hv
 import holoviews.operation.datashader as hd
-import hvplot.pandas  # noqa
 import numpy as np
 import pandas as pd
 import panel as pn
@@ -57,7 +56,7 @@ def make_muts_panel(log_y, tsm):
     shaded = hd.datashade(filtered, width=400, height=400, streams=streams)
 
     main = (shaded * hover).opts(
-        hv.opts.Points(tools=["hover"], alpha=0.1, hover_alpha=0.2, size=10)
+        hv.opts.Points(tools=["hover"], alpha=0.1, hover_alpha=0.2, size=10),
     )
 
     time_hist = hv.DynamicMap(
@@ -78,25 +77,76 @@ def make_muts_panel(log_y, tsm):
         }
     )
     trees_hist = hv.DynamicMap(selected_hist(bp_df), streams=streams)
+    trees_hist.opts(
+        width=config.PLOT_WIDTH,
+        height=100,
+        hooks=[customise_ticks],
+        xlabel="tree density",
+    )
     trees_line = hv.Segments(bp_df, ["position", "y0", "x1", "y1"])
+    trees_line.opts(
+        width=config.PLOT_WIDTH,
+        height=100,
+        yaxis=None,
+        ylabel=None,
+        xlabel="tree breakpoints",
+        color="grey",
+        line_width=0.5,
+        alpha=0.5,
+    )
+
+    if config.GENES_FILE:
+        genes_df = tsm.genes_df(config.GENES_FILE)
+        annot_track = make_annotation_plot(tsm, genes_df)
+        layout = (
+            ((main << time_hist << site_hist) + trees_hist + trees_line + annot_track)
+            .opts(shared_axes=True)
+            .cols(1)
+        )
 
-    return pn.Column(
-        (
-            (main << time_hist << site_hist)
-            + trees_hist.opts(
-                width=config.PLOT_WIDTH, height=100, hooks=[customise_ticks]
-            )
-            + trees_line.opts(
-                width=config.PLOT_WIDTH,
-                height=100,
-                yaxis=None,
-                ylabel=None,
-                xlabel="Tree breakpoints",
-            )
+    else:
+        layout = (
+            ((main << time_hist << site_hist) + trees_hist + trees_line)
+            .opts(shared_axes=True)
+            .cols(1)
         )
-        .opts(shared_axes=True)
-        .cols(1)
+
+    return pn.Column(layout)
+
+
+def make_annotation_plot(tsm, genes_df):
+    min_y = tsm.mutations_df["time"].min()
+    max_y = tsm.mutations_df["time"].max()
+    genes_df["y0"] = min_y + 0.3 * (max_y - min_y)
+    genes_df["y1"] = max_y - 0.3 * (max_y - min_y)
+    genes_rects = hv.Rectangles(
+        genes_df, kdims=["position", "y0", "end", "y1"], vdims=["name", "id", "strand"]
+    )
+    hover_tool = HoverTool(
+        tooltips=[
+            ("gene name", "@name"),
+            ("ensembl id", "@id"),
+            ("strand", "@strand"),
+        ]
+    )
+    genes_rects.opts(
+        ylabel=None,
+        line_color=None,
+        shared_axes=True,
+        color="maroon",
+        hooks=[customise_ticks],
+        width=config.PLOT_WIDTH,
+        height=100,
+        yaxis=None,
+        xlabel="genes",
+        tools=[hover_tool],
+    )
+
+    genes_rects = (
+        hv.HLine(min_y + (max_y - min_y) / 2).opts(color="black", line_width=0.7)
+        * genes_rects
     )
+    return genes_rects
 
 
 def page(tsm):