Merge branch 'release/0.3.4'

scepia/cli.py

@@ -0,0 +1,53 @@
+#!/usr/bin/env python
+import click
+import scepia
+import sys
+import os
+from scepia import _version
+from scepia import generate_signal, link_it_up
+from scepia.sc import full_analysis
+
+CONTEXT_SETTINGS = dict(help_option_names=["-h", "--help"])
+
+@click.group(context_settings=CONTEXT_SETTINGS)
+@click.version_option(_version)
+def cli():
+    """SCEPIA - Single Cell Epigenome-based Inference of Activity
+
+    Version: {}""".format(
+        _version
+    )
+    pass
+
+
+@click.command("area27", short_help="Determine the enhancer-based regulatory potential (ERP) score.")
+@click.argument("bamfile")
+@click.argument("outfile")
+@click.option("-w", "--window", help="window")
+@click.option("-N", "--nthreads", help="number of threads")
+def area27(bamfile, outfile, window=2000, nthreads=4):
+    """
+    Determine the enhancer-based regulatory potential (ERP) score per gene. This
+    approach is based on the method developed by Wang et al., 2016. There is one 
+    difference. In this implementation the score is calculated based only on 
+    H3K27ac signal in enhancers. We use log-transformed, z-score normalized 
+    H3K27ac read counts in 2kb windows centered at enhancer locations. 
+    """
+    signal = generate_signal(bamfile, window=2000, nthreads=nthreads)
+    link_it_up(outfile, signal)
+
+@click.command("infer_motifs", short_help="Run SCEPIA motif inference on an h5ad file.")
+@click.argument("infile")
+@click.argument("outfile")
+def infer_motifs(infile, outfile):
+    """
+    Infer motifs.
+    """
+    full_analysis(infile, outfile)   
+
+
+cli.add_command(area27)
+cli.add_command(infer_motifs)
+
+if __name__ == "__main__":
+    cli()

scepia/sc.py

@@ -302,7 +302,7 @@ def annotate_with_k27(
             if use_neighbors:
                 my_neighbors = (
                     pd.DataFrame(
-                        (adata.uns["neighbors"]["connectivities"][i] != 0).todense()
+                        (adata.obsp["connectivities"][i] != 0).todense()
                     )
                     .iloc[0]
                     .values
@@ -559,16 +559,18 @@ def infer_motifs(
         run_maelstrom(
             fname,
             "hg38",
-            "tmp.lala",
-            methods=["bayesianridge", "lightningregressor", "xgboost"],
+            "scepia.maelstrom",
+            center=True,
+            filter_redundant=False,
         )
 
         motif_act = pd.read_csv(
-            os.path.join("tmp.lala", "final.out.csv"),
+            os.path.join("scepia.maelstrom", "final.out.txt"),
             sep="\t",
             comment="#",
             index_col=0,
         )
+        motif_act.columns = motif_act.columns.str.replace("z-score\s+", "")
     else:
         motif_act = moap(
             fname,
@@ -956,9 +958,9 @@ def plot_volcano_corr(
     n_anno: Optional[int] = 40,
     size_anno: Optional[float] = 6,
     palette: Optional[str] = None,
-    alpha: Optional[float] = 0.6,
+    alpha: Optional[float] = 0.8,
     linewidth: Optional[float] = 0,
-    sizes: Optional[Tuple[int, int]] = (1, 20),
+    sizes: Optional[Tuple[int, int]] = (5, 25),
     **kwargs,
 ) -> Axes:
     """Volcano plot of significance of motif-TF correlations.
@@ -975,15 +977,15 @@ def plot_volcano_corr(
             "Motif-TF correlation data not found. Did you run `determine_significance()`?"
         )
 
-    if palette is None:
-        n_colors = len(
-            sns.utils.categorical_order(adata.uns["scepia"]["correlation"]["std"])
-        )
-        cmap = LinearSegmentedColormap.from_list(
-            name="grey_black", colors=["grey", "black"]
-        )
-        palette = sns.color_palette([cmap(i) for i in np.arange(0, 1, 1 / n_colors)])
-
+#    if palette is None:
+#        n_colors = len(
+#            sns.utils.categorical_order(adata.uns["scepia"]["correlation"]["std"])
+#        )
+#        cmap = LinearSegmentedColormap.from_list(
+#            name="grey_black", colors=["grey", "black"]
+#        )
+#        palette = sns.color_palette([cmap(i) for i in np.arange(0, 1, 1 / n_colors)])
+#
     sns.set_style("ticks")
     g = sns.scatterplot(
         data=adata.uns["scepia"]["correlation"],
@@ -1042,3 +1044,13 @@ def test(n_rounds=100):
 
     pool.close()
     return result
+
+
+def full_analysis(infile, outfile):
+    """
+    Run full SCEPIA analysis on h5ad infile.
+    """
+    adata = sc.read(infile)
+    adata = infer_motifs(adata, dataset="ENCODE")
+    determine_significance(adata)
+    adata.write(outfile)

setup.py

@@ -7,6 +7,8 @@
 with open("README.md") as f:
     long_description = f.read()
 
+entry_points = {"console_scripts": ["scepia=scepia.cli:cli"]}
+
 setup(
     name="scepia",
     version=versioneer.get_version(),
@@ -18,7 +20,7 @@
     url="https://github.com/vanheeringen-lab/scepia/",
     license="MIT",
     packages=find_packages(),
-    scripts=["scripts/scepia"],
+    entry_points=entry_points,
     include_package_data=True,
     zip_safe=False,
     classifiers=[