180_单细胞doublet检测¶

import scrublet as scr
import scanpy as sc
import numpy as np
import matplotlib.pyplot as plt

# 加载数据
adata = sc.read_10x_h5("filtered_feature_bc_matrix.h5")

# 运行Scrublet
scrub = scr.Scrublet(adata.X, expected_doublet_rate=0.06)
doublet_scores, predicted_doublets = scrub.scrub_doublets(
    min_counts=2,
    min_cells=3,
    min_gene_variability_pctl=85,
    n_prin_comps=30
)

# 查看结果
print(f"检测到的doublet数: {sum(predicted_doublets)}")
print(f"Doublet比例: {sum(predicted_doublets)/len(predicted_doublets)*100:.1f}%")

# 可视化doublet分数分布
scrub.plot_histogram()
plt.savefig("scrublet_histogram.png", dpi=300)

# 将结果添加到adata
adata.obs['doublet_score'] = doublet_scores
adata.obs['predicted_doublet'] = predicted_doublets

# 如果自动阈值不合适，手动调整
# scrub.call_doublets(threshold=0.25)

# UMAP上展示doublet
sc.pp.normalize_total(adata)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata)
sc.tl.pca(adata)
sc.pp.neighbors(adata)
sc.tl.umap(adata)

sc.pl.umap(adata, color=['doublet_score', 'predicted_doublet'], 
           save='_doublets.png')

# 过滤doublet
adata_clean = adata[~adata.obs['predicted_doublet']].copy()
print(f"过滤前: {adata.n_obs} cells, 过滤后: {adata_clean.n_obs} cells")

library(Seurat)
library(DoubletFinder)

# 加载并预处理数据
pbmc <- Read10X(data.dir = "filtered_feature_bc_matrix/")
pbmc <- CreateSeuratObject(counts = pbmc, min.cells = 3, min.features = 200)
pbmc <- NormalizeData(pbmc)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
pbmc <- ScaleData(pbmc)
pbmc <- RunPCA(pbmc)
pbmc <- RunUMAP(pbmc, dims = 1:20)
pbmc <- FindNeighbors(pbmc, dims = 1:20)
pbmc <- FindClusters(pbmc, resolution = 0.5)

# Step1: 优化pK参数
sweep.res <- paramSweep(pbmc, PCs = 1:20, sct = FALSE)
sweep.stats <- summarizeSweep(sweep.res, GT = FALSE)
bcmvn <- find.pK(sweep.stats)

# 选择最优pK
optimal_pk <- as.numeric(as.character(bcmvn$pK[which.max(bcmvn$BCmetric)]))
cat("Optimal pK:", optimal_pk, "\n")

# Step2: 估计同型doublet比例
homotypic.prop <- modelHomotypic(pbmc$seurat_clusters)
nExp_poi <- round(0.08 * nrow(pbmc@meta.data))  # 假设8% doublet率
nExp_poi.adj <- round(nExp_poi * (1 - homotypic.prop))

# Step3: 运行DoubletFinder
pbmc <- doubletFinder(
  pbmc,
  PCs = 1:20,
  pN = 0.25,
  pK = optimal_pk,
  nExp = nExp_poi.adj,
  sct = FALSE
)

# 查看结果
df_col <- grep("DF.classifications", colnames(pbmc@meta.data), value = TRUE)
table(pbmc@meta.data[[df_col]])

# 可视化
DimPlot(pbmc, group.by = df_col, cols = c("Singlet" = "grey", "Doublet" = "red"))

# 过滤doublet
pbmc_clean <- subset(pbmc, cells = colnames(pbmc)[pbmc@meta.data[[df_col]] == "Singlet"])

library(scDblFinder)
library(SingleCellExperiment)

# 将Seurat对象转为SCE
sce <- as.SingleCellExperiment(pbmc)

# 运行scDblFinder
set.seed(42)
sce <- scDblFinder(sce, clusters = TRUE)

# 查看结果
table(sce$scDblFinder.class)
# singlet doublet
#   8500     680

# 将结果转回Seurat
pbmc$scDblFinder_class <- sce$scDblFinder.class
pbmc$scDblFinder_score <- sce$scDblFinder.score

# 可视化
FeaturePlot(pbmc, features = "scDblFinder_score")
DimPlot(pbmc, group.by = "scDblFinder_class")

# 过滤
pbmc_clean <- subset(pbmc, scDblFinder_class == "singlet")

# 使用Demuxlet（基于SNP的方法）
# 需要参考VCF文件
demuxlet \
    --sam possorted_bam.bam \
    --vcf donors.vcf \
    --field GP \
    --out demuxlet_result

# 或使用Vireo（更灵活）
cellsnp-lite \
    -s possorted_bam.bam \
    -b barcodes.tsv \
    -O cellsnp_output \
    -R genome1K.phase3.SNP_AF5e2.chr1toX.hg38.vcf.gz \
    -p 20 \
    --minMAF 0.1 \
    --minCOUNT 20

vireo -c cellsnp_output \
    -N 4 \
    -o vireo_output \
    --randSeed 42

# Scrublet一行运行
import scrublet as scr
scores, preds = scr.Scrublet(adata.X, expected_doublet_rate=0.06).scrub_doublets()

# Scanpy内置（基于Scrublet）
sc.pp.scrublet(adata, expected_doublet_rate=0.06)

# scDblFinder一行运行
sce <- scDblFinder(as.SingleCellExperiment(seurat_obj))

# 批量处理多个样本
library(scDblFinder)
sample_list <- SplitObject(merged_seurat, split.by = "sample")
for (s in names(sample_list)) {
  sce <- as.SingleCellExperiment(sample_list[[s]])
  sce <- scDblFinder(sce)
  sample_list[[s]]$doublet <- sce$scDblFinder.class
}