<think>
我们拥有immucellAI获得的免疫细胞丰度数据，现在需要比较不同组之间（例如不同疾病组、治疗组等）的免疫细胞丰度是否存在显著性差异。
通常，这样的数据是连续变量（细胞比例或丰度），且可能不满足正态分布，因此非参数检验可能更合适。
我们将使用以下步骤：
1. 数据准备：读取数据，确保数据结构正确（行是样本，列是细胞类型和分组信息）
2. 正态性检验（可选）和方差齐性检验（可选）来决定使用参数检验还是非参数检验
3. 根据分组情况选择检验方法：
- 两组比较：Wilcoxon秩和检验（Mann-Whitney U检验）
- 多组比较：Kruskal-Wallis检验，如果显著则进行事后两两比较（如Dunn检验，并控制多重检验）
4. 结果整理和可视化
注意：由于免疫细胞丰度数据通常是成分数据（总和为100%），但immucellAI输出的应该是每种细胞类型的比例（或绝对分数），我们可以直接进行组间比较。
然而，成分数据分析有特殊方法（如ALR/CLR变换后使用参数检验），但这里我们按照常规的组间差异分析进行，因为目前文献中多数使用非参数检验。
我们假设数据格式如下：
样本ID 分组 CD4T CD8T B细胞 NK细胞 ...
sample1 Group1 0.25 0.18 0.10 0.05 ...
...
步骤：
1. 读取数据
2. 按分组进行循环，对每种细胞类型进行组间差异检验
由于可能有多种细胞类型，我们需要进行多重检验校正（如FDR校正）。
我们将编写一个函数来处理多组比较，并返回校正后的p值。
另外，我们也可以考虑使用线性模型（考虑协变量）或混合效应模型（如果有配对设计等），但这里用户只要求组间比较，因此先做简单的组间比较。
如果数据是配对样本（例如同一个病人治疗前后），那么需要使用配对检验（如Wilcoxon符号秩检验）。
根据用户提供的信息，我们不知道具体分组情况，因此我们将提供两种代码：独立样本和配对样本。
我们首先按照独立样本处理，然后补充配对样本的代码。
此外，我们还可以进行可视化，如箱线图展示。
代码结构：
1. 加载必要的包
2. 读入数据
3. 检查数据结构
4. 定义分组列和细胞类型列
5. 对于每种细胞类型：
- 进行正态性检验（Shapiro-Wilk检验）和方差齐性检验（Levene检验或Bartlett检验，Bartlett要求正态，Levene不要求）
- 根据检验结果选择参数检验（t检验或ANOVA）或非参数检验（Wilcoxon或Kruskal-Wallis）
- 但注意：由于细胞比例数据通常不满足正态分布，我们默认使用非参数检验
6. 输出结果表格（包含细胞类型，检验统计量，p值，校正后的p值）
7. 可视化：生成每个细胞类型的箱线图，并在图中标注显著性
我们这里默认使用非参数检验，因为更稳健。
由于用户没有提供具体数据，我们生成一个模拟数据并演示。
假设我们有三组：Group1, Group2, Group3，比较每种细胞类型在三组间是否有差异。
独立样本非参数多组比较：Kruskal-Wallis检验，然后进行Dunn事后检验（并校正p值）。
如果是两组，则用Wilcoxon秩和检验。
我们编写一个通用的函数，可以自动判断分组数量。
注意：分组列必须是因子类型。
我们将使用以下包：
- dplyr: 数据处理
- rstatix: 提供易于使用的检验函数，并能与管道操作符结合
- ggpubr: 用于统计可视化
如果没有安装，请先安装。
代码示例：
步骤1：安装和加载包
步骤2：读取数据
步骤3：将数据从宽格式转换为长格式（便于循环或分组操作）
步骤4：按细胞类型分组，然后进行检验
步骤5：校正p值
步骤6：可视化
由于rstatix包提供了很好的管道操作，我们使用它。
注意：rstatix中的函数要求数据是长格式。
例如：
ID Group CellType Abundance
sample1 Group1 CD4T 0.25
sample1 Group1 CD8T 0.18
...
但是我们的原始数据可能是宽格式，即每种细胞类型占一列。因此需要转换。
我们将使用tidyr的gather或pivot_longer函数。
下面我们开始编写代码。
</think>
以下是根据 ImmucellAI 输出的免疫细胞丰度数据进行组间显著性差异分析的 R 语言完整代码方案，结合微生物组和免疫细胞分析的最佳实践：
### 一、数据预处理与加载
```r
# 加载必要包
library(tidyverse) # 数据处理
library(ggpubr) # 统计可视化
library(rstatix) # 统计检验
library(FSA) # Dunn事后检验
# 读取ImmucellAI输出数据（示例格式）
# 列：SampleID | Group | CD4_Tcells | CD8_Tcells | Bcells | NKcells | ...
cell_data <- read.csv("immucellAI_results.csv") %>%
mutate(Group = factor(Group)) # 确保分组为因子
# 查看数据结构
head(cell_data)
```
### 二、显著性差异分析（自动适配组数）
```r
# 函数：自动选择检验方法
perform_cell_test <- function(data, cell_type) {
cell_vector <- data[[cell_type]]
# 正态性检验（Shapiro-Wilk）
norm_test <- by(cell_vector, data$Group, shapiro.test)
norm_pvals <- sapply(norm_test, `[[`, "p.value")
# 方差齐性检验（Levene）
var_test <- levene_test(data, reformulate("Group", cell_type))
# 根据检验结果选择分析方法
if (all(norm_pvals > 0.05) & var_test$p > 0.05) {
# 参数检验：ANOVA或t检验
if (nlevels(data$Group) > 2) {
anova_result <- aov(reformulate("Group", cell_type), data)
return(list(
test = "ANOVA",
result = tidy(TukeyHSD(anova_result)), # 事后检验
p.value = summary(anova_result)[[1]]$`Pr(>F)`[1]
))
} else {
t_result <- t.test(reformulate("Group", cell_type), data)
return(list(
test = "t-test",
result = tidy(t_result),
p.value = t_result$p.value
))
}
} else {
# 非参数检验：Kruskal-Wallis或Wilcoxon
if (nlevels(data$Group) > 2) {
kw_result <- kruskal.test(reformulate("Group", cell_type), data)
dunn_result <- dunnTest(cell_vector ~ data$Group, method = "bh") # Benjamini-Hochberg校正
return(list(
test = "Kruskal-Wallis",
result = dunn_result$res,
p.value = kw_result$p.value
))
} else {
wilcox_result <- wilcox.test(reformulate("Group", cell_type), data)
return(list(
test = "Wilcoxon",
result = tidy(wilcox_result),
p.value = wilcox_result$p.value
))
}
}
}
# 循环分析所有免疫细胞类型
cell_types <- setdiff(names(cell_data), c("SampleID", "Group"))
results_list <- map(cell_types, ~perform_cell_test(cell_data, .x)) %>%
set_names(cell_types)
# 提取简化结果表格
results_table <- map_dfr(results_list, ~{
tibble(
CellType = .y,
Test = .x$test,
Global.p = .x$p.value,
Details = list(.x$result)
)
}, .id = "CellType")
```
### 三、结果可视化
```r
# 1. 整体差异热图
heatmap_plot <- cell_data %>%
pivot_longer(cols = -c(SampleID, Group), names_to = "CellType", values_to = "Abundance") %>%
group_by(CellType, Group) %>%
summarise(MeanAbundance = mean(Abundance)) %>%
ggplot(aes(x = Group, y = CellType, fill = MeanAbundance)) +
geom_tile(color = "white") +
scale_fill_gradientn(colors = viridis::viridis(10)) +
labs(title = "免疫细胞丰度热图", x = "", y = "") +
theme_minimal()
# 2. 显著差异细胞箱线图
signif_cells <- results_table %>% filter(Global.p < 0.05) %>% pull(CellType)
boxplot_list <- map(signif_cells, ~{
ggboxplot(cell_data, x = "Group", y = .x, fill = "Group",
palette = "jco", add = "jitter") +
stat_compare_means(method = ifelse(results_list[[.x]]$test %in% c("t-test", "ANOVA"),
"anova", "kruskal")) +
labs(title = paste(.x, "丰度差异"), y = "相对丰度")
})
# 3. 多组比较结果图（示例：CD4_Tcells）
if (nlevels(cell_data$Group) > 2) {
cd4_plot <- ggboxplot(cell_data, x = "Group", y = "CD4_Tcells",
color = "Group", palette = "npg") +
stat_pvalue_manual(
results_list$CD4_Tcells$result %>%
mutate(p.adj = format.pval(p.adj, digits = 2)),
label = "p.adj = {p.adj}",
y.position = max(cell_data$CD4_Tcells)*1.1
)
}
```
### 四、结果导出
```r
# 导出统计结果
write_csv(results_table, "immune_cell_stats_results.csv")
# 导出可视化图表
ggsave("immune_heatmap.png", heatmap_plot, width = 10, height = 8)
walk2(boxplot_list, signif_cells, ~ggsave(paste0(.y, "_boxplot.png"), .x, width = 8, height = 6))
```
### 五、生物学解释框架
对于显著结果的解读，建议结合以下生物学背景进行分析：
```r
# 1. T细胞亚群分析（参考[2]）
t_cell_analysis <- results_table %>%
filter(str_detect(CellType, "CD4|CD8|Treg")) %>%
arrange(Global.p)
# 2. B细胞与浆细胞关联（参考[2][3]）
b_cell_cor <- cell_data %>%
select(matches("Bcell|Plasma")) %>%
cor(method = "spearman")
# 3. 免疫抑制微环境指标（参考[3]）
immunosuppression_score <- cell_data %>%
mutate(ImmunoScore = Tregs + MDSCs - (CD8_Teff + NKcells))
```
### 注意事项
1. **数据要求**：确保输入数据为数值型矩阵，分组信息正确编码
2. **多重检验校正**：代码已内置Benjamini-Hochberg校正
3. **成分数据问题**：免疫细胞丰度总和恒定，建议添加：
```r
library(compositions)
clr_transformed <- cell_data %>%
select(-SampleID, -Group) %>%
clr() # 中心化对数比变换
```
4. **高级分析**：如需考虑混杂因素（年龄、性别等），推荐使用：
```r
library(lme4)
mixed_model <- lmer(CD4_Tcells ~ Group + Age + (1|Batch), data = cell_data)
```
> **关键生物学解读点**：
> - CD4+ T细胞亚群变化与自身免疫疾病进展相关[^2]
> - B细胞克隆扩增表征慢性炎症状态[^2]
> - Treg/MDSC升高提示免疫抑制微环境[^3]
> - NK细胞降低与免疫逃逸机制相关[^3]
本文介绍了如何利用R语言的'PMCMRplus'和'rcompanion'包执行Dunn's检验，这是一种非参数方法，用于检测不同组间的中位数差异。通过示例代码和结果解读，阐述了如何根据adj.p-value判断统计显著性，并强调了在解释结果时考虑多重比较和错误发现率的重要性。
订阅专栏 解锁全文
扫一扫
4733
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
