185_单细胞数据质控标准¶

library(Seurat)
library(ggplot2)
library(patchwork)

# 加载数据
data <- Read10X(data.dir = "filtered_feature_bc_matrix/")
seurat_obj <- CreateSeuratObject(counts = data, project = "sample1",
                                  min.cells = 3, min.features = 200)

# 计算QC指标
# 线粒体基因比例
seurat_obj[["percent.mt"]] <- PercentageFeatureSet(seurat_obj, pattern = "^MT-")

# 核糖体基因比例
seurat_obj[["percent.ribo"]] <- PercentageFeatureSet(seurat_obj, pattern = "^RP[SL]")

# 血红蛋白基因比例（血液样本）
seurat_obj[["percent.hb"]] <- PercentageFeatureSet(seurat_obj, pattern = "^HB[^(P)]")

# 计算复杂度（log10GenesPerUMI）
seurat_obj[["log10GenesPerUMI"]] <- log10(seurat_obj$nFeature_RNA) / log10(seurat_obj$nCount_RNA)

# 查看QC指标分布
head(seurat_obj@meta.data[, c("nFeature_RNA", "nCount_RNA", "percent.mt",
                                "percent.ribo", "log10GenesPerUMI")])

# 小提琴图
p1 <- VlnPlot(seurat_obj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
               ncol = 3, pt.size = 0)
print(p1)

# 散点图：UMI vs 基因数（应该呈正相关）
p2 <- FeatureScatter(seurat_obj, feature1 = "nCount_RNA", feature2 = "nFeature_RNA") +
  geom_smooth(method = "lm", se = FALSE) +
  NoLegend()

# 散点图：UMI vs 线粒体比例
p3 <- FeatureScatter(seurat_obj, feature1 = "nCount_RNA", feature2 = "percent.mt") +
  NoLegend()

p2 | p3

# 密度图
library(ggplot2)
ggplot(seurat_obj@meta.data, aes(x = nFeature_RNA)) +
  geom_density(fill = "steelblue", alpha = 0.5) +
  geom_vline(xintercept = c(200, 5000), linetype = "dashed", color = "red") +
  theme_minimal() +
  labs(title = "Gene Count Distribution")

ggplot(seurat_obj@meta.data, aes(x = percent.mt)) +
  geom_density(fill = "salmon", alpha = 0.5) +
  geom_vline(xintercept = 20, linetype = "dashed", color = "red") +
  theme_minimal() +
  labs(title = "Mitochondrial % Distribution")

# 方法1：基于MAD的自适应阈值
calc_mad_threshold <- function(values, nmads = 3, type = "both") {
  med <- median(values)
  mad_val <- mad(values)
  lower <- med - nmads * mad_val
  upper <- med + nmads * mad_val
  if (type == "lower") return(lower)
  if (type == "upper") return(upper)
  return(c(lower = lower, upper = upper))
}

# 计算各指标阈值
nfeature_lower <- calc_mad_threshold(log10(seurat_obj$nFeature_RNA), 3, "lower")
nfeature_upper <- calc_mad_threshold(log10(seurat_obj$nFeature_RNA), 3, "upper")
mt_upper <- calc_mad_threshold(seurat_obj$percent.mt, 3, "upper")

cat("nFeature阈值:", 10^nfeature_lower, "-", 10^nfeature_upper, "\n")
cat("percent.mt上限:", mt_upper, "\n")

# 方法2：scater包的自动QC
library(scater)
sce <- as.SingleCellExperiment(seurat_obj)
sce <- addPerCellQCMetrics(sce, subsets = list(mito = grep("^MT-", rownames(sce))))
qc_results <- quickPerCellQC(colData(sce),
                              sub.fields = "subsets_mito_percent",
                              nmads = 3)
cat("自动过滤的细胞数:", sum(qc_results$discard), "\n")

# 方法3：常用的经验阈值（参考值）
# nFeature_RNA: 200-5000 (或 200-8000 for complex tissues)
# nCount_RNA: 500-50000
# percent.mt: <10% (default) 或 <20% (for certain tissues)
# log10GenesPerUMI: >0.8

# 记录过滤前信息
cat("过滤前细胞数:", ncol(seurat_obj), "\n")

# 过滤
seurat_obj <- subset(seurat_obj,
                      subset = nFeature_RNA > 200 &
                        nFeature_RNA < 5000 &
                        nCount_RNA > 500 &
                        nCount_RNA < 50000 &
                        percent.mt < 15 &
                        log10GenesPerUMI > 0.8)

cat("过滤后细胞数:", ncol(seurat_obj), "\n")

# 过滤后重新可视化确认
VlnPlot(seurat_obj, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
        ncol = 3, pt.size = 0)

# 使用SoupX去除ambient RNA
library(SoupX)

# 加载数据（需要raw和filtered矩阵）
toc <- Read10X("filtered_feature_bc_matrix/")
tod <- Read10X("raw_feature_bc_matrix/")

# 创建SoupChannel
sc <- SoupChannel(tod, toc)

# 需要先做聚类
seurat_temp <- CreateSeuratObject(counts = toc)
seurat_temp <- NormalizeData(seurat_temp) %>%
  FindVariableFeatures() %>%
  ScaleData() %>%
  RunPCA() %>%
  FindNeighbors(dims = 1:20) %>%
  FindClusters(resolution = 0.5)

sc <- setClusters(sc, setNames(seurat_temp$seurat_clusters,
                                 colnames(seurat_temp)))

# 估计污染比例
sc <- autoEstCont(sc)
cat("估计的污染比例:", sc$fit$rhoEst, "\n")

# 校正
corrected_counts <- adjustCounts(sc)

# 使用校正后的数据
seurat_obj <- CreateSeuratObject(counts = corrected_counts)

import scanpy as sc
import numpy as np

adata = sc.read_10x_h5("filtered_feature_bc_matrix.h5")

# 计算QC指标
adata.var['mt'] = adata.var_names.str.startswith('MT-')
adata.var['ribo'] = adata.var_names.str.startswith(('RPS', 'RPL'))
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt', 'ribo'], percent_top=None,
                            log1p=False, inplace=True)

# 可视化
sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],
             jitter=0.4, multi_panel=True, save='_qc.png')

sc.pl.scatter(adata, x='total_counts', y='n_genes_by_counts',
              color='pct_counts_mt', save='_qc_scatter.png')

# MAD过滤
def mad_filter(adata, metric, nmads=3):
    med = np.median(adata.obs[metric])
    mad = np.median(np.abs(adata.obs[metric] - med))
    upper = med + nmads * 1.4826 * mad
    lower = med - nmads * 1.4826 * mad
    return lower, upper

nf_low, nf_high = mad_filter(adata, 'n_genes_by_counts')
mt_low, mt_high = mad_filter(adata, 'pct_counts_mt')

print(f"Gene filter: {nf_low:.0f} - {nf_high:.0f}")
print(f"MT filter: < {mt_high:.1f}%")

# 过滤
print(f"Before filtering: {adata.n_obs} cells")
sc.pp.filter_cells(adata, min_genes=200)
adata = adata[adata.obs.n_genes_by_counts < 5000, :]
adata = adata[adata.obs.pct_counts_mt < 15, :]
print(f"After filtering: {adata.n_obs} cells")

# Seurat快速QC
obj[["percent.mt"]] <- PercentageFeatureSet(obj, pattern = "^MT-")
obj <- subset(obj, nFeature_RNA > 200 & nFeature_RNA < 5000 & percent.mt < 15)

# 自适应MAD阈值
library(scater); qc <- quickPerCellQC(colData(sce), sub.fields = "subsets_mito_percent")

# SoupX ambient RNA去除
sc <- SoupChannel(raw, filtered); sc <- autoEstCont(sc); out <- adjustCounts(sc)