Fano

dynamo/configuration.py

@@ -9,7 +9,7 @@
 from cycler import cycler
 from matplotlib import cm, colors, rcParams
 
-from .dynamo_logger import main_info, main_warning
+from .dynamo_logger import main_debug, main_info
 
 
 class DynamoAdataKeyManager:
@@ -786,5 +786,5 @@ def set_pub_style_mpltex():
 
 # initialize DynamoSaveConfig and DynamoVisConfig mode defaults
 DynamoAdataConfig.update_data_store_mode("full")
-main_info("setting visualization default mode in dynamo. Your customized matplotlib settings might be overritten.")
+main_debug("setting visualization default mode in dynamo. Your customized matplotlib settings might be overwritten.")
 DynamoVisConfig.set_default_mode()

dynamo/dynamo_logger.py

@@ -154,8 +154,8 @@ def error(self, message, indent_level=1, *args, **kwargs):
 
     def info_insert_adata(self, key, adata_attr="obsm", indent_level=1, *args, **kwargs):
         message = "<insert> %s to %s in AnnData Object." % (key, adata_attr)
-        message = format_logging_message(message, logging.INFO, indent_level=indent_level)
-        return self.logger.error(message, *args, **kwargs)
+        message = format_logging_message(message, logging.DEBUG, indent_level=indent_level)
+        return self.logger.debug(message, *args, **kwargs)
 
     def info_insert_adata_var(self, key, indent_level=1, *args, **kwargs):
         return self.info_insert_adata(self, key, adata_attr="var", indent_level=1, *args, **kwargs)
@@ -189,18 +189,18 @@ def report_progress(self, percent=None, count=None, total=None, progress_name=""
 
     def finish_progress(self, progress_name="", time_unit="s", indent_level=1):
         self.log_time()
-        self.report_progress(percent=100, progress_name=progress_name)
+        #self.report_progress(percent=100, progress_name=progress_name)
 
         saved_terminator = self.logger_stream_handler.terminator
         self.logger_stream_handler.terminator = ""
-        self.logger.info("\n")
+        #self.logger.info("\n")
         self.logger_stream_handler.flush()
         self.logger_stream_handler.terminator = saved_terminator
 
         if time_unit == "s":
-            self.info("[%s] finished [%.4fs]" % (progress_name, self.time_passed), indent_level=indent_level)
+            self.info("[%s] completed [%.4fs]" % (progress_name, self.time_passed), indent_level=indent_level)
         elif time_unit == "ms":
-            self.info("[%s] finished [%.4fms]" % (progress_name, self.time_passed * 1e3), indent_level=indent_level)
+            self.info("[%s] completed [%.4fms]" % (progress_name, self.time_passed * 1e3), indent_level=indent_level)
         else:
             raise NotImplementedError
         # self.logger.info("|")

dynamo/external/sctransform.py

@@ -8,7 +8,6 @@
 # =================================================================
 
 import os
-from multiprocessing import Manager, Pool
 
 import numpy as np
 import pandas as pd
@@ -128,6 +127,8 @@ def sctransform_core(
     """
     A re-implementation of SCTransform from the Satija lab.
     """
+    import multiprocessing
+
     main_info("sctransform adata on layer: %s" % (layer))
     X = DKM.select_layer_data(adata, layer).copy()
     X = sp.sparse.csr_matrix(X)
@@ -139,10 +140,8 @@ def sctransform_core(
     genes_ix = genes.copy()
 
     X = X[:, genes]
-    Xraw = X.copy()
     gene_names = gene_names[genes]
     genes = np.arange(X.shape[1])
-    genes_cell_count = X.sum(0).A.flatten()
 
     genes_log_gmean = np.log10(gmean(X, axis=0, eps=gmean_eps))
 
@@ -188,7 +187,10 @@ def sctransform_core(
     bin_ind = np.ceil(np.arange(1, genes_step1.size + 1) / bin_size)
     max_bin = max(bin_ind)
 
-    ps = Manager().dict()
+    ps = multiprocessing.Manager().dict()
+
+    # create a process context of fork that copy a Python process from an existing process.
+    ctx = multiprocessing.get_context("fork")
 
     for i in range(1, int(max_bin) + 1):
         genes_bin_regress = genes_step1[bin_ind == i]
@@ -197,7 +199,9 @@ def sctransform_core(
         mm = np.vstack((np.ones(data_step1.shape[0]), data_step1["log_umi"].values.flatten())).T
 
         pc_chunksize = umi_bin.shape[1] // os.cpu_count() + 1
-        pool = Pool(os.cpu_count(), _parallel_init, [genes_bin_regress, umi_bin, gene_names, mm, ps])
+
+        pool = ctx.Pool(os.cpu_count(), _parallel_init, [genes_bin_regress, umi_bin, gene_names, mm, ps])
+
         try:
             pool.map(_parallel_wrapper, range(umi_bin.shape[1]), chunksize=pc_chunksize)
         finally:
@@ -254,10 +258,6 @@ def sctransform_core(
     full_model_pars["theta"] = theta
     del full_model_pars["dispersion"]
 
-    model_pars_outliers = outliers
-
-    regressor_data = np.vstack((np.ones(cell_attrs.shape[0]), cell_attrs["log_umi"].values)).T
-
     d = X.data
     x, y = X.nonzero()
     mud = np.exp(full_model_pars.values[:, 0][y] + full_model_pars.values[:, 1][y] * cell_attrs["log_umi"].values[x])

dynamo/plot/preprocess.py

@@ -9,7 +9,7 @@
 from ..configuration import DynamoAdataKeyManager
 from ..dynamo_logger import main_warning
 from ..preprocessing import preprocess as pp
-from ..preprocessing.preprocess_monocle_utils import top_table
+from ..preprocessing.gene_selection import get_prediction_by_svr
 from ..preprocessing.utils import detect_experiment_datatype
 from ..tools.utils import get_mapper, update_dict
 from .utils import save_fig
@@ -649,10 +649,9 @@ def feature_genes(
     save_show_or_return: str = "show",
     save_kwargs: dict = {},
 ):
-    """Plot selected feature genes on top of the mean vs. dispersion scatterplot.
+    """Plot selected feature genes on top of the mean vs. dispersion scatter plot.
 
-    Parameters
-    ----------
+    Args:
         adata: :class:`~anndata.AnnData`
             AnnData object
         layer: `str` (default: `X`)
@@ -664,32 +663,22 @@ def feature_genes(
         save_show_or_return: {'show', 'save', 'return'} (default: `show`)
             Whether to save, show or return the figure.
         save_kwargs: `dict` (default: `{}`)
-            A dictionary that will passed to the save_fig function. By default it is an empty dictionary and the
+            A dictionary that will be passed to the save_fig function. By default, it is an empty dictionary and the
             save_fig function will use the {"path": None, "prefix": 'feature_genes', "dpi": None, "ext": 'pdf',
-            "transparent": True, "close": True, "verbose": True} as its parameters. Otherwise you can provide a
+            "transparent": True, "close": True, "verbose": True} as its parameters. Otherwise, you can provide a
             dictionary that properly modify those keys according to your needs.
 
-    Returns
-    -------
+    Returns:
         Nothing but plots the selected feature genes via the mean, CV plot.
     """
 
     import matplotlib.pyplot as plt
 
     mode = adata.uns["feature_selection"] if mode is None else mode
-
     layer = DynamoAdataKeyManager.get_available_layer_keys(adata, layer, include_protein=False)[0]
-
     uns_store_key = None
-    if mode == "dispersion":
-        uns_store_key = "dispFitInfo" if layer in ["raw", "X"] else layer + "_dispFitInfo"
 
-        table = top_table(adata, layer)
-        x_min, x_max = (
-            np.nanmin(table["mean_expression"]),
-            np.nanmax(table["mean_expression"]),
-        )
-    elif mode == "SVR":
+    if "_dispersion" in mode:  # "cv_dispersion", "fano_dispersion"
         prefix = "" if layer == "X" else layer + "_"
         uns_store_key = "velocyto_SVR" if layer == "raw" or layer == "X" else layer + "_velocyto_SVR"
 
@@ -709,11 +698,12 @@ def feature_genes(
     ordering_genes = adata.var["use_for_pca"] if "use_for_pca" in adata.var.columns else None
 
     mu_linspace = np.linspace(x_min, x_max, num=1000)
-    fit = (
-        adata.uns[uns_store_key]["disp_func"](mu_linspace)
-        if mode == "dispersion"
-        else adata.uns[uns_store_key]["SVR"](mu_linspace.reshape(-1, 1))
-    )
+    if "_dispersion" in mode:
+        mean = adata.uns[uns_store_key]["mean"]
+        cv = adata.uns[uns_store_key]["cv"]
+        svr_gamma = adata.uns[uns_store_key]["svr_gamma"]
+        fit, _ = get_prediction_by_svr(mean, cv, svr_gamma)
+        fit = fit(mu_linspace.reshape(-1, 1))
 
     plt.figure(figsize=figsize)
     plt.plot(mu_linspace, fit, alpha=0.4, color="r")
@@ -724,15 +714,7 @@ def feature_genes(
     )
 
     valid_disp_table = table.iloc[valid_ind, :]
-    if mode == "dispersion":
-        ax = plt.scatter(
-            valid_disp_table["mean_expression"],
-            valid_disp_table["dispersion_empirical"],
-            s=3,
-            alpha=1,
-            color="xkcd:red",
-        )
-    elif mode == "SVR":
+    if "_dispersion" in mode:
         ax = plt.scatter(
             valid_disp_table[prefix + "log_m"],
             valid_disp_table[prefix + "log_cv"],
@@ -743,15 +725,7 @@ def feature_genes(
 
     neg_disp_table = table.iloc[~valid_ind, :]
 
-    if mode == "dispersion":
-        ax = plt.scatter(
-            neg_disp_table["mean_expression"],
-            neg_disp_table["dispersion_empirical"],
-            s=3,
-            alpha=0.5,
-            color="xkcd:grey",
-        )
-    elif mode == "SVR":
+    if "_dispersion" in mode:
         ax = plt.scatter(
             neg_disp_table[prefix + "log_m"],
             neg_disp_table[prefix + "log_cv"],
@@ -760,9 +734,6 @@ def feature_genes(
             color="xkcd:grey",
         )
 
-    # plt.xlim((0, 100))
-    if mode == "dispersion":
-        plt.xscale("log")
     plt.yscale("log")
     plt.xlabel("Mean (log)")
     plt.ylabel("Dispersion (log)") if mode == "dispersion" else plt.ylabel("CV (log)")
@@ -1080,7 +1051,7 @@ def highest_frac_genes(
 
         else:
             main_warning(
-                "%s not in adata.var, ignoring the gene annotation key when plotting",
+                "%s not in adata.var, ignoring the gene annotation key when plotting" % gene_annotation_key,
                 indent_level=2,
             )