函数
import scanpy as sc
sc.pp.highly_variable_genes
功能
取出高可变基因,默认使用log的数据,当使用flavor=seurat_v3的时候,采用count data。(这里一定要注意,如果你先对数据做了标准化,再选择seurat_v3将会报错)
flavor参数可以选择是使用Seurat,Cell ranger还是seurat v3的算法。
Seurat and Cellranger中,使用的是dispersion-based方法,获得归一化的方差。先对基因按照表达量平均值进行分bin,然后计算落在每个bin的基因的离散度(dispersion)的均值和SD,最终获得归一化的dispersion。对于每个表达量的bin,选择不同的高可变表达基因。
而Seurat3的算法,计算每个基因的方差进行归一化。首先对数据在规范化标准偏差下(a regularized standard deviation)进行标准化(使用z标准化),之后计算每个基因的归一化的方差,并且进行排序,获得高可变基因。
# 重要参数 """ * adata:输入的数据,每行是一个细胞,每列是一个特征 * layer:使用的是哪一个layer * n_top_genes:如果是使用seurate_v3的方法,那么需要指定该参数。 * min_mean:默认0.0125 ;max_mean:默认是3 ;min_disp: 默认0.5, max_disp: 默认是inf。如果指定了n_top_genes , 这个和其他所有mean和disp参数都会无效,因此设置了 flavor='seurat_v3' 该参数无用。 * span:默认是0.3;当flavor=seurat_v3的时候,用loess模型来估计variance的数据的比例。 * n_bins : 默认是20,对表达量分bin的数目,对每个bin里的数据进行归一化,如果只有一个基因落到bin里,那么该bin的dispersion会设为1。 * flavor: {‘seurat’, ‘cell_ranger’, ‘seurat_v3’} (default: 'seurat') * subset:默认是false,只是返回高可变基因,否则就原位替换 * inplace:默认是True,在var中进行存储矩阵 * batch_key: If specified, highly-variable genes are selected within each batch separately and merged. This simple process avoids the selection of batch-specific genes and acts as a lightweight batch correction method. For all flavors, genes are first sorted by how many batches they are a HVG. For dispersion-based flavors ties are broken by normalized dispersion. If flavor = 'seurat_v3', ties are broken by the median (across batches) rank based on within-batch normalized variance * check_values:True,在seurat_v3模式下有用,检测每个count是不是为整型 """
代码
## _highly_variable_genes.py mean, var = materialize_as_ndarray(_get_mean_var(X)) # now actually compute the dispersion mean[mean == 0] = 1e-12 # set entries equal to zero to small value dispersion = var / mean df['dispersions_norm'] = ( df['dispersions'].values # use values here as index differs - disp_mean_bin[df['mean_bin'].values].values ) / disp_std_bin[df['mean_bin'].values].values
获得每个基因的dispersion值,并进行排序
mean, var = _get_mean_var(X_batch) not_const = var > 0 estimat_var = np.zeros(X.shape[1], dtype=np.float64) y = np.log10(var[not_const]) x = np.log10(mean[not_const]) model = loess(x, y, span=span, degree=2) ### 对mean和var进行loess回归 model.fit() estimat_var[not_const] = model.outputs.fitted_values reg_std = np.sqrt(10 ** estimat_var) batch_counts = X_batch.astype(np.float64).copy()