Source code for dynamo.preprocessing.utils

import warnings
from typing import Callable, Iterable, List, Tuple, Union

try:
    from typing import Literal
except ImportError:
    from typing_extensions import Literal

import anndata
import numpy as np
import numpy.typing as npt
import pandas as pd
import scipy
import scipy.sparse
import statsmodels.api as sm
from anndata import AnnData
from scipy.sparse import csr_matrix, issparse
from sklearn.decomposition import PCA, TruncatedSVD

from ..configuration import DKM, DynamoAdataKeyManager
from ..dynamo_logger import (
    LoggerManager,
    main_debug,
    main_exception,
    main_info,
    main_info_insert_adata_var,
    main_warning,
)
from ..utils import areinstance


# ---------------------------------------------------------------------------------------------------
# symbol conversion related
def convert2gene_symbol(input_names: List[str], scopes: Union[List[str], None] = "ensembl.gene") -> pd.DataFrame:
    """Convert ensemble gene id to official gene names using mygene package.

    Args:
        input_names: the ensemble gene id names that you want to convert to official gene names. All names should come
            from the same species.
        scopes: scopes are needed when you use non-official gene name as your gene indices (or adata.var_name).
            This arugument corresponds to type of types of identifiers, either a list or a comma-separated fields to
            specify type of input qterms, e.g. “entrezgene”, “entrezgene,symbol”, [“ensemblgene”, “symbol”]. Refer to
            official MyGene.info docs (https://docs.mygene.info/en/latest/doc/query_service.html#available_fields) for
            full list of fields. Defaults to "ensembl.gene".

    Raises:
        ImportError: fail to import `mygene`.

    Returns:
        A pandas dataframe that includes the following columns:
            query: the input ensmble ids
            _id: identified id from mygene
            _score: confidence of the retrieved official gene name.
            symbol: retrieved official gene name
    """

    try:
        import mygene
    except ImportError:
        raise ImportError(
            "You need to install the package `mygene` (pip install mygene --user) "
            "See https://pypi.org/project/mygene/ for more details."
        )

    mg = mygene.MyGeneInfo()
    main_info("Storing myGene name info into local cache db: mygene_cache.sqlite.")
    mg.set_caching()

    ensemble_names = [i.split(".")[0] for i in input_names]
    var_pd = mg.querymany(
        ensemble_names,
        scopes=scopes,
        fields="symbol",
        as_dataframe=True,
        df_index=True,
    )
    # var_pd.drop_duplicates(subset='query', inplace=True) # use when df_index is not True
    var_pd = var_pd.loc[~var_pd.index.duplicated(keep="first")]

    return var_pd


[docs]def convert2symbol(adata: AnnData, scopes: Union[str, Iterable, None] = None, subset=True) -> AnnData: """This helper function converts unofficial gene names to official gene names. Args: adata: an AnnData object. scopes: scopes are needed when you use non-official gene name as your gene indices (or adata.var_name). This arugument corresponds to type of types of identifiers, either a list or a comma-separated fields to specify type of input qterms, e.g. “entrezgene”, “entrezgene,symbol”, [“ensemblgene”, “symbol”]. Refer to official MyGene.info docs (https://docs.mygene.info/en/latest/doc/query_service.html#available_fields) for full list of fields. Defaults to None. subset: whether to inplace subset the results. Defaults to True. Raises: Exception: gene names in adata.var_names are invalid. Returns: The updated AnnData object. """ if np.all(adata.var_names.str.startswith("ENS")) or scopes is not None: logger = LoggerManager.gen_logger("dynamo-utils") logger.info("convert ensemble name to official gene name", indent_level=1) prefix = adata.var_names[0] if scopes is None: if prefix[:4] == "ENSG" or prefix[:7] == "ENSMUSG": scopes = "ensembl.gene" elif prefix[:4] == "ENST" or prefix[:7] == "ENSMUST": scopes = "ensembl.transcript" else: raise Exception( "Your adata object uses non-official gene names as gene index. \n" "Dynamo finds those IDs are neither from ensembl.gene or ensembl.transcript and thus cannot " "convert them automatically. \n" "Please pass the correct scopes or first convert the ensemble ID to gene short name " "(for example, using mygene package). \n" "See also dyn.pp.convert2gene_symbol" ) adata.var["query"] = [i.split(".")[0] for i in adata.var.index] if scopes is str: adata.var[scopes] = adata.var.index else: adata.var["scopes"] = adata.var.index logger.warning( "Your adata object uses non-official gene names as gene index. \n" "Dynamo is converting those names to official gene names." ) official_gene_df = convert2gene_symbol(adata.var_names, scopes) merge_df = adata.var.merge(official_gene_df, left_on="query", right_on="query", how="left").set_index( adata.var.index ) adata.var = merge_df valid_ind = np.where(merge_df["notfound"] != True)[0] # noqa: E712 if subset is True: adata._inplace_subset_var(valid_ind) adata.var.index = adata.var["symbol"].values.copy() else: indices = np.array(adata.var.index) indices[valid_ind] = adata.var.loc[valid_ind, "symbol"].values.copy() adata.var.index = indices if np.sum(adata.var_names.isnull()) > 0: main_info( "Subsetting adata object and removing Nan columns from adata when converting gene names.", indent_level=1, ) adata._inplace_subset_var(adata.var_names.notnull()) return adata
def compute_gene_exp_fraction(X: scipy.sparse.spmatrix, threshold: float = 0.001) -> Tuple[np.ndarray, np.ndarray]: """Calculate fraction of each gene's count to total counts across cells and identify high fraction genes. Args: X: a sparse matrix containing gene expression data. threshold: lower bound for valid data. Defaults to 0.001. Returns: A tuple (frac, valid_ids) where frac is the fraction of each gene's count to total count and valid_ids is the indices of valid genes. """ frac = X.sum(0) / X.sum() if issparse(X): frac = frac.A.reshape(-1, 1) valid_ids = np.where(frac > threshold)[0] return frac, valid_ids # --------------------------------------------------------------------------------------------------- # implmentation of Cooks' distance (but this is for Poisson distribution fitting) # https://stackoverflow.com/questions/47686227/poisson-regression-in-statsmodels-and-r # from __future__ import division, print_function # https://stats.stackexchange.com/questions/356053/the-identity-link-function-does-not-respect-the-domain-of-the-gamma- # family def _weight_matrix(fitted_model: sm.Poisson) -> np.ndarray: """Calculates weight matrix in Poisson regression. Args: fitted_model: a fitted Poisson model Returns: A diagonal weight matrix in Poisson regression. """ return np.diag(fitted_model.fittedvalues) def _hessian(X: np.ndarray, W: np.ndarray) -> np.ndarray: """Hessian matrix calculated as -X'*W*X. Args: X: the matrix of covariates. W: the weight matrix. Returns: The result Hessian matrix. """ return -np.dot(X.T, np.dot(W, X)) def _hat_matrix(X: np.ndarray, W: np.ndarray) -> np.ndarray: """Calculate hat matrix = W^(1/2) * X * (X'*W*X)^(-1) * X'*W^(1/2) Args: X: the matrix of covariates. W: the diagonal weight matrix Returns: The result hat matrix """ # W^(1/2) Wsqrt = W ** (0.5) # (X'*W*X)^(-1) XtWX = -_hessian(X=X, W=W) XtWX_inv = np.linalg.inv(XtWX) # W^(1/2)*X WsqrtX = np.dot(Wsqrt, X) # X'*W^(1/2) XtWsqrt = np.dot(X.T, Wsqrt) return np.dot(WsqrtX, np.dot(XtWX_inv, XtWsqrt)) def cook_dist(model: sm.Poisson, X: np.ndarray, good: npt.ArrayLike) -> np.ndarray: """calculate Cook's distance Args: model: a fitted Poisson model. X: the matrix of covariates. good: the dispersion table for MSE calculation. Returns: The result Cook's distance. """ # Weight matrix W = _weight_matrix(model) # Hat matrix H = _hat_matrix(X, W) hii = np.diag(H) # Diagonal values of hat matrix # fit.get_influence().hat_matrix_diag # Pearson residuals r = model.resid_pearson # Cook's distance (formula used by R = (res/(1 - hat))^2 * hat/(dispersion * p)) # Note: dispersion is 1 since we aren't modeling overdispersion resid = good.disp - model.predict(good) rss = np.sum(resid**2) MSE = rss / (good.shape[0] - 2) # use the formula from: https://www.mathworks.com/help/stats/cooks-distance.html cooks_d = r**2 / (2 * MSE) * hii / (1 - hii) ** 2 # (r / (1 - hii)) ** 2 * / (1 * 2) return cooks_d # --------------------------------------------------------------------------------------------------- # preprocess utilities
[docs]def filter_genes_by_pattern( adata: anndata.AnnData, patterns: Tuple[str] = ("MT-", "RPS", "RPL", "MRPS", "MRPL", "ERCC-"), drop_genes: bool = False, ) -> Union[List[bool], None]: """Utility function to filter mitochondria, ribsome protein and ERCC spike-in genes, etc. Args: adata: an AnnData object. patterns: the patterns used to filter genes. Defaults to ("MT-", "RPS", "RPL", "MRPS", "MRPL", "ERCC-"). drop_genes: whether inplace drop the genes from the AnnData object. Defaults to False. Returns: A list of indices of matched genes if `drop_genes` is False. Otherwise, returns none. """ logger = LoggerManager.gen_logger("dynamo-utils") matched_genes = pd.Series(adata.var_names).str.startswith(patterns).to_list() logger.info( "total matched genes is " + str(sum(matched_genes)), indent_level=1, ) if sum(matched_genes) > 0: if drop_genes: gene_bools = np.ones(adata.n_vars, dtype=bool) gene_bools[matched_genes] = False logger.info( "inplace subset matched genes ... ", indent_level=1, ) # let us ignore the `inplace` parameter in pandas.Categorical.remove_unused_categories warning. with warnings.catch_warnings(): warnings.simplefilter("ignore") adata._inplace_subset_var(gene_bools) logger.finish_progress(progress_name="filter_genes_by_pattern") return None else: logger.finish_progress(progress_name="filter_genes_by_pattern") return matched_genes
def basic_stats(adata: anndata.AnnData) -> None: """Generate basic stats of the adata, including number of genes, number of cells, and number of mitochondria genes. Args: adata: an AnnData object. """ adata.obs["nGenes"], adata.obs["nCounts"] = np.array((adata.X > 0).sum(1)), np.array((adata.X).sum(1)) adata.var["nCells"], adata.var["nCounts"] = np.array((adata.X > 0).sum(0).T), np.array((adata.X).sum(0).T) if adata.var_names.inferred_type == "bytes": adata.var_names = adata.var_names.astype("str") mito_genes = adata.var_names.str.upper().str.startswith("MT-") if sum(mito_genes) > 0: try: adata.obs["pMito"] = np.array(adata.X[:, mito_genes].sum(1) / adata.obs["nCounts"].values.reshape((-1, 1))) except: # noqa E722 main_exception( "no mitochondria genes detected; looks like your var_names may be corrupted (i.e. " "include nan values). If you don't believe so, please report to us on github or " "via xqiu@wi.mit.edu" ) else: adata.obs["pMito"] = 0 def unique_var_obs_adata(adata: anndata.AnnData) -> anndata.AnnData: """Function to make the obs and var attribute's index unique Args: adata: an AnnData object. Returns: The updated annData object. """ adata.obs_names_make_unique() adata.var_names_make_unique() return adata def convert_layers2csr(adata: anndata.AnnData) -> anndata.AnnData: """Function to convert a layer of sparse matrix to compressed csr_matrix. Args: adata: an AnnData object. Returns: The updated annData object. """ for key in adata.layers.keys(): adata.layers[key] = csr_matrix(adata.layers[key]) if not issparse(adata.layers[key]) else adata.layers[key] return adata def merge_adata_attrs(adata_ori: AnnData, adata: AnnData, attr: Literal["var", "obs"]) -> AnnData: """Merge two adata objects. Args: adata_ori: an AnnData object to be merged into. adata: the AnnData object to be merged. attr: the attribution of adata to be merged, either "var" or "obs". Returns: The merged AnnData object. """ def _merge_by_diff(origin_df: pd.DataFrame, diff_df: pd.DataFrame) -> pd.DataFrame: """Merge two DatafFames. Args: origin_df: the DataFrame to be merged into. diff_df: the DataFrame to be merged. Returns: The merged DataFrame. """ _columns = set(diff_df.columns).difference(origin_df.columns) new_df = origin_df.merge(diff_df[_columns], how="left", left_index=True, right_index=True) return new_df.loc[origin_df.index, :] if attr == "var": adata_ori.var = _merge_by_diff(adata_ori.var, adata.var) elif attr == "obs": obs_df = _merge_by_diff(adata_ori.obs, adata.obs) if obs_df.shape[0] > adata_ori.n_obs: raise ValueError( "Left join generates more rows. Please check if you obs names are unique before calling this fucntion." ) adata_ori.obs = obs_df return adata_ori def get_inrange_shared_counts_mask( adata: anndata.AnnData, layers: List[str], min_shared_count: int, count_by: Literal["gene", "cells"] = "gene" ) -> np.ndarray: """Generate the mask showing the genes having counts more than the provided minimal count. Args: adata: an AnnData object. layers: the layers to be operated on. min_shared_count: the minimal shared number of counts for each genes across cell between layers. count_by: the count type of the data, either "gene: or "cells". Defaults to "gene". Raises: ValueError: invalid count type. Returns: The result mask showing the genes having counts more than the provided minimal count. """ layers = list(set(layers).difference(["X", "matrix", "ambiguous", "spanning"])) # choose shared counts sum by row or columns based on type: `gene` or `cells` sum_dim_index = None ret_dim_index = None if count_by == "gene": sum_dim_index = 0 ret_dim_index = 1 elif count_by == "cells": sum_dim_index = 1 ret_dim_index = 0 else: raise ValueError("Not supported shared account type") if len(np.array(layers)) == 0: main_warning("No layers exist in adata, skipp filtering by shared counts") return np.repeat(True, adata.shape[ret_dim_index]) layers = np.array(layers)[~pd.DataFrame(layers)[0].str.startswith("X_").values] _nonzeros, _sum = None, None # TODO fix bug: when some layers are sparse and some others are not (mixed sparse and ndarray), if the first one happens to be sparse, # dimension mismatch error will be raised; if the first layer (layers[0]) is not sparse, then the following loop works fine. # also check if layers2csr() function works for layer in layers: main_debug(adata.layers[layer].shape) main_debug("layer: %s" % layer) if issparse(adata.layers[layers[0]]): main_debug("when sparse, layer type:" + str(type(adata.layers[layer]))) _nonzeros = adata.layers[layer] > 0 if _nonzeros is None else _nonzeros.multiply(adata.layers[layer] > 0) else: main_debug("when not sparse, layer type:" + str(type(adata.layers[layer]))) _nonzeros = adata.layers[layer] > 0 if _nonzeros is None else _nonzeros * (adata.layers[layer] > 0) for layer in layers: if issparse(adata.layers[layers[0]]): _sum = ( _nonzeros.multiply(adata.layers[layer]) if _sum is None else _sum + _nonzeros.multiply(adata.layers[layer]) ) else: _sum = ( np.multiply(_nonzeros, adata.layers[layer]) if _sum is None else _sum + np.multiply(_nonzeros, adata.layers[layer]) ) return ( np.array(_sum.sum(sum_dim_index).A1 >= min_shared_count) if issparse(adata.layers[layers[0]]) else np.array(_sum.sum(sum_dim_index) >= min_shared_count) ) def clusters_stats( U: pd.DataFrame, S: pd.DataFrame, clusters_uid: np.ndarray, cluster_ix: np.ndarray, size_limit: int = 40 ) -> Tuple[np.ndarray, np.ndarray]: """Calculate the averages per cluster for unspliced and spliced data. Args: U: the unspliced DataFrame. S: the spliced DataFrame. clusters_uid: the uid of the clusters. cluster_ix: the indices of the clusters in adata.obs. size_limit: the max number of members to be considered in a cluster during calculation. Defaults to 40. Returns: U_avgs: the average of clusters for unspliced data. S_avgs: the average of clusters for spliced data. """ U_avgs = np.zeros((S.shape[1], len(clusters_uid))) S_avgs = np.zeros((S.shape[1], len(clusters_uid))) # avgU_div_avgS = np.zeros((S.shape[1], len(clusters_uid))) # slopes_by_clust = np.zeros((S.shape[1], len(clusters_uid))) for i, uid in enumerate(clusters_uid): cluster_filter = cluster_ix == i n_cells = np.sum(cluster_filter) if n_cells > size_limit: U_avgs[:, i], S_avgs[:, i] = ( U[cluster_filter, :].mean(0), S[cluster_filter, :].mean(0), ) else: U_avgs[:, i], S_avgs[:, i] = U.mean(0), S.mean(0) return U_avgs, S_avgs def get_svr_filter( adata: anndata.AnnData, layer: str = "spliced", n_top_genes: int = 3000, return_adata: bool = False ) -> Union[anndata.AnnData, np.ndarray]: """Generate the mask showing the genes with valid svr scores. Args: adata: an AnnData object. layer: the layer to operate on. Defaults to "spliced". n_top_genes: number of top genes to be filtered. Defaults to 3000. return_adata: whether return an updated AnnData or the mask as an array. Defaults to False. Returns: adata: updated adata object with the mask. filter_bool: the filter mask as a bool array. """ score_name = "score" if layer in ["X", "all"] else layer + "_score" valid_idx = np.where(np.isfinite(adata.var.loc[:, score_name]))[0] valid_table = adata.var.iloc[valid_idx, :] nth_score = np.sort(valid_table.loc[:, score_name])[::-1][np.min((n_top_genes - 1, valid_table.shape[0] - 1))] feature_gene_idx = np.where(valid_table.loc[:, score_name] >= nth_score)[0][:n_top_genes] feature_gene_idx = valid_idx[feature_gene_idx] if return_adata: adata.var.loc[:, "use_for_pca"] = False adata.var.loc[adata.var.index[feature_gene_idx], "use_for_pca"] = True res = adata else: filter_bool = np.zeros(adata.n_vars, dtype=bool) filter_bool[feature_gene_idx] = True res = filter_bool return res def sz_util( adata: anndata.AnnData, layer: str, round_exprs: bool, method: Literal["mean-geometric-mean-total", "geometric", "median"], locfunc: Callable, total_layers: List[str] = None, CM: pd.DataFrame = None, scale_to: Union[float, None] = None, ) -> Tuple[pd.Series, pd.Series]: """Calculate the size factor for a given layer. Args: adata: an AnnData object. layer: the layer to operate on. round_exprs: whether the gene expression should be rounded into integers. method: the method used to calculate the expected total reads / UMI used in size factor calculation. Only `mean-geometric-mean-total` / `geometric` and `median` are supported. When `median` is used, `locfunc` will be replaced with `np.nanmedian`. locfunc: the function to normalize the data. total_layers: the layer(s) that can be summed up to get the total mRNA. For example, ["spliced", "unspliced"], ["uu", "ul", "su", "sl"] or ["new", "old"], etc. Defaults to None. CM: the data to operate on, overriding the layer. Defaults to None. scale_to: the final total expression for each cell that will be scaled to. Defaults to None. Raises: NotImplementedError: method is invalid. Returns: A tuple (sfs, cell_total) where sfs is the size factors and cell_total is the initial cell size. """ adata = adata.copy() if layer == "_total_" and "_total_" not in adata.layers.keys(): if total_layers is not None: if not isinstance(total_layers, list): total_layers = [total_layers] if len(set(total_layers).difference(adata.layers.keys())) == 0: total = None for t_key in total_layers: total = adata.layers[t_key] if total is None else total + adata.layers[t_key] adata.layers["_total_"] = total if layer == "raw": CM = adata.raw.X if CM is None else CM elif layer == "X": CM = adata.X if CM is None else CM elif layer == "protein": if "protein" in adata.obsm_keys(): CM = adata.obsm["protein"] if CM is None else CM else: return None, None else: CM = adata.layers[layer] if CM is None else CM if round_exprs: main_info("rounding expression data of layer: %s during size factor calculation" % (layer)) if issparse(CM): CM.data = np.round(CM.data, 0) else: CM = CM.round().astype("int") cell_total = CM.sum(axis=1).A1 if issparse(CM) else CM.sum(axis=1) cell_total += cell_total == 0 # avoid infinity value after log (0) if method in ["mean-geometric-mean-total", "geometric"]: sfs = cell_total / (np.exp(locfunc(np.log(cell_total))) if scale_to is None else scale_to) elif method == "median": sfs = cell_total / (np.nanmedian(cell_total) if scale_to is None else scale_to) elif method == "mean": sfs = cell_total / (np.nanmean(cell_total) if scale_to is None else scale_to) else: raise NotImplementedError(f"This method {method} is not supported!") return sfs, cell_total def get_sz_exprs( adata: anndata.AnnData, layer: str, total_szfactor: Union[str, None] = None ) -> Tuple[np.ndarray, npt.ArrayLike]: """Get the size factor from an AnnData object. Args: adata: an AnnData object. layer: the layer for which to get the size factor. total_szfactor: the key-name for total size factor entry in `adata.obs`. If not None, would override the layer selected. Defaults to None. Returns: A tuple (szfactors, CM), where szfactors is the queried size factor and CM is the data of the layer corresponding to the size factor. """ if layer == "raw": CM = adata.raw.X szfactors = adata.obs[layer + "Size_Factor"].values[:, None] elif layer == "X": CM = adata.X szfactors = adata.obs["Size_Factor"].values[:, None] elif layer == "protein": if "protein" in adata.obsm_keys(): CM = adata.obsm[layer] szfactors = adata.obs["protein_Size_Factor"].values[:, None] else: CM, szfactors = None, None else: CM = adata.layers[layer] szfactors = adata.obs[layer + "_Size_Factor"].values[:, None] if total_szfactor is not None and total_szfactor in adata.obs.keys(): szfactors = adata.obs[total_szfactor][:, None] elif total_szfactor is not None: main_warning("`total_szfactor` is not `None` and it is not in adata object.") return szfactors, CM def normalize_mat_monocle( mat: np.ndarray, szfactors: np.ndarray, relative_expr: bool, pseudo_expr: int, norm_method: Callable = np.log1p ) -> np.ndarray: """Normalize the given array for monocle recipe. Args: mat: the array to operate on. szfactors: the size factors corresponding to the array. relative_expr: whether we need to divide gene expression values first by size factor before normalization. pseudo_expr: a pseudocount added to the gene expression value before log/log2 normalization. norm_method: the method used to normalize data. Defaults to np.log1p. Returns: The normalized array. """ if norm_method == np.log1p: pseudo_expr = 0 if relative_expr: mat = mat.multiply(csr_matrix(1 / szfactors)) if issparse(mat) else mat / szfactors if pseudo_expr is None: pseudo_expr = 1 if issparse(mat): mat.data = norm_method(mat.data + pseudo_expr) if norm_method is not None else mat.data if norm_method is not None and norm_method.__name__ == "Freeman_Tukey": mat.data -= 1 else: mat = norm_method(mat + pseudo_expr) if norm_method is not None else mat return mat def Freeman_Tukey(X: np.ndarray, inverse=False) -> np.ndarray: """perform Freeman-Tukey transform or inverse transform on the given array. Args: X: a matrix. inverse: whether to perform inverse Freeman-Tukey transform. Defaults to False. Returns: The transformed array. """ if inverse: res = np.sqrt(X) + np.sqrt((X + 1)) else: res = (X**2 - 1) ** 2 / (4 * X**2) return res def anndata_bytestring_decode(adata_item: pd.DataFrame) -> None: """Decode contents of an annotation of an AnnData object inplace. Args: adata_item: an annotation of an AnnData object. """ for key in adata_item.keys(): df = adata_item[key] if df.dtype.name == "category" and areinstance(df.cat.categories, bytes): cat = [c.decode() for c in df.cat.categories] df.cat.rename_categories(cat, inplace=True) def decode_index(adata_item: pd.DataFrame) -> None: """Decode indices of an annotation of an AnnData object inplace. Args: adata_item: an annotation of an AnnData object. """ if areinstance(adata_item.index, bytes): index = {i: i.decode() for i in adata_item.index} adata_item.rename(index, inplace=True) def decode(adata: anndata.AnnData) -> None: """Decode an AnnData object. Args: adata: an AnnData object. """ decode_index(adata.obs) decode_index(adata.var) anndata_bytestring_decode(adata.obs) anndata_bytestring_decode(adata.var) # --------------------------------------------------------------------------------------------------- # pca
[docs]def pca_monocle( adata: AnnData, X_data: np.ndarray = None, n_pca_components: int = 30, pca_key: str = "X", pcs_key: str = "PCs", genes_to_append: Union[List[str], None] = None, layer: Union[List[str], str, None] = None, return_all: bool = False, ) -> Union[AnnData, Tuple[AnnData, Union[PCA, TruncatedSVD], np.ndarray]]: """Perform PCA reduction for monocle recipe. Args: adata: an AnnData object. X_data: the data to perform dimension reduction on. Defaults to None. n_pca_components: number of PCA components reduced to. Defaults to 30. pca_key: the key to store the reduced data. Defaults to "X". pcs_key: the key to store the principle axes in feature space. Defaults to "PCs". genes_to_append: a list of genes should be inspected. Defaults to None. layer: the layer(s) to perform dimension reduction on. Would be overrided by X_data. Defaults to None. return_all: whether to return the PCA fit model and the reduced array together with the updated AnnData object. Defaults to False. Raises: ValueError: layer provided is not invalid. ValueError: list of genes to append is invalid. Returns: The the updated AnnData object with reduced data if `return_all` is False. Otherwise, a tuple (adata, fit, X_pca), where adata is the updated AnnData object, fit is the fit model for dimension reduction, and X_pca is the reduced array, will be returned. """ # only use genes pass filter (based on use_for_pca) to perform dimension reduction. if X_data is None: if "use_for_pca" not in adata.var.keys(): adata.var["use_for_pca"] = True if layer is None: X_data = adata.X[:, adata.var.use_for_pca.values] else: if "X" in layer: X_data = adata.X[:, adata.var.use_for_pca.values] elif "total" in layer: X_data = adata.layers["X_total"][:, adata.var.use_for_pca.values] elif "spliced" in layer: X_data = adata.layers["X_spliced"][:, adata.var.use_for_pca.values] elif "protein" in layer: X_data = adata.obsm["X_protein"] elif type(layer) is str: X_data = adata.layers["X_" + layer][:, adata.var.use_for_pca.values] else: raise ValueError( f"your input layer argument should be either a `str` or a list that includes one of `X`, " f"`total`, `protein` element. `Layer` currently is {layer}." ) cm_genesums = X_data.sum(axis=0) valid_ind = np.logical_and(np.isfinite(cm_genesums), cm_genesums != 0) valid_ind = np.array(valid_ind).flatten() bad_genes = np.where(adata.var.use_for_pca)[0][~valid_ind] if genes_to_append is not None and len(adata.var.index[bad_genes].intersection(genes_to_append)) > 0: raise ValueError( f"The gene list passed to argument genes_to_append contains genes with no expression " f"across cells or non finite values. Please check those genes:" f"{set(bad_genes).intersection(genes_to_append)}!" ) adata.var.iloc[bad_genes, adata.var.columns.tolist().index("use_for_pca")] = False X_data = X_data[:, valid_ind] USE_TRUNCATED_SVD_THRESHOLD = 100000 if adata.n_obs < USE_TRUNCATED_SVD_THRESHOLD: pca = PCA( n_components=min(n_pca_components, X_data.shape[1] - 1), svd_solver="arpack", random_state=0, ) fit = pca.fit(X_data.toarray()) if issparse(X_data) else pca.fit(X_data) X_pca = fit.transform(X_data.toarray()) if issparse(X_data) else fit.transform(X_data) adata.obsm[pca_key] = X_pca adata.uns[pcs_key] = fit.components_.T adata.uns["explained_variance_ratio_"] = fit.explained_variance_ratio_ else: # unscaled PCA fit = TruncatedSVD( n_components=min(n_pca_components + 1, X_data.shape[1] - 1), random_state=0, ) # first columns is related to the total UMI (or library size) X_pca = fit.fit_transform(X_data)[:, 1:] adata.obsm[pca_key] = X_pca adata.uns[pcs_key] = fit.components_.T adata.uns["explained_variance_ratio_"] = fit.explained_variance_ratio_[1:] adata.uns["explained_variance_ratio_"] = fit.explained_variance_ratio_[1:] adata.uns["pca_mean"] = fit.mean_ if hasattr(fit, "mean_") else None if return_all: return adata, fit, X_pca else: return adata
def pca_genes(PCs: list, n_top_genes: int = 100) -> np.ndarray: """For each gene, if the gene is n_top in some principle component then it is valid. Return all such valid genes. Args: PCs: principle components(PC) of PCA n_top_genes: number of gene candidates in EACH PC. Defaults to 100. Returns: A bool array indicating whether the gene is valid. """ valid_genes = np.zeros(PCs.shape[0], dtype=bool) for q in PCs.T: sorted_q = np.sort(np.abs(q))[::-1] is_pc_top_n = np.abs(q) > sorted_q[n_top_genes] valid_genes = np.logical_or(is_pc_top_n, valid_genes) return valid_genes
[docs]def top_pca_genes( adata: AnnData, pc_key: str = "PCs", n_top_genes: int = 100, pc_components: Union[int, None] = None, adata_store_key: str = "top_pca_genes", ) -> AnnData: """Define top genes as any gene that is ``n_top_genes`` in some principle component. Args: adata: an AnnData object. pc_key: component key stored in adata.uns. Defaults to "PCs". n_top_genes: number of top genes as valid top genes in each component. Defaults to 100. pc_components: number of top principle components to use. Defaults to None. adata_store_key: the key for storing pca genes. Defaults to "top_pca_genes". Raises: Exception: invalid pc_key. Returns: The AnnData object with top genes stored as values of adata.var[adata_store_key]. """ if pc_key in adata.uns.keys(): Q = adata.uns[pc_key] elif pc_key in adata.varm.keys(): Q = adata.varm[pc_key] else: raise Exception(f"No PC matrix {pc_key} found in neither .uns nor .varm.") if pc_components is not None: if type(pc_components) == int: Q = Q[:, :pc_components] elif type(pc_components) == list: Q = Q[:, pc_components] pcg = pca_genes(Q, n_top_genes=n_top_genes) genes = np.zeros(adata.n_vars, dtype=bool) if DKM.VAR_USE_FOR_PCA in adata.var.keys(): genes[adata.var[DKM.VAR_USE_FOR_PCA]] = pcg else: genes = pcg main_info_insert_adata_var(adata_store_key, indent_level=2) adata.var[adata_store_key] = genes return adata
def add_noise_to_duplicates(adata: anndata.AnnData, basis: str = "pca") -> None: """Add noise to duplicated elements of the reduced array inplace. Args: adata: an AnnData object. basis: the type of dimension redduction. Defaults to "pca". """ X_data = adata.obsm["X_" + basis] min_val = abs(X_data).min() n_obs, n_var = X_data.shape while True: _, index = np.unique(X_data, axis=0, return_index=True) duplicated_idx = np.setdiff1d(np.arange(n_obs), index) if len(duplicated_idx) == 0: adata.obsm["X_" + basis] = X_data break else: X_data[duplicated_idx, :] += np.random.normal(0, min_val / 1000, (len(duplicated_idx), n_var)) # --------------------------------------------------------------------------------------------------- # labeling related def collapse_species_adata(adata: anndata.AnnData) -> None: """Function to collapse the four species data, will be generalized to handle dual-datasets. Args: adata: an AnnData object. """ ( only_splicing, only_labeling, splicing_and_labeling, ) = DKM.allowed_layer_raw_names() if np.all([name in adata.layers.keys() for name in splicing_and_labeling]): if only_splicing[0] not in adata.layers.keys(): adata.layers[only_splicing[0]] = adata.layers["su"] + adata.layers["sl"] if only_splicing[1] not in adata.layers.keys(): adata.layers[only_splicing[1]] = adata.layers["uu"] + adata.layers["ul"] if only_labeling[0] not in adata.layers.keys(): adata.layers[only_labeling[0]] = adata.layers["ul"] + adata.layers["sl"] if only_labeling[1] not in adata.layers.keys(): adata.layers[only_labeling[1]] = adata.layers[only_labeling[0]] + adata.layers["uu"] + adata.layers["su"] return adata def detect_experiment_datatype(adata: anndata.AnnData) -> Tuple[bool, bool, bool, bool]: """Tells what kinds of experiment data are stored in an AnnData object. Args: adata: an AnnData object. Returns: A tuple (has_splicing, has_labeling, splicing_labeling, has_protein), where has_splicing represents whether the object containing unspliced and spliced data, has_labeling represents whether the object containing new expression and total expression (i.e. labelling) data, splicing_labeling represents whether the object containing both splicing and labelling data, and has_protein represents whether the object containing protein data. """ has_splicing, has_labeling, splicing_labeling, has_protein = ( False, False, False, False, ) layers = adata.layers.keys() if ( len({"ul", "sl", "uu", "su"}.difference(layers)) == 0 or len({"X_ul", "X_sl", "X_uu", "X_su"}.difference(layers)) == 0 ): has_splicing, has_labeling, splicing_labeling = True, True, True elif ( len({"unspliced", "spliced", "new", "total"}.difference(layers)) == 0 or len({"X_unspliced", "X_spliced", "X_new", "X_total"}.difference(layers)) == 0 ): has_splicing, has_labeling = True, True elif ( len({"unspliced", "spliced"}.difference(layers)) == 0 or len({"X_unspliced", "X_spliced"}.difference(layers)) == 0 ): has_splicing = True elif len({"new", "total"}.difference(layers)) == 0 or len({"X_new", "X_total"}.difference(layers)) == 0: has_labeling = True if "protein" in adata.obsm.keys(): has_protein = True return has_splicing, has_labeling, splicing_labeling, has_protein def default_layer(adata: anndata.AnnData) -> str: """Returns the defualt layer preferred in a given AnnData object. Args: adata: an AnnData object. Returns: The key of the default layer. """ has_splicing, has_labeling, splicing_labeling, _ = detect_experiment_datatype(adata) if has_splicing: if has_labeling: if len(set(adata.layers.keys()).intersection(["new", "total", "spliced", "unspliced"])) == 4: adata = collapse_species_adata(adata) default_layer = ( "M_t" if "M_t" in adata.layers.keys() else "X_total" if "X_total" in adata.layers.keys() else "total" ) else: default_layer = ( "M_s" if "M_s" in adata.layers.keys() else "X_spliced" if "X_spliced" in adata.layers.keys() else "spliced" ) else: default_layer = ( "M_t" if "M_t" in adata.layers.keys() else "X_total" if "X_total" in adata.layers.keys() else "total" ) return default_layer def calc_new_to_total_ratio(adata: anndata.AnnData) -> Union[Tuple[np.ndarray, np.ndarray], Tuple[None, None]]: """Calculate the new to total ratio across cells. Note that NTR for the first time point in degradation approximates gamma/beta. Args: adata: an AnnData object. Returns: ntr: the new to total ratio of all genes for each cell. Returned if the object has labelling or splicing layers. var_ntr: the new to total ratio of all cells for each gene. Returned if the object has labelling or splicing layers. """ if len({"new", "total"}.intersection(adata.layers.keys())) == 2: ntr = adata.layers["new"].sum(1) / adata.layers["total"].sum(1) ntr = ntr.A1 if issparse(adata.layers["new"]) else ntr var_ntr = adata.layers["new"].sum(0) / adata.layers["total"].sum(0) var_ntr = var_ntr.A1 if issparse(adata.layers["new"]) else var_ntr elif len({"uu", "ul", "su", "sl"}.intersection(adata.layers.keys())) == 4: new = adata.layers["ul"].sum(1) + adata.layers["sl"].sum(1) total = new + adata.layers["uu"].sum(1) + adata.layers["su"].sum(1) ntr = new / total ntr = ntr.A1 if issparse(adata.layers["uu"]) else ntr new = adata.layers["ul"].sum(0) + adata.layers["sl"].sum(0) total = new + adata.layers["uu"].sum(0) + adata.layers["su"].sum(0) var_ntr = new / total var_ntr = var_ntr.A1 if issparse(adata.layers["uu"]) else var_ntr elif len({"unspliced", "spliced"}.intersection(adata.layers.keys())) == 2: with warnings.catch_warnings(): warnings.simplefilter("ignore") ntr = adata.layers["unspliced"].sum(1) / (adata.layers["unspliced"] + adata.layers["spliced"]).sum(1) var_ntr = adata.layers["unspliced"].sum(0) / (adata.layers["unspliced"] + adata.layers["spliced"]).sum(0) ntr = ntr.A1 if issparse(adata.layers["unspliced"]) else ntr var_ntr = var_ntr.A1 if issparse(adata.layers["unspliced"]) else var_ntr else: ntr, var_ntr = None, None return ntr, var_ntr
[docs]def scale( adata: anndata.AnnData, layers: Union[List[str], str, None] = None, scale_to_layer: Union[str, None] = None, scale_to: float = 1e6, ) -> anndata.AnnData: """Scale layers to a particular total expression value, similar to `normalize_expr_data` function. Args: adata: an AnnData object. layers: the layers to scale. Defaults to None. scale_to_layer: use which layer to calculate a global scale factor. If None, calculate each layer's own scale factor and scale all layers to same total value. Defaults to None. scale_to: the total expression value that layers are scaled to. Defaults to 1e6. Returns: The scaled AnnData object. """ layers = DynamoAdataKeyManager.get_available_layer_keys(adata, layers) has_splicing, has_labeling, _ = detect_experiment_datatype(adata) if scale_to_layer is None: scale_to_layer = "total" if has_labeling else None scale = scale_to / adata.layers[scale_to_layer].sum(1) else: scale = None for layer in layers: if scale is None: scale = scale_to / adata.layers[layer].sum(1) adata.layers[layer] = csr_matrix(adata.layers[layer].multiply(scale)) return adata
# --------------------------------------------------------------------------------------------------- # ERCC related def relative2abs( adata: anndata.AnnData, dilution: float, volume: float, from_layer: Union[str, None] = None, to_layers: Union[str, List[str], None] = None, mixture_type: Literal[1, 2] = 1, ERCC_controls: Union[np.ndarray, None] = None, ERCC_annotation: Union[pd.DataFrame, None] = None, ) -> anndata.AnnData: """Converts FPKM/TPM data to transcript counts using ERCC spike-in. This is based on the relative2abs function from monocle 2 (Qiu, et. al, Nature Methods, 2017). Args: adata: an Annodata object. dilution: the dilution of the spikein transcript in the lysis reaction mix. Default is 40, 000. The number of spike-in transcripts per single-cell lysis reaction was calculated from. volume: the approximate volume of the lysis chamber (nanoliters). from_layer: the layer in which the ERCC TPM values will be used as the covariate for the ERCC based linear regression. Defaults to None. to_layers: the layers that our ERCC based transformation will be applied to. Defaults to None. mixture_type: the type of spikein transcripts from the spikein mixture added in the experiments. Note that m/c we inferred are also based on mixture 1. Defaults to 1. ERCC_controls: the FPKM/TPM matrix for each ERCC spike-in transcript in the cells if user wants to perform the transformation based on their spike-in data. Note that the row and column names should match up with the ERCC_annotation and relative_exprs_matrix respectively. Defaults to None. ERCC_annotation: the ERCC_annotation matrix from illumina USE GUIDE which will be ued for calculating the ERCC transcript copy number for performing the transformation. Defaults to None. Raises: Exception: the number of ERCC gene in `ERCC_annotation["ERCC ID"]` is not enough. Exception: the layers specified in to_layers are invalid. Returns: An adata object with the data specified in the to_layers transformed into absolute counts. """ if ERCC_annotation is None: ERCC_annotation = pd.read_csv( "https://www.dropbox.com/s/cmiuthdw5tt76o5/ERCC_specification.txt?dl=1", sep="\t", ) ERCC_id = ERCC_annotation["ERCC ID"] ERCC_id = adata.var_names.intersection(ERCC_id) if len(ERCC_id) < 10 and ERCC_controls is None: raise Exception("The adata object you provided has less than 10 ERCC genes.") if to_layers is not None: to_layers = [to_layers] if to_layers is str else to_layers to_layers = list(set(adata.layers.keys()).intersection(to_layers)) if len(to_layers) == 0: raise Exception( f"The layers {to_layers} that will be converted to absolute counts doesn't match any layers" f"from the adata object." ) mixture_name = ( "concentration in Mix 1 (attomoles/ul)" if mixture_type == 1 else "concentration in Mix 2 (attomoles/ul)" ) ERCC_annotation["numMolecules"] = ERCC_annotation.loc[:, mixture_name] * ( volume * 10 ** (-3) * 1 / dilution * 10 ** (-18) * 6.02214129 * 10 ** (23) ) ERCC_annotation["rounded_numMolecules"] = ERCC_annotation["numMolecules"].astype(int) if from_layer in [None, "X"]: X, X_ercc = ( adata.X, adata[:, ERCC_id].X if ERCC_controls is None else ERCC_controls, ) else: X, X_ercc = ( adata.layers[from_layer], adata[:, ERCC_id] if ERCC_controls is None else ERCC_controls, ) logged = False if X.max() > 100 else True if not logged: X, X_ercc = ( np.log1p(X.A) if issparse(X_ercc) else np.log1p(X), np.log1p(X_ercc.A) if issparse(X_ercc) else np.log1p(X_ercc), ) else: X, X_ercc = ( X.A if issparse(X_ercc) else X, X_ercc.A if issparse(X_ercc) else X_ercc, ) y = np.log1p(ERCC_annotation["numMolecules"]) for i in range(adata.n_obs): X_i, X_ercc_i = X[i, :], X_ercc[i, :] X_i, X_ercc_i = sm.add_constant(X_i), sm.add_constant(X_ercc_i) res = sm.RLM(y, X_ercc_i).fit() k, b = res.params[::-1] if to_layers is None: X = adata.X logged = False if X.max() > 100 else True if not logged: X_i = np.log1p(X[i, :].A) if issparse(X) else np.log1p(X[i, :]) else: X_i = X[i, :].A if issparse(X) else X[i, :] res = k * X_i + b res = res if logged else np.expm1(res) adata.X[i, :] = csr_matrix(res) if issparse(X) else res else: for cur_layer in to_layers: X = adata.layers[cur_layer] logged = False if X.max() > 100 else True if not logged: X_i = np.log1p(X[i, :].A) if issparse(X) else np.log1p(X[i, :]) else: X_i = X[i, :].A if issparse(X) else X[i, :] res = k * X_i + b if logged else np.expm1(k * X_i + b) adata.layers[cur_layer][i, :] = csr_matrix(res) if issparse(X) else res # --------------------------------------------------------------------------------------------------- # coordinate/vector space operations def affine_transform(X: npt.ArrayLike, A: npt.ArrayLike, b: npt.ArrayLike) -> np.ndarray: """Perform affine trasform on an array. Args: X: the array to tranform. A: the scaling/rotation/shear matrix. b: the transformation matrix. Returns: The result array. """ X = np.array(X) A = np.array(A) b = np.array(b) return (A @ X.T).T + b def gen_rotation_2d(degree: float) -> np.ndarray: """Calculate the 2D rotation transform matrix for given rotation in degrees. Args: degree: the degrees to rotate. Returns: The rotation matrix. """ from math import cos, radians, sin rad = radians(degree) R = [ [cos(rad), -sin(rad)], [sin(rad), cos(rad)], ] return np.array(R)