Merge pull request #10 from Mye-InfoBank/serafina/scvi_tools_dgea

Serafina/scvi tools dgea

Author: nictru

app.py

@@ -5,14 +5,15 @@
 
 from composition import composition_server, composition_ui
 from export import export_ui, export_server
-from dgea.dgea import dgea_server, dgea_ui
+from dgea.dgea_scvi import dgea_server, dgea_ui
 from tree import tree_server, tree_ui
 
 with open("data/config.json") as f:
     config = json.load(f)
     adata = sc.read_h5ad("data/" + config["adata"])
     tree = pickle.load(open("data/" + config["tree"], "rb")) if "tree" in config else None
     name = config["name"]
+    model_path = config["model_path"]
 
 categorical_columns = adata.obs.select_dtypes(include="category").columns.to_list()
 
@@ -34,6 +35,7 @@ def server(input, output, session):
     _dataframe = reactive.value(adata.obs)
     _adata = reactive.value(adata)
     _tree = reactive.value(tree)
+    _model = reactive.value(model_path)
     composition_server("composition", _dataframe)
     export_server("export")
     dgea_server("dgea", _adata)

data/atlas.h5ad

@@ -1,5 +1,6 @@
 {
     "adata": "merged.h5ad",
     "tree": "scarches.tree.pkl",
-    "name": "possible_atlas"
+    "name": "possible_atlas",
+    "model_path": "./data"
 }

data/config.json

@@ -0,0 +1,40 @@
+from shiny import reactive, ui, render, module
+import anndata as ad
+
+from dgea.run_dgea_scvi import run_dgea_ui, run_dgea_server
+from dgea.filter_dgea_scvi import filter_dgea_ui, filter_dgea_server
+from dgea.plot_dgea_scvi import plot_dgea_ui, plot_dgea_server
+
+@module.ui
+def dgea_ui():
+    return ui.layout_sidebar(
+                     ui.sidebar(
+                        run_dgea_ui("run_dgea"),
+                        filter_dgea_ui("filter_dgea"),
+                        title="Select covariates"
+                     ),
+                     *plot_dgea_ui("plot_dgea")
+                    )
+
+@module.server
+def dgea_server(input, output, session, _adata: reactive.Value[ad.AnnData]):
+    _counts = reactive.value(None)
+    _uniques = reactive.value([])
+
+    _contrast = reactive.value(None)
+    _reference = reactive.value(None)
+    _alternative = reactive.value(None)
+    _log10_p = reactive.value(0.05)
+    _lfc = reactive.value(1)
+
+    _result = reactive.value(None)
+    _filtered_result = reactive.value(None)
+    _filtered_genes = reactive.value(None)
+    _filtered_counts = reactive.value(None)
+
+    run_dgea_server("run_dgea", _adata, _result, _counts, _reference, _alternative, _uniques, _contrast)
+    filter_dgea_server("filter_dgea", _adata, _counts, _uniques,
+                       _result, _filtered_result, _filtered_genes,
+                       _filtered_counts, _reference, _alternative, _contrast, _log10_p, _lfc)
+    plot_dgea_server("plot_dgea", _filtered_counts, _contrast, _reference, _alternative,
+                     _result, _log10_p, _lfc)

data/model.pt

@@ -0,0 +1,50 @@
+import numpy as np
+import pandas as pd
+import scanpy as sc
+import anndata as ad
+
+from scvi.model import SCANVI, SCVI
+
+def scanvi_dgea(adata:ad.AnnData, groupby:str, reference:str, alternative:str, directory_model:str):
+          
+    if 'cell_type' in adata.obs.columns:
+        model_type = SCANVI
+        print('is scavi')
+        
+    else:
+        model_type = SCVI
+        print('is scanvi')
+  
+    model_type.prepare_query_anndata(adata = adata, reference_model=directory_model)
+    
+    model = model_type.load_query_data(adata, directory_model)
+    
+    groups = np.array(adata.obs[groupby].unique())
+
+    idx1 = adata.obs[groupby] == reference
+    idx2 = adata.obs[groupby] == alternative
+    
+    dge_change = model.differential_expression(adata=adata, groupby=groupby, idx1=idx1, idx2=idx2, mode="change")
+    
+    epsilon = 1e-10
+    dge_change['proba_not_de'] = np.maximum(dge_change["proba_not_de"], epsilon)
+    dge_change["log10_pscore"] = np.log10(dge_change["proba_not_de"])
+    dge_change["-log10_pscore"] = -np.log10(dge_change["proba_not_de"])	
+    
+    return dge_change
+
+def get_normalized_counts(adata):
+    print(adata.shape)
+    sc.pp.normalize_total(adata, target_sum=1e4)
+    sc.pp.log1p(adata)
+    adata.layers["counts"] = adata.X.copy().tocsr()
+    counts = adata.layers["counts"]
+    dense_matrix = counts.toarray()
+    df_counts = pd.DataFrame(dense_matrix, index=adata.obs_names, columns=adata.var_names)
+    return df_counts
+
+if __name__ == '__main__':
+    print('Running DGEA test')
+    adata = sc.read_h5ad('/workspaces/SIMBA-Downstream_1/data/atlas.h5ad')
+    dge_test = scanvi_dgea(adata, "cell_type", "Endothelial", "Epithelial", './data')
+    print(dge_test.head())

dgea/dgea_scvi.py

@@ -0,0 +1,95 @@
+from shiny import reactive, ui, render, module
+import anndata as ad
+import pandas as pd
+from dgea.dgea_scvi_helpers import scanvi_dgea, get_normalized_counts
+
+@module.ui
+def filter_dgea_ui():
+    return ui.div(
+        ui.output_ui("select_reference"),
+        ui.output_ui("select_alternative"),
+        ui.input_slider("log10_pscore", "Ropability in Reference (significance threshold)", min=0, max=20, step=0.01, value=3),
+        ui.input_slider("lfc", "Log2 fold change", min=0, max=10, step=0.1, value=1),
+        ui.output_ui("open_gprofiler")
+    )
+
+@module.server
+def filter_dgea_server(input, output, session, 
+                       _adata, _counts, _uniques,
+                       _result, _filtered_result, _filtered_genes, _filtered_counts,
+                       _reference, _alternative, _contrast,
+                       _log10_p, _lfc
+                       ):
+
+    @output
+    @render.ui
+    def select_reference():
+        uniques = _uniques.get()
+
+        if not uniques or len(uniques) < 2:
+            return ui.p("Run analysis to see options")
+
+        return ui.input_select("reference", "Reference", choices=uniques, selected=uniques[0])
+    
+    @output
+    @render.ui
+    def select_alternative():
+        uniques = _uniques.get()
+
+        if not uniques or len(uniques) < 2:
+            return ui.p("Run analysis to see options")
+        
+        print(uniques)
+
+        return ui.input_select("alternative", "Alternative", choices=uniques, selected=uniques[1])
+
+
+    @reactive.effect
+    def update_filters():
+        _reference.set(input["reference"].get())
+        _alternative.set(input["alternative"].get())
+        _log10_p.set(input["log10_pscore"].get())
+        _lfc.set(input["lfc"].get())
+
+    @reactive.effect
+    def update_result():
+        adata = _adata.get()
+        reference = _reference.get()
+        alternative = _alternative.get()
+        contrast = _contrast.get()
+
+        if None in (reference, alternative, contrast):
+            return
+
+        res_df = scanvi_dgea(adata, contrast, reference, alternative)
+        res_counts = get_normalized_counts(adata)
+        _result.set(res_df)
+        _counts.set(res_counts)
+    
+    @reactive.effect
+    def filter_result():
+        result = _result.get()
+        log10_p = input["log10_pscore"].get()
+        lfc = input["lfc"].get()
+        counts = _counts.get()
+
+        if result is None:
+            return None
+
+        result = result[(result["-log10_pscore"] < log10_p) & (result["lfc_mean"].abs() > lfc)]
+        _filtered_result.set(result)
+        genes = result.index.tolist()
+        genes_not_found = [gene for gene in genes if gene not in counts.columns]
+        if genes_not_found:
+            print(f"Genes not found in the DataFrame: {genes_not_found}")
+        _filtered_genes.set(genes)
+        _filtered_counts.set(counts.loc[:, genes])
+
+    @render.ui
+    def open_gprofiler():
+        genes = _filtered_genes.get()
+        if not genes:
+            return None
+        
+        return ui.input_action_button("gprofiler", label="Open g:Profiler",
+                              onclick=f"window.open('https://biit.cs.ut.ee/gprofiler/gost?organism=hsapiens&query={'%0A'.join(genes)}', '_blank')")

dgea/dgea_scvi_helpers.py

@@ -0,0 +1,116 @@
+from shiny import reactive, ui, render, module
+import numpy as np
+import seaborn as sns
+import shinywidgets as sw
+import tempfile
+import plotly.express as px
+
+@module.ui
+def plot_dgea_ui():
+    return [
+        ui.card(
+            ui.card_header("Volcano plot"),
+            sw.output_widget("plot_volcano")
+        ),
+        ui.card(
+            ui.card_header("Heatmap"),
+            ui.output_plot("plot_heatmap"),
+            ui.card_footer(
+                ui.download_button("download_dgea", "Download DGEA matrix"),
+                ui.download_button("download_plot", "Download plot")
+            )
+        )
+    ]
+
+@module.server
+def plot_dgea_server(input, output, session,
+                     _filtered_counts,
+                     _contrast,
+                     _reference,
+                     _alternative,
+                     _result,
+                     _log10_p,
+                     _lfc
+                     ):
+    _heatmap = reactive.value(None)
+
+    @output
+    @render.plot
+    def plot_heatmap():
+        counts_df = _filtered_counts.get()
+        contrast = _contrast.get()
+        reference = _reference.get()
+        alternative = _alternative.get()
+
+        if counts_df is None:
+            return None
+
+        if counts_df.empty:
+            return None
+
+        plot = sns.clustermap(counts_df.T, cmap="viridis", figsize=(10, 10))
+        _heatmap.set(plot)
+
+        return plot
+
+    @render.download(
+            filename=lambda: f"dgea_matrix_{_contrast.get()}-{_reference.get()}:{_alternative.get()}.csv"
+    )
+    def download_dgea():
+        scanvi_results = _result.get()
+        if scanvi_results is None:
+            return None
+        with tempfile.NamedTemporaryFile(suffix=".csv", delete=False) as temp:
+            scanvi_results.to_csv(temp.name)
+            return temp.name
+        
+    @render.download(
+            filename=lambda: f"heatmap_{_contrast.get()}-{_reference.get()}:{_alternative.get()}.png"
+    )
+    def download_plot():
+        plot = _heatmap.get()
+        if plot is None:
+            return None
+        with tempfile.NamedTemporaryFile(suffix=".png", delete=False) as temp:
+            plot.savefig(temp.name, bbox_inches="tight")
+            return temp.name
+    
+    @output
+    @sw.render_plotly
+    def plot_volcano():
+        scanvi_results = _result.get()
+        log10_p = _log10_p.get()
+        lfc = _lfc.get()
+
+        if scanvi_results is None:
+            return None
+        
+        df_plot = scanvi_results.copy()
+
+        df_plot["category"] = "Not significant"
+        df_plot["gene"] = df_plot.index
+        df_plot.loc[(df_plot["lfc_mean"] < -lfc) & (df_plot["-log10_pscore"] < log10_p), "category"] = f"High in {_alternative.get()}"
+        df_plot.loc[(df_plot["lfc_mean"] > lfc) & (df_plot["-log10_pscore"] < log10_p), "category"] = f"High in {_reference.get()}"
+
+        colormap = {"Not significant": "grey", f"High in {_alternative.get()}": "blue", f"High in {_reference.get()}": "red"}
+        
+        hover_data = {
+            "lfc_mean": True,
+            "-log10_pscore": True,
+            "category": True,
+        }
+        
+        fig = px.scatter(df_plot, x="lfc_mean", y="-log10_pscore", color="category",
+                          color_discrete_map=colormap,
+                          hover_name="gene",
+                          hover_data=hover_data,
+                          labels={"lfc_mean": "Log2 fold change mean",
+                                  "-log10_pscore": "Negative log 10 P-value",
+                                  "category": "Category"})
+        
+        num_points = len(df_plot)
+        fig.add_annotation(text=f"number of counts = {num_points}", x=0.5, y=-0.15, showarrow=False, font=dict(size=12), xanchor="center")
+        fig.update_layout(margin=dict(l=40, r=40, b=80, t=40), height=600)
+        
+        return fig
+