155_单细胞数据降噪去批¶

import scanpy as sc
import numpy as np

adata = sc.read_h5ad("pbmc3k.h5ad")

# 计算零值比例
total_entries = adata.X.shape[0] * adata.X.shape[1]
zero_entries = (adata.X == 0).sum() if hasattr(adata.X, 'toarray') else np.sum(adata.X.toarray() == 0)
sparsity = zero_entries / total_entries
print(f"数据稀疏度: {sparsity:.1%}")  # 通常 85-95%

# 每个细胞检测到的基因数
genes_per_cell = (adata.X > 0).sum(axis=1)
print(f"平均每个细胞检测到: {np.mean(genes_per_cell):.0f} 基因")
print(f"总基因数: {adata.n_vars}")

# 安装: pip install magic-impute
import magic
import scanpy as sc
import scprep

# 1. 加载数据
adata = sc.read_h5ad("raw_counts.h5ad")

# 2. 预处理 (MAGIC需要归一化后的数据)
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

# 3. 运行MAGIC
magic_op = magic.MAGIC(
    n_pca=20,           # PCA维度(用于KNN图)
    n_jobs=-1,          # 并行
    knn=10,             # K近邻
    t='auto',           # 扩散步数(auto自动选择)
    solver='exact'      # 精确求解(小数据集)
)

# 插补
adata_magic = magic_op.fit_transform(adata.to_df())

# 4. 将结果存回anndata
adata.layers["magic"] = adata_magic.values

# 5. 可视化对比
import matplotlib.pyplot as plt
fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# 原始数据
sc.pl.scatter(adata, x="CD4", y="CD8A", ax=axes[0], title="Original")
# MAGIC后
adata_viz = adata.copy()
adata_viz.X = adata.layers["magic"]
sc.pl.scatter(adata_viz, x="CD4", y="CD8A", ax=axes[1], title="MAGIC")
plt.savefig("magic_comparison.png", dpi=150)

# 6. 调节扩散程度
# t=1: 最少平滑; t=auto: 自适应; t=7: 较强平滑
magic_mild = magic.MAGIC(t=2).fit_transform(adata.to_df())    # 轻微平滑
magic_strong = magic.MAGIC(t=7).fit_transform(adata.to_df())  # 较强平滑

# 安装: pip install dca
from dca.api import dca
import scanpy as sc

# 1. 加载原始counts数据
adata = sc.read_h5ad("raw_counts.h5ad")

# 2. 运行DCA
# DCA修改adata.X为去噪后的数据，同时保存分布参数
adata_dca = dca(
    adata,
    mode='denoise',          # 'denoise'或'latent'
    ae_type='zinb-conddisp', # 模型类型: nb, zinb, zinb-conddisp
    hidden_size=(64, 32, 64),# 隐藏层大小
    hidden_dropout=0.0,      # dropout正则化
    batchnorm=True,          # 批归一化
    epochs=300,              # 训练轮数
    batch_size=32,           # 批次大小
    learning_rate=0.001,     # 学习率
    early_stop=10,           # 早停耐心值
    threads=8,               # CPU线程
    verbose=True
)

# 3. DCA结果
# adata.X 已被替换为去噪后的数据
# adata.obsm['X_dca'] 为瓶颈层的潜在表示
print(f"去噪后零值比例: {(adata.X == 0).sum() / adata.X.size:.1%}")

# 4. 使用DCA潜在表示进行聚类
sc.pp.neighbors(adata, use_rep='X_dca')
sc.tl.umap(adata)
sc.tl.leiden(adata)
sc.pl.umap(adata, color='leiden')

# 安装: install.packages("scImpute") 或 devtools::install_github("Vivianstats/scImpute")
library(scImpute)

# 1. 准备输入 (CSV格式, 基因×细胞)
# scImpute读取文件路径
write.csv(count_matrix, "raw_counts.csv")

# 2. 运行scImpute
scimpute(
  count_path = "raw_counts.csv",       # 输入文件
  infile = "csv",                       # 输入格式
  outfile = "csv",                      # 输出格式
  out_dir = "scimpute_output/",         # 输出目录
  labeled = FALSE,                      # 是否有细胞类型标签
  drop_thre = 0.5,                      # dropout概率阈值(>0.5才插补)
  Kcluster = 5,                         # 预聚类数
  ncores = 8                            # 并行核数
)

# 如果有细胞类型标签，可以利用:
scimpute(
  count_path = "raw_counts.csv",
  infile = "csv",
  outfile = "csv",
  out_dir = "scimpute_output/",
  labeled = TRUE,
  labels = cell_type_labels,    # 细胞类型向量
  drop_thre = 0.5,
  ncores = 8
)

# 3. 读取结果
imputed <- read.csv("scimpute_output/scimpute_count.csv", row.names = 1)

# 4. 对比插补前后
# 插补了多少值
original <- read.csv("raw_counts.csv", row.names = 1)
n_imputed <- sum(original == 0 & imputed != 0)
cat("插补的零值数:", n_imputed, "\n")
cat("插补比例:", n_imputed / sum(original == 0), "\n")

# 安装: pip install deepimpute
from deepimpute.multinet import MultiNet
import pandas as pd

# 1. 加载数据 (细胞×基因 的DataFrame)
data = pd.read_csv("raw_counts.csv", index_col=0)

# 2. 运行DeepImpute
model = MultiNet(
    n_cores=8,                # 并行核数
    sub_outputdim=512,        # 每个子网络的目标基因数
    hidden_dim=(256,),        # 隐藏层维度
    dropout_rate=0.2,         # dropout正则化
    batch_size=64,            # 批次大小
    max_epochs=100,           # 最大训练轮数
    learning_rate=0.0001      # 学习率
)

# 训练并插补
imputed_data = model.fit(data).predict(data)

# 3. 保存结果
imputed_data.to_csv("deepimpute_result.csv")

# 4. 检查效果
original_zeros = (data == 0).sum().sum()
remaining_zeros = (imputed_data == 0).sum().sum()
print(f"原始零值: {original_zeros}")
print(f"插补后零值: {remaining_zeros}")
print(f"恢复比例: {(original_zeros - remaining_zeros) / original_zeros:.1%}")

def choose_imputation(purpose, data_size, need_conservative):
    if purpose == "visualization":
        return "MAGIC (效果最直观)"
    if purpose == "clustering_de":
        if need_conservative:
            return "scImpute (只补确定的dropout)"
        return "DCA (同时做去噪和降维)"
    if data_size > 50000:
        return "DeepImpute (快速可扩展)"
    return "DCA (综合性能好)"

# === Python环境 ===
pip install magic-impute dca deepimpute scanpy

# === R环境 ===
Rscript -e '
install.packages("scImpute")
BiocManager::install("SingleCellExperiment")
'

# === MAGIC命令行 ===
python -c "
import magic, scanpy as sc
adata = sc.read_h5ad('data.h5ad')
sc.pp.normalize_total(adata); sc.pp.log1p(adata)
magic_op = magic.MAGIC(t='auto')
adata.layers['magic'] = magic_op.fit_transform(adata.to_df()).values
adata.write('data_magic.h5ad')
"

# === DCA命令行 ===
dca raw_counts.h5ad dca_output/ --threads 8

# === 批量对比 ===
python benchmark_imputation.py \
  --input raw_counts.h5ad \
  --methods magic,dca,deepimpute \
  --output benchmark_results/