Select a result to preview

生信分析 - 单细胞 bbknn 去批次（python 和 r）

在 python 中用 bbknn

（一）库和环境

安装 bbknn 库：pip install bbknn

import numpy as np
import pandas as pd
import anndata as ad #导入 anndata 库，并将其命名为 ad，类似 Seurat

import scanpy as sc
import seaborn as sns
import matplotlib.pyplot as pt
import bbknn

# A. 环境设置
sc.settings.verbosity = 1  # verbosity 日志等级: errors (0), warnings (1), info (2), hints (3)
sc.logging.print_versions()
sc.settings.set_figure_params(dpi=80, frameon=False, figsize=(3, 3), facecolor='white')

（二）导入数据集

使用离线数据集构建，数据集下载: scanpy.datasets.pbmc68k_reduced — scanpy

# 若使用在线数据集：
# adata_ref = sc.datasets.pbmc3k_processed()
# adata = sc.datasets.pbmc3k_reduced()
adata_ref = sc.read_10x_mtx('./data/pbmc68k_processed') # ref 批次
adata = sc.read_10x_mtx('./data/pbmc68k_reduced') # new 批次
adata_concat = adata_ref.concatenate(adata, batch_categories=['ref', 'new']) #合并批次
print(adata_concat.obs)

（三）数据预处理

sc.pp.normalize_total(adata_concat, target_sum=1e4) #正则
sc.pp.log1p(adata_concat) #取对数, 在 Seruat 与上一步是合并的
# Use sc.pp.highly_variable_genes
sc.pp.filter_genes(adata_concat,min_cells=1)
sc.pp.filter_genes_dispersion(adata_concat, n_top_genes =2000) #top 2000 gene
sc.pp.scale(adata_concat) #标准化

（四）降维、kkbnn 去批次、聚类和可视化

# 降维
sc.tl.pca(adata_concat, svd_solver="arpack")

# 去批次
import bbknn
sc.external.pp.bbknn(adata_concat, batch_key='batch')

# 聚类和可视化
sc.pp.neighbors(adata_concat, n_neighbors=10, n_pcs=40)
sc.tl.leiden(
    adata_concat,
    resolution=0.9,
    random_state=0,
    flavor="igraph",
    n_iterations=2,
    directed=False,
)
print(adata_concat.obs)
print(adata_concat.obsm)
sc.tl.umap(adata_concat)
sc.pl.umap(adata_concat, color=["leiden"])

在 R 中用 bbknn（基于 reticulate 调用 python 执行）

参考说明：Teichlab/bbknn: Batch balanced KNN

（一）库和环境

原理是通过 reticulate 这个 R 的程序包，在 R 中调用 python 代码执行。因此，需要：

  1. 安装 reticulate R 包
  2. 具备 python 环境及其 bbknn 库（即在 python 中使用 bbknn 是正常的）
  3. 设置 use_python 路径（可在终端通过 which 命令查询）

（二）导入数据集和数据处理

使用 ifnb.data = LoadData("ifnb") 数据集

ifnb.data = LoadData("ifnb")
ifnb = UpdateSeuratObject(ifnb)
ifnb <- NormalizeData(ifnb)
ifnb <- FindVariableFeatures(ifnb)
ifnb <- ScaleData(ifnb)

（三）调用 python 的 pca 降维 和 bbknn 去批次

特别注意！Python Anndata 和 R Seurat 对象存储有区别，因此需要先构建正确的 Anndata 对象，即针对不同的数据结构正确提取 X 数据矩阵和 obs 细胞元信息（包含有批次信息），从而构建 Anndata 对象。

# python 库
anndata = import("anndata",convert=FALSE)
bbknn = import("bbknn", convert=FALSE)
sc = import("scanpy",convert=FALSE)
np <- import("numpy",convert=FALSE)

# 基于 ifnb Seurat 构建 Anndata 对象
X_ <- as.matrix(t(ifnb.data@assays$RNA@layers$data))  # 确保是矩阵格式
obs_ <- as.data.frame(ifnb.data@meta.data)
adata = anndata$AnnData(X=X_, obs=obs_)

# 运行 pca、bbknn 和 umap
sc$tl$pca(adata)
sc$external$pp$bbknn(adata,batch_key="orig.ident")
sc$tl$umap(adata)
umap = py_to_r(adata$obsm[["X_umap"]])

（四）把 umap 结果在 r 中使用

# 提取细胞名称（r 中是列，python 中是行）
cell_names <- colnames(ifnb.data)
# 检查 UMAP 数据的维度
if (nrow(umap) != length(cell_names)) {
  stop("UMAP 数据的行数与 Seurat 对象中的细胞数不匹配。")
}
# 将细胞名称设置为 umap 的行名
rownames(umap) <- cell_names
ifnb.data[["umap"]] <- CreateDimReducObject(embeddings = umap, key = "UMAP_", assay = DefaultAssay(ifnb.data))
p1 <- DimPlot(ifnb.data, reduction = "umap", group.by = "stim")
p2 <- DimPlot(ifnb.data, reduction = "umap", group.by = "seurat_annotations", label = TRUE)
P.total = wrap_plots(p1,p2)
P.total