单细胞细胞通讯分析¶

# CellChatDB 数据库结构
library(CellChat)
CellChatDB <- CellChatDB.human
str(CellChatDB)
# 包含: interaction, complex, cofactor, geneInfo
table(CellChatDB$interaction$annotation)
# Secreted Signaling, ECM-Receptor, Cell-Cell Contact

# 从Seurat对象提取数据用于CellChat
library(Seurat)
library(CellChat)

# 假设已有注释好的Seurat对象 seurat_obj
# 提取归一化表达矩阵
data.input <- GetAssayData(seurat_obj, assay = "RNA", slot = "data")

# 提取元数据
meta <- data.frame(
  labels = seurat_obj$celltype,    # 细胞类型标签
  samples = seurat_obj$sample,      # 样本信息（可选）
  row.names = colnames(seurat_obj)
)

# 创建CellChat对象
cellchat <- createCellChat(object = data.input, meta = meta, group.by = "labels")

# CellPhoneDB 输入准备（Python）
import pandas as pd
import scanpy as sc

# 读取 h5ad 文件
adata = sc.read_h5ad("scRNA_annotated.h5ad")

# 导出 CellPhoneDB 所需文件
# 1. counts矩阵（转置，基因×细胞）
counts = adata.to_df().T
counts.to_csv("cellphonedb_count.txt", sep="\t")

# 2. meta文件（细胞ID + 细胞类型）
meta = pd.DataFrame({
    'Cell': adata.obs.index,
    'cell_type': adata.obs['celltype'].values
})
meta.to_csv("cellphonedb_meta.txt", sep="\t", index=False)

# === CellChat 完整分析流程 ===

# 1. 设置数据库
CellChatDB <- CellChatDB.human
# 可选择子集，如只看分泌型信号
# CellChatDB.use <- subsetDB(CellChatDB, search = "Secreted Signaling")
cellchat@DB <- CellChatDB

# 2. 预处理：识别过表达基因
cellchat <- subsetData(cellchat)
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)

# 3. 推断细胞通讯网络
cellchat <- computeCommunProb(cellchat, type = "triMean")
# type参数：triMean(默认截断均值), truncatedMean, median
# 过滤低细胞数的通讯
cellchat <- filterCommunication(cellchat, min.cells = 10)

# 4. 通路水平通讯推断
cellchat <- computeCommunProbPathway(cellchat)

# 5. 聚合通讯网络
cellchat <- aggregateNet(cellchat)

# 6. 可视化：通讯网络圆形图
groupSize <- as.numeric(table(cellchat@idents))
par(mfrow = c(1, 2), xpd = TRUE)
netVisual_circle(cellchat@net$count, vertex.weight = groupSize,
                 weight.scale = TRUE, label.edge = FALSE,
                 title.name = "Number of interactions")
netVisual_circle(cellchat@net$weight, vertex.weight = groupSize,
                 weight.scale = TRUE, label.edge = FALSE,
                 title.name = "Interaction weights/strength")

# 7. 查看特定通路
netVisual_aggregate(cellchat, signaling = c("WNT"), layout = "circle")
netVisual_aggregate(cellchat, signaling = c("WNT"), layout = "chord")
netVisual_heatmap(cellchat, signaling = c("WNT"), color.heatmap = "Reds")

# 8. 配体-受体对贡献度
netAnalysis_contribution(cellchat, signaling = "WNT")

# 9. 气泡图展示所有显著通讯
netVisual_bubble(cellchat, sources.use = c(1, 2), targets.use = c(3, 4),
                 remove.isolate = FALSE)

# CellPhoneDB v5 命令行运行
# 安装
pip install cellphonedb

# 方法1：统计分析（置换检验）
cellphonedb method statistical_analysis \
    cellphonedb_meta.txt \
    cellphonedb_count.txt \
    --counts-data hgnc_symbol \
    --threads 8 \
    --iterations 1000 \
    --threshold 0.1 \
    --output-path ./cpdb_results

# 方法2：DEG分析（基于差异基因）
cellphonedb method degs_analysis \
    cellphonedb_meta.txt \
    cellphonedb_count.txt \
    degs_file.txt \
    --counts-data hgnc_symbol \
    --threshold 0.1 \
    --output-path ./cpdb_degs_results

# CellPhoneDB 结果解析与可视化
import pandas as pd
from cellphonedb.utils import db_utils

# 读取结果
means = pd.read_csv("cpdb_results/means.txt", sep="\t")
pvalues = pd.read_csv("cpdb_results/pvalues.txt", sep="\t")
significant = pd.read_csv("cpdb_results/significant_means.txt", sep="\t")

# 过滤显著互作 (p < 0.05)
# pvalues的列为细胞类型对（如 T_cells|Macrophages）
celltype_cols = [c for c in pvalues.columns if '|' in c]
sig_mask = pvalues[celltype_cols] < 0.05

# 统计每个细胞类型对的显著互作数
sig_counts = sig_mask.sum()
print(sig_counts.sort_values(ascending=False).head(20))

# CellPhoneDB v5 内置可视化
from cellphonedb.utils import vizutils
vizutils.dot_plot(
    means_path="cpdb_results/means.txt",
    pvalues_path="cpdb_results/pvalues.txt",
    rows=["PDGFA_PDGFRA", "VEGFA_FLT1", "TNF_TNFRSF1A"],
    columns=["Fibroblast|Endothelial", "Macrophage|T_cell"],
    output_path="cpdb_dotplot.pdf"
)

# === NicheNet 分析流程 ===
library(nichenetr)
library(Seurat)
library(tidyverse)

# 加载先验模型
ligand_target_matrix <- readRDS("ligand_target_matrix.rds")
lr_network <- readRDS("lr_network.rds")
weighted_networks <- readRDS("weighted_networks.rds")

# 定义发送和接收细胞
sender_celltypes <- c("Macrophage", "Dendritic_cell")
receiver_celltype <- "CD8_T_cell"

# 获取接收细胞的差异基因作为靶基因集
seurat_obj_receiver <- subset(seurat_obj, celltype == receiver_celltype)
DE_genes <- FindMarkers(seurat_obj, ident.1 = "Disease", ident.2 = "Control",
                         subset.ident = receiver_celltype)
geneset_oi <- rownames(DE_genes %>% filter(p_val_adj < 0.05 & abs(avg_log2FC) > 0.5))

# 获取接收细胞表达的背景基因
background_expressed_genes <- rownames(seurat_obj_receiver)[
  Matrix::rowMeans(GetAssayData(seurat_obj_receiver) > 0) >= 0.10
]

# 获取发送细胞表达的配体
sender_cells <- subset(seurat_obj, celltype %in% sender_celltypes)
expressed_ligands <- rownames(sender_cells)[
  Matrix::rowMeans(GetAssayData(sender_cells) > 0) >= 0.10
]
ligands <- intersect(expressed_ligands, lr_network$from)

# 获取接收细胞表达的受体
expressed_receptors <- rownames(seurat_obj_receiver)[
  Matrix::rowMeans(GetAssayData(seurat_obj_receiver) > 0) >= 0.10
]
receptors <- intersect(expressed_receptors, lr_network$to)

# 过滤合法的配体-受体对
lr_network_filtered <- lr_network %>%
  filter(from %in% ligands & to %in% receptors)
potential_ligands <- unique(lr_network_filtered$from)

# 预测配体活性
ligand_activities <- predict_ligand_activities(
  geneset = geneset_oi,
  background_expressed_genes = background_expressed_genes,
  ligand_target_matrix = ligand_target_matrix,
  potential_ligands = potential_ligands
)

# 排序取top配体
best_upstream_ligands <- ligand_activities %>%
  arrange(-aupr_corrected) %>%
  pull(test_ligand) %>%
  head(30)

# 配体-靶基因热图
active_ligand_target <- get_active_ligand_target_links(
  best_upstream_ligands, geneset_oi, ligand_target_matrix, n = 200
)
vis_ligand_target <- prepare_ligand_target_visualization(
  active_ligand_target, ligand_target_matrix, cutoff = 0.33
)
p <- make_heatmap_ggplot(vis_ligand_target, "Prioritized ligands",
                          "Predicted target genes",
                          legend_title = "Regulatory potential")

# === CellChat 多条件比较 ===

# 分别创建两个条件的CellChat对象
cellchat_ctrl <- createCellChat(data_ctrl, meta = meta_ctrl, group.by = "celltype")
cellchat_disease <- createCellChat(data_disease, meta = meta_disease, group.by = "celltype")

# 分别运行完整分析（同上流程）
# ... 省略中间步骤 ...

# 合并用于比较
object.list <- list(Control = cellchat_ctrl, Disease = cellchat_disease)
cellchat_merged <- mergeCellChat(object.list, add.names = names(object.list))

# 比较互作数量和强度
gg1 <- compareInteractions(cellchat_merged, show.legend = FALSE, group = c(1, 2))
gg2 <- compareInteractions(cellchat_merged, show.legend = FALSE,
                            group = c(1, 2), measure = "weight")
gg1 + gg2

# 差异互作网络
netVisual_diffInteraction(cellchat_merged, weight.scale = TRUE)
netVisual_diffInteraction(cellchat_merged, weight.scale = TRUE, measure = "weight")

# 差异热图
netVisual_heatmap(cellchat_merged)

# 识别上调/下调的信号通路
gg <- rankNet(cellchat_merged, mode = "comparison", stacked = TRUE, do.stat = TRUE)
gg

# === 通讯模式分析（NMF） ===
# 选择模式数（通过Cophenetic和Silhouette指标）
selectK(cellchat, pattern = "outgoing")

# 识别传出通讯模式
cellchat <- identifyCommunicationPatterns(cellchat, pattern = "outgoing", k = 3)
netAnalysis_river(cellchat, pattern = "outgoing")
netAnalysis_dot(cellchat, pattern = "outgoing")

# 识别传入通讯模式
cellchat <- identifyCommunicationPatterns(cellchat, pattern = "incoming", k = 4)
netAnalysis_river(cellchat, pattern = "incoming")

# 网络中心性分析
cellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP")

# 热图展示信号通路的中心性
netAnalysis_signalingRole_network(cellchat, signaling = "WNT",
                                   width = 8, height = 2.5,
                                   font.size = 10)

# 散点图：传出vs传入信号强度
gg <- netAnalysis_signalingRole_scatter(cellchat)
gg

# 识别主要的sender/receiver/mediator/influencer
gg1 <- netAnalysis_signalingRole_heatmap(cellchat, pattern = "outgoing")
gg2 <- netAnalysis_signalingRole_heatmap(cellchat, pattern = "incoming")
gg1 + gg2

# 信号通路相似性聚类
cellchat <- computeNetSimilarity(cellchat, type = "functional")
cellchat <- netEmbedding(cellchat, type = "functional")
cellchat <- netClustering(cellchat, type = "functional")
netVisual_embedding(cellchat, type = "functional", label.size = 3.5)

# === 环境配置 ===
# R包安装
# install.packages("devtools")
# CellChat v2已迁移至新仓库（旧仓库sqjin/CellChat已停止维护）
# devtools::install_github("jinworks/CellChat")
# devtools::install_github("saeyslab/nichenetr")

# Python环境
pip install cellphonedb scanpy

# === CellPhoneDB 命令行 ===
cellphonedb method statistical_analysis meta.txt counts.txt \
    --counts-data hgnc_symbol --threads 8 --output-path ./results

# === CellPhoneDB v5 数据库更新 ===
cellphonedb database download --version v5.0.0

# === LIANA+（多方法整合框架） ===
# LIANA+已发展为Python包（liana-py），是scverse生态的一部分
# pip install liana
# R版liana仍可用：devtools::install_github("saezlab/liana")
library(liana)
liana_results <- liana_wrap(seurat_obj, method = c("cellchat", "natmi", "sca"))
liana_aggregate <- liana_results %>% liana_aggregate()
# 注：LIANA+ Python版（liana-py）功能更全，支持空间数据、多模态、代谢物通讯等