数据清洗与批次效应处理（Data Cleaning & Batch Effect Correction）¶

# ===== R语言：检查缺失值 =====
library(tidyverse)  # 加载数据处理包

# 读取表达矩阵（行=基因，列=样本）
expr <- read.csv("expression_matrix.csv", row.names = 1)

# 查看每列缺失值数量
colSums(is.na(expr))  # 按样本统计NA数量

# 查看缺失比例
missing_rate <- colMeans(is.na(expr))  # 每个样本的缺失比例
print(missing_rate)

# 策略1：删除缺失值过多的样本（>20%缺失就删）
expr_clean <- expr[, missing_rate < 0.2]

# 策略2：用中位数填充（适合表达数据）
expr_filled <- expr %>%
  mutate(across(everything(), ~ifelse(is.na(.), median(., na.rm = TRUE), .)))
  # across(everything()) = 对每一列操作
  # ifelse(is.na(.), median(.), .) = 如果是NA就用中位数替代，否则保持原值

# 策略3：用KNN填充（更智能，考虑相邻样本的值）
library(impute)  # Bioconductor的填充包
expr_knn <- impute.knn(as.matrix(expr), k = 10)$data
# k=10 表示用最近的10个基因的值来估算缺失值

# ===== 异常值检测方法 =====

# 方法1：箱线图法（IQR法则）
detect_outliers_iqr <- function(x) {
  Q1 <- quantile(x, 0.25, na.rm = TRUE)  # 第一四分位数
  Q3 <- quantile(x, 0.75, na.rm = TRUE)  # 第三四分位数
  IQR <- Q3 - Q1                          # 四分位距
  lower <- Q1 - 1.5 * IQR                 # 下界
  upper <- Q3 + 1.5 * IQR                 # 上界
  return(x < lower | x > upper)           # 返回TRUE/FALSE标记
}

# 方法2：Z-score法（适合正态分布数据）
detect_outliers_zscore <- function(x, threshold = 3) {
  z <- (x - mean(x, na.rm = TRUE)) / sd(x, na.rm = TRUE)
  # z = (每个值 - 均值) / 标准差
  return(abs(z) > threshold)  # |z| > 3 视为异常
}

# 方法3：MAD法（中位数绝对偏差，对偏态数据更稳健）
detect_outliers_mad <- function(x, threshold = 3) {
  med <- median(x, na.rm = TRUE)
  mad_val <- mad(x, na.rm = TRUE)  # MAD = median(|x - median(x)|) * 1.4826
  modified_z <- 0.6745 * (x - med) / mad_val
  return(abs(modified_z) > threshold)
}

# ===== 重复样本/重复基因检查 =====

# 检查是否有重复的样本名
duplicated_samples <- duplicated(colnames(expr))
cat("重复样本数：", sum(duplicated_samples), "\n")

# 检查是否有重复的基因名
duplicated_genes <- duplicated(rownames(expr))
cat("重复基因数：", sum(duplicated_genes), "\n")

# 处理重复基因：取表达量最大的那个（常见策略）
if (sum(duplicated_genes) > 0) {
  expr$gene <- rownames(expr)
  expr$max_expr <- rowMeans(expr[, -ncol(expr)])  # 计算每行均值
  expr <- expr %>%
    group_by(gene) %>%           # 按基因名分组
    slice_max(max_expr, n = 1) %>%  # 每组保留均值最大的一行
    ungroup()
  rownames(expr) <- expr$gene
  expr <- expr %>% select(-gene, -max_expr)
}

# ===== 数据类型验证 =====

# 检查表达矩阵是否全是数值
str(expr)  # 查看每列的数据类型

# 强制转换为数值型（有时读入的数据会变成字符型）
expr_numeric <- apply(expr, 2, as.numeric)
# apply(矩阵, 2, 函数) = 对每一列应用函数
# 2 表示按列操作（1表示按行）

# 检查转换后是否产生了新的NA（说明原数据有非数字字符）
new_na <- sum(is.na(expr_numeric)) - sum(is.na(expr))
if (new_na > 0) {
  warning(paste("转换产生了", new_na, "个新的NA，请检查原始数据"))
}

# ===== 常见标准化方法 =====

# 方法1：log2转换（最常用，压缩数据范围）
expr_log <- log2(expr + 1)
# +1 防止 log(0) = -Inf

# 方法2：Z-score标准化（每个基因均值=0，标准差=1）
expr_zscore <- t(scale(t(expr)))
# scale() 默认对列标准化，所以先转置再转回来
# 效果：让不同基因的表达量可以直接比较

# 方法3：分位数标准化（quantile normalization，让样本间分布一致）
library(preprocessCore)  # Bioconductor包
expr_qnorm <- normalize.quantiles(as.matrix(expr))
# 原理：让每个样本的表达值分布完全相同

# 方法4：TPM/FPKM/CPM（RNA-seq专用，考虑基因长度和测序深度）
# CPM = Counts Per Million
cpm <- t(t(expr) / colSums(expr) * 1e6)
# 每个值 / 该样本总reads数 × 100万

# ===== 检测样本混淆 =====

# 方法1：性别检查（如果有性别信息）
# XIST基因在女性高表达，Y染色体基因在男性高表达
xist_expr <- expr["XIST", ]  # 提取XIST基因的表达
rps4y1_expr <- expr["RPS4Y1", ]  # 提取Y染色体基因

# 画散点图，正常情况下应该分成两群（男/女）
plot(as.numeric(xist_expr), as.numeric(rps4y1_expr),
     xlab = "XIST (女性标志)", ylab = "RPS4Y1 (男性标志)",
     main = "性别检查：发现样本混淆")
# 如果有样本落在"错误"的群里，说明标签可能贴反了

# 方法2：SNP指纹检测（更可靠）
# 用样本的遗传变异谱做"身份证"比对

# ===== 低质量样本筛查 =====

# 指标1：文库大小（总reads数太少 = 测序深度不够）
lib_size <- colSums(expr)  # 每个样本的总reads数
low_depth <- lib_size < quantile(lib_size, 0.05)  # 低于5%分位数

# 指标2：检测到的基因数太少
genes_detected <- colSums(expr > 0)  # 每个样本表达>0的基因数
low_genes <- genes_detected < quantile(genes_detected, 0.05)

# 指标3：线粒体基因比例过高（常用于单细胞，比例高=细胞破损）
mt_genes <- grep("^MT-", rownames(expr), value = TRUE)  # 找线粒体基因
mt_pct <- colSums(expr[mt_genes, ]) / colSums(expr) * 100
high_mt <- mt_pct > 20  # 线粒体比例>20%视为低质量

# 综合标记低质量样本
bad_samples <- low_depth | low_genes | high_mt
cat("低质量样本数：", sum(bad_samples), "\n")

# ===== GC偏差检查与校正 =====

# 检查：按GC含量分组，看表达量是否均匀
# 通常用EDASeq或cqn包校正

library(EDASeq)  # Bioconductor的GC校正包
# 获取基因的GC含量信息
# gene_gc <- getGeneLengthAndGCContent(rownames(expr), "hsa")

# 校正GC偏差
# expr_gc_corrected <- withinLaneNormalization(expr, gene_gc$gc, which = "full")
# withinLaneNormalization = 在同一个样本内，消除GC含量对表达量的影响

# ===== PCA着色检测批次效应 =====
library(ggplot2)

# 准备数据：log转换 + 标准化
expr_log <- log2(expr + 1)

# 做PCA
pca_result <- prcomp(t(expr_log), scale. = TRUE)
# prcomp = 主成分分析函数
# t() = 转置（PCA需要样本在行、基因在列）
# scale. = TRUE 表示先标准化

# 提取PC1和PC2
pca_df <- data.frame(
  PC1 = pca_result$x[, 1],   # 第一主成分
  PC2 = pca_result$x[, 2],   # 第二主成分
  Batch = metadata$batch,     # 批次信息
  Group = metadata$condition  # 生物学分组（如T2D vs 健康）
)

# 按批次着色
ggplot(pca_df, aes(x = PC1, y = PC2, color = Batch, shape = Group)) +
  geom_point(size = 3) +  # 画散点
  theme_minimal() +         # 简洁主题
  labs(title = "PCA: 按批次着色 (校正前)")
# 如果样本主要按批次聚类而非按生物学分组聚类 → 有批次效应
# 理想情况：同一生物组的样本聚在一起，不管它们来自哪个批次

# ===== MDS检测批次效应 =====
library(limma)

# plotMDS是limma包自带的快速MDS可视化
plotMDS(expr_log,
        col = as.numeric(factor(metadata$batch)),  # 按批次着色
        pch = 16,  # 实心圆点
        main = "MDS: 批次效应检测")
legend("topright", levels(factor(metadata$batch)),
       col = 1:nlevels(factor(metadata$batch)), pch = 16)
# MDS和PCA类似，但基于样本间距离而非方差分解

# ===== kBET：定量检测批次效应 =====
# kBET = k-nearest neighbour Batch Effect Test
# 原理：检查每个样本的k个最近邻里，各批次的比例是否符合总体比例

library(kBET)  # 需要从GitHub安装

# 运行kBET
kbet_result <- kBET(
  df = t(expr_log),          # 表达矩阵（样本在行）
  batch = metadata$batch,     # 批次标签
  k0 = 25                     # k值（检查最近的25个邻居）
)

# 查看结果
cat("kBET拒绝率：", kbet_result$summary$kBET.observed[1], "\n")
# 拒绝率高 = 批次效应严重
# 拒绝率接近0 = 批次混合良好

# ===== LISI：评估批次混合程度 =====
# LISI = Local Inverse Simpson Index
# 由Harmony的作者提出，越高表示混合越好

library(lisi)  # 安装：install.packages("lisi")

# 计算LISI
lisi_result <- compute_lisi(
  X = pca_result$x[, 1:20],     # 用PCA前20个主成分
  meta_data = metadata,           # 元数据
  label_colnames = c("batch")     # 要评估的标签列
)

cat("平均LISI：", mean(lisi_result$batch), "\n")
# LISI = 1 → 邻居全来自同一批次（批次效应严重）
# LISI = 批次数 → 邻居来自所有批次（完美混合）

# ===== ComBat校正 =====
library(sva)  # 安装：BiocManager::install("sva")

# 准备模型矩阵（保留你关心的生物学差异）
mod <- model.matrix(~condition, data = metadata)
# ~condition 表示"保留condition（如T2D vs 健康）的差异"
# 这很关键：不加mod，ComBat可能把生物学差异也消掉

# 运行ComBat
expr_combat <- ComBat(
  dat = as.matrix(expr_log),  # 表达矩阵（已log转换）
  batch = metadata$batch,      # 批次标签
  mod = mod                    # 模型矩阵（保留生物学差异）
)
# ComBat的原理：
# 1. 估计每个批次的均值偏移和方差缩放
# 2. 用贝叶斯方法"借力"其他基因的信息，让估计更稳定
# 3. 把每个批次的数据调整到统一水平

# 校正后再做PCA验证
pca_after <- prcomp(t(expr_combat), scale. = TRUE)
# 重新画PCA图，确认样本不再按批次聚类

# ===== limma::removeBatchEffect =====
library(limma)  # Bioconductor核心包

# removeBatchEffect用线性回归去除批次效应
expr_limma <- removeBatchEffect(
  x = expr_log,                           # 表达矩阵
  batch = metadata$batch,                  # 批次变量
  design = model.matrix(~condition, data = metadata)
  # design矩阵保护生物学差异不被消除
)

# 注意：removeBatchEffect产生的校正数据仅用于可视化和探索
# 做差异分析时，应在limma的线性模型里直接加入batch作为协变量：
# design <- model.matrix(~condition + batch, data = metadata)
# 这样更统计正确

# ===== Harmony批次校正（单细胞） =====
library(Seurat)   # 单细胞分析框架
library(harmony)  # 安装：install.packages("harmony")

# 假设已有Seurat对象 seurat_obj，包含多个批次的单细胞数据
# 先做标准预处理
seurat_obj <- NormalizeData(seurat_obj)        # 归一化
seurat_obj <- FindVariableFeatures(seurat_obj) # 找高变基因
seurat_obj <- ScaleData(seurat_obj)            # 标准化
seurat_obj <- RunPCA(seurat_obj)               # 降维

# 运行Harmony
seurat_obj <- RunHarmony(
  seurat_obj,
  group.by.vars = "batch",  # 按哪个变量校正
  dims.use = 1:30           # 用前30个PC
)
# Harmony的原理：
# 1. 在PCA空间中迭代聚类
# 2. 让每个聚类内的批次比例趋于一致
# 3. 调整PCA坐标，输出校正后的embedding

# 用Harmony的结果做后续分析
seurat_obj <- RunUMAP(seurat_obj, reduction = "harmony", dims = 1:30)
seurat_obj <- FindNeighbors(seurat_obj, reduction = "harmony", dims = 1:30)
seurat_obj <- FindClusters(seurat_obj, resolution = 0.5)

# 可视化对比
DimPlot(seurat_obj, group.by = "batch")  # 按批次着色
DimPlot(seurat_obj, group.by = "celltype")  # 按细胞类型着色

# ===== MNN校正（单细胞） =====
library(batchelor)  # Bioconductor包，BiocManager::install("batchelor")

# 假设有两个批次的SingleCellExperiment对象
# batch1_sce, batch2_sce

# 运行MNN
mnn_result <- fastMNN(
  batch1_sce,
  batch2_sce,
  d = 50  # 保留前50个主成分
)
# fastMNN比原始MNN快很多
# 原理：
# 1. 在两个批次中找互近邻对（A的最近邻是B中的某个细胞，反之亦然）
# 2. 用这些配对估计批次间的偏移向量
# 3. 校正PCA空间中的坐标

# 提取校正后的低维表示
corrected_pcs <- reducedDim(mnn_result, "corrected")

# ===== scVI批次校正（Python） =====
import scvi  # pip install scvi-tools
import scanpy as sc  # pip install scanpy

# 读取数据
adata = sc.read_h5ad("single_cell_data.h5ad")

# 预处理
sc.pp.highly_variable_genes(adata, n_top_genes=2000)
# 选择2000个高变基因用于建模

# 设置scVI模型
scvi.model.SCVI.setup_anndata(
    adata,
    batch_key="batch"  # 指定批次变量
)

# 训练模型
model = scvi.model.SCVI(
    adata,
    n_latent=30,     # 潜在空间维度
    n_layers=2,      # 神经网络层数
    gene_likelihood="nb"  # 负二项分布（适合count数据）
)
model.train(max_epochs=200)  # 训练200轮

# 获取校正后的潜在表示
adata.obsm["X_scVI"] = model.get_latent_representation()
# 这个潜在表示已经消除了批次效应

# 用scVI的结果做聚类和可视化
sc.pp.neighbors(adata, use_rep="X_scVI")  # 基于scVI结果找邻居
sc.tl.umap(adata)  # UMAP降维
sc.tl.leiden(adata)  # Leiden聚类

# 可视化
sc.pl.umap(adata, color=["batch", "cell_type"])

数据到手
   ↓
PCA着色(按batch) → 无明显批次分离 → 不需要校正，直接分析
   ↓ 有批次效应
判断混杂程度
   ↓
batch和group不混杂 → ComBat / removeBatchEffect / 加协变量
   ↓
batch和group部分混杂 → 加协变量 + 验证结论稳健性
   ↓
batch和group完全混杂 → 无法可靠分离，需要补实验