148_表达矩阵反卷积方法总结¶

# 简单线性模型演示
set.seed(42)
# 模拟3种细胞类型的参考谱 (100基因 × 3类型)
S <- matrix(abs(rnorm(300)), nrow = 100, ncol = 3)
colnames(S) <- c("T_cell", "Macrophage", "Fibroblast")

# 真实比例
true_fractions <- c(0.5, 0.3, 0.2)

# 生成模拟bulk表达 (加噪声)
bulk <- S %*% true_fractions + rnorm(100, sd = 0.1)

# 用非负最小二乘求解
library(nnls)
fit <- nnls(S, bulk)
estimated <- coef(fit) / sum(coef(fit))  # 归一化
names(estimated) <- colnames(S)
print(round(estimated, 3))

# CIBERSORTx 主要通过网页界面操作: https://cibersortx.stanford.edu/
# 但也可以使用Docker版本本地运行

# 1. 准备参考谱：从scRNA-seq生成签名矩阵
# 输入文件格式：基因 × 细胞类型的平均表达矩阵(tab-separated)
# 列名为细胞类型标签

# 2. 使用Docker运行 (命令行方式)
# docker pull cibersortx/fractions
# docker run -v /path/to/data:/src/data \
#   -v /path/to/output:/src/outdir \
#   cibersortx/fractions \
#   --username your_email \
#   --token your_token \
#   --single_cell TRUE \
#   --refsample sc_reference.txt \
#   --mixture bulk_mixture.txt \
#   --perm 100

# 3. 结果解读
results <- read.delim("CIBERSORTx_Results.txt")
head(results)
# 每行一个样本，每列一个细胞类型的比例
# 最后3列: P-value, Correlation, RMSE

# 关键参数
# --perm 100         : 置换检验次数(评估统计显著性)
# --QN TRUE          : 分位数归一化(推荐用于微阵列数据)
# --absolute FALSE   : FALSE=相对比例; TRUE=绝对模式(需要参考)
# --rmbatchSmode TRUE: S-mode批次校正(跨平台时推荐)
# --rmbatchBmode TRUE: B-mode批次校正(参考谱与bulk来自不同技术)

# 安装MuSiC
# devtools::install_github('xuranw/MuSiC')
library(MuSiC)
library(Biobase)

# 1. 准备数据：需要ExpressionSet格式
# 单细胞参考数据
sc.eset <- ExpressionSet(
  assayData = as.matrix(sc_counts),  # 基因 × 细胞 的counts矩阵
  phenoData = AnnotatedDataFrame(sc_meta)  # 需要包含 cellType 和 SubjectName 列
)

# bulk数据
bulk.eset <- ExpressionSet(
  assayData = as.matrix(bulk_counts),
  phenoData = AnnotatedDataFrame(bulk_meta)
)

# 2. 运行MuSiC反卷积
music_result <- music_prop(
  bulk.eset = bulk.eset,
  sc.eset = sc.eset,
  clusters = "cellType",      # 细胞类型列名
  samples = "SubjectName",    # 个体/样本列名
  select.ct = c("T_cell", "B_cell", "Macrophage", "Fibroblast"),  # 选择细胞类型
  verbose = TRUE
)

# 3. 提取结果
estimated_props <- music_result$Est.prop.weighted
head(estimated_props)

# 4. 可视化
library(ggplot2)
library(tidyr)
props_df <- as.data.frame(estimated_props)
props_df$Sample <- rownames(props_df)
props_long <- pivot_longer(props_df, -Sample, names_to = "CellType", values_to = "Proportion")

ggplot(props_long, aes(x = Sample, y = Proportion, fill = CellType)) +
  geom_bar(stat = "identity") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "MuSiC Deconvolution Results")

# 安装
# devtools::install_github('cozygene/bisque')
library(BisqueRNA)
library(Biobase)

# === Reference-based模式 ===
# 准备单细胞和bulk的ExpressionSet
sc.eset <- ExpressionSet(
  assayData = as.matrix(sc_counts),
  phenoData = AnnotatedDataFrame(data.frame(
    cellType = sc_labels,
    SubjectName = sc_subjects,
    row.names = colnames(sc_counts)
  ))
)

bulk.eset <- ExpressionSet(
  assayData = as.matrix(bulk_counts),
  phenoData = AnnotatedDataFrame(bulk_meta)
)

# 运行reference-based反卷积
res_ref <- ReferenceBasedDecomposition(
  bulk.eset = bulk.eset,
  sc.eset = sc.eset,
  cell.types = "cellType",
  subject.names = "SubjectName",
  use.overlap = TRUE  # 使用bulk和sc共有基因
)
print(res_ref$bulk.props)

# === Marker-based模式 ===
# 仅需标记基因列表
markers <- list(
  T_cell = c("CD3D", "CD3E", "CD4", "CD8A"),
  B_cell = c("CD19", "CD79A", "MS4A1"),
  Macrophage = c("CD68", "CD163", "CSF1R")
)

res_marker <- MarkerBasedDecomposition(
  bulk.eset = bulk.eset,
  markers = markers
)
print(res_marker$bulk.props)

# 安装
# devtools::install_github('meichendong/SCDC')
library(SCDC)
library(Biobase)

# 准备多个参考数据集
sc.eset1 <- ExpressionSet(assayData = as.matrix(sc_counts1),
  phenoData = AnnotatedDataFrame(sc_meta1))  # 参考数据集1
sc.eset2 <- ExpressionSet(assayData = as.matrix(sc_counts2),
  phenoData = AnnotatedDataFrame(sc_meta2))  # 参考数据集2

bulk.eset <- ExpressionSet(assayData = as.matrix(bulk_counts),
  phenoData = AnnotatedDataFrame(bulk_meta))

# SCDC集成反卷积
scdc_result <- SCDC_ENSEMBLE(
  bulk.eset = bulk.eset,
  sc.eset.list = list(ref1 = sc.eset1, ref2 = sc.eset2),
  ct.varname = "cellType",
  sample = "SubjectName",
  ct.sub = c("T_cell", "B_cell", "Macrophage", "Fibroblast"),
  search.length = 0.01  # 网格搜索步长
)

# 查看集成权重和比例
print(scdc_result$w_table)      # 各参考数据集的权重
print(scdc_result$prop.est.mvw) # 最终估计比例

# 安装
# devtools::install_github("kkang7/CDSeq_R_Package")
library(CDSeq)

# 输入：bulk counts矩阵 (基因 × 样本)
bulk_matrix <- as.matrix(bulk_counts)

# 运行CDSeq
cdseq_result <- CDSeq(
  bulk_data = bulk_matrix,
  cell_type_number = 5,       # 预估细胞类型数
  mcmc_iterations = 700,      # MCMC迭代次数
  dilution_factor = 1,        # 稀释因子
  block_number = 6,           # 块数(并行)
  gene_length = gene_lengths, # 基因长度(用于TPM校正, 可选)
  reference_gep = NULL        # NULL表示无参考
)

# 结果
estimated_gep <- cdseq_result$estGEP       # 估计的细胞类型表达谱
estimated_prop <- cdseq_result$estProp      # 估计的细胞比例

# 细胞类型注释：将估计的GEP与已知marker比对
# 这一步需要手动或利用已知marker进行注释

# 安装
# install.packages("dtangle")
library(dtangle)

# 输入数据
# pure_samples: 纯化细胞类型的表达矩阵
# mixture: bulk混合样本的表达矩阵
# pure_samples和mixture需要在同一个矩阵中

# 合并数据
Y <- rbind(pure_samples, mixture_samples)

# 定义哪些行是纯参考
pure_indices <- list(
  T_cell = 1:10,         # 行1-10是纯T细胞
  B_cell = 11:20,        # 行11-20是纯B细胞
  Macrophage = 21:30     # 行21-30是纯巨噬细胞
)

# 运行dtangle
dt_result <- dtangle(
  Y = Y,
  pure_samples = pure_indices,
  n_markers = 50,        # 每种细胞类型的标记基因数
  data_type = "rna-seq", # 数据类型
  marker_method = "ratio" # 标记基因选择方法
)

# 查看结果
head(dt_result$estimates)  # 估计的细胞比例
dt_result$markers          # 选择的标记基因

# 方法选择建议逻辑
choose_deconv_method <- function(has_scrna, n_refs, need_gep, n_samples) {
  if (!has_scrna) {
    if (n_samples >= 20) return("CDSeq (无参考)")
    else return("dtangle (标记基因) 或 CIBERSORT (LM22内置谱)")
  }
  if (n_refs > 1) return("SCDC (多参考集成)")
  if (need_gep) return("CIBERSORTx (高分辨率模式)")
  return("MuSiC 或 BisqueRNA (单参考)")
}

# 示例
choose_deconv_method(has_scrna = TRUE, n_refs = 1, need_gep = FALSE, n_samples = 50)
# [1] "MuSiC 或 BisqueRNA (单参考)"

# 从scRNA-seq构建pseudo-bulk (已知比例)
library(Seurat)

# 假设已有Seurat对象 sc_obj, 细胞类型在 sc_obj$celltype
create_pseudobulk <- function(sc_obj, n_samples = 20, n_cells = 500) {
  cell_types <- unique(sc_obj$celltype)
  pseudo_list <- list()
  true_props <- matrix(0, nrow = n_samples, ncol = length(cell_types))
  colnames(true_props) <- cell_types

  for (i in 1:n_samples) {
    # 随机生成比例 (Dirichlet分布)
    props <- rgamma(length(cell_types), shape = 2)
    props <- props / sum(props)
    true_props[i, ] <- props

    # 按比例抽取细胞
    cells_selected <- c()
    for (j in seq_along(cell_types)) {
      ct_cells <- WhichCells(sc_obj, expression = celltype == cell_types[j])
      n_select <- round(props[j] * n_cells)
      if (n_select > 0 & length(ct_cells) > 0) {
        cells_selected <- c(cells_selected, sample(ct_cells, min(n_select, length(ct_cells))))
      }
    }

    # 汇总为pseudo-bulk
    pseudo_list[[i]] <- rowSums(GetAssayData(sc_obj, slot = "counts")[, cells_selected])
  }

  pseudo_bulk <- do.call(cbind, pseudo_list)
  colnames(pseudo_bulk) <- paste0("Pseudo_", 1:n_samples)
  return(list(bulk = pseudo_bulk, true_props = true_props))
}

# 创建pseudo-bulk
pb <- create_pseudobulk(sc_obj, n_samples = 30)

# 用各方法反卷积后与 pb$true_props 比较
# 评估指标：Pearson相关系数、RMSE、Lin's CCC
library(DescTools)
cor(estimated_props, pb$true_props)          # Pearson相关
sqrt(mean((estimated_props - pb$true_props)^2))  # RMSE
CCC(estimated_props[,1], pb$true_props[,1])  # 一致性相关系数

# === 环境准备 ===
# R包安装
Rscript -e '
install.packages(c("nnls", "dtangle", "DescTools"))
if (!require("BiocManager")) install.packages("BiocManager")
BiocManager::install(c("Biobase", "SingleCellExperiment"))
devtools::install_github("xuranw/MuSiC")
devtools::install_github("cozygene/bisque")
devtools::install_github("meichendong/SCDC")
devtools::install_github("kkang7/CDSeq_R_Package")
'

# === CIBERSORTx Docker版本 ===
docker pull cibersortx/fractions
docker run -v $(pwd)/data:/src/data -v $(pwd)/output:/src/outdir \
  cibersortx/fractions \
  --username your@email.com \
  --token YOUR_TOKEN \
  --single_cell TRUE \
  --refsample sc_reference.txt \
  --mixture bulk_mixture.txt \
  --perm 100

# === 批量运行MuSiC反卷积脚本 ===
Rscript run_music.R --sc_data sc_eset.rds \
  --bulk_data bulk_eset.rds \
  --cell_types "T_cell,B_cell,Macrophage" \
  --output music_results.csv

# === Pseudo-bulk基准评估 ===
Rscript benchmark_deconv.R --sc_obj seurat_obj.rds \
  --n_pseudo 50 \
  --methods "MuSiC,BisqueRNA,dtangle" \
  --output benchmark_results.pdf