pbmc <- SetIdent(pbmc, value = pbmc$score_group) > DEGs <- FindMarkers(pbmc, + ident.1 = "High", + ident.2 = "Low", + test.use = "wilcox", # 可选方法："wilcox", "bimod", "roc", "t", "negbinom", "poisson", "LR", "MAST" + min.pct = 0.1, # 基因在至少一组中表达的细胞比例 + logfc.threshold = 0.25 # 对数倍变化阈值（log2 fold change） + )对pbmc进行细胞打分后，根据得分将所有细胞分为两组（High、Low），对两组进行FindMarkers差异分析，然后获得DEGs结果，下面进行结果可视化，要求做散点图

<think>我们正在处理单细胞RNA测序数据，已经使用FindMarkers函数找到了高分组和低分组之间的差异表达基因（DEGs）。现在需要对这些DEGs结果进行散点图可视化。 常见的散点图可视化方法有两种： 1. 平均表达量散点图：横轴和纵轴分别表示基因在两组（如Low组和High组）的平均表达量，可以直观展示基因在两组间的表达差异。 2. 火山图：横轴是log2FC（对数倍变化），纵轴是-log10(p-value)，用于同时展示差异倍数和统计显著性。 根据用户需求，这里要求绘制散点图，我们选择第一种：平均表达量散点图。 步骤： 1. 计算每个基因在两组（High和Low）的平均表达量。 2. 将FindMarkers得到的统计结果（如p_val, p_val_adj, avg_log2FC）与平均表达量合并。 3. 绘制散点图，并标记显著差异表达的基因。 注意：在单细胞数据中，平均表达量通常使用对数归一化后的值（如log1p(CP10K)），而Seurat的AverageExpression函数默认返回的是经过指数转换后的数据（即counts的指数转换，但实际是log-normalized值的反转换？）。我们需要明确： 在Seurat中，通常使用NormalizeData进行归一化（log1p(CP10K)），然后存储在RNA assay的data槽中。AverageExpression默认使用data槽（即归一化后的数据），但我们可以通过参数指定。 然而，AverageExpression返回的是按组平均的表达式，默认情况下，对于log-normalized数据，它返回的是未取对数的值（即先取平均，再指数减1？）。但为了与差异分析中的logFC一致，我们更希望使用对数空间的值。因此，我们可以设置`return.seurat = TRUE`来获得一个包含平均表达量的Seurat对象，然后提取数据，或者直接使用差异分析中使用的数据。 另一种做法：由于差异分析（FindMarkers）使用的是归一化后的数据（在log空间），我们可以直接计算每个组在log空间的平均表达量。 这里我们采用： avg_expr <- as.data.frame(AverageExpression(pbmc, group.by = "group", assays = "RNA", slot = "data")$RNA) 注意：slot="data"表示使用归一化后的数据（即log1p(CP10K)），这样每个基因的平均表达量也是在对数空间。 然后，我们将这个平均表达量数据框与FindMarkers的结果合并（通过基因名）。 绘制散点图时，横轴为Low组的平均表达量，纵轴为High组的平均表达量，并按照基因的显著性（如p_val_adj）和logFC来标记颜色。 另外，我们还可以标记一些重要的基因（如Top DEGs）。 代码步骤： 1. 计算平均表达量 2. 合并DEG结果 3. 绘制散点图 注意：由于平均表达量已经是log空间的值（但注意，log1p(CP10K)是自然对数），而差异分析中的logFC是自然对数还是2为底？在Seurat中，默认使用自然对数（即loge），但FindMarkers返回的avg_log2FC实际上是自然对数的FC（但命名为log2FC，这可能是Seurat的一个历史遗留问题，在Seurat v4中，FindMarkers返回的avg_log2FC是log2(avg_exp_group1/avg_exp_group2)）。因此，我们需要确保一致性。 实际上，AverageExpression(slot="data")返回的是log1p(CP10K)的平均值，而FindMarkers计算的logFC也是基于这个归一化后的数据，所以是匹配的。 具体代码： 假设我们已经有一个Seurat对象pbmc，其中有一个分组列"group"，包含"High"和"Low"两个分组。 首先，我们进行差异分析（如果还没做）： markers <- FindMarkers(pbmc, ident.1 = "High", ident.2 = "Low", group.by = "group") 然后，计算平均表达量： avg_expr <- as.data.frame(AverageExpression(pbmc, group.by = "group", assays = "RNA", slot = "data")$RNA) 现在，avg_expr是一个数据框，行名为基因，列名为组名（如High和Low）。 将markers结果（行名为基因）与avg_expr合并： avg_expr$gene <- rownames(avg_expr) markers$gene <- rownames(markers) plot_data <- merge(avg_expr, markers, by = "gene") 然后，我们绘制散点图： library(ggplot2) ggplot(plot_data, aes(x = Low, y = High)) + geom_point(aes(color = p_val_adj < 0.05 & abs(avg_log2FC) > log(1.5))) + # 例如logFC阈值设为log(1.5)（自然对数），但注意avg_log2FC实际上是log2，所以阈值应为log2(1.5) # 注意：FindMarkers返回的avg_log2FC是log2(平均表达量High / 平均表达量Low) # 因此，我们使用abs(avg_log2FC) > log2(1.5) 作为阈值 scale_color_manual(values = c("grey", "red")) + geom_abline(slope = 1, intercept = 0, linetype = "dashed") + # 由于数据是log空间，我们可以用对数坐标轴，但这里已经是log值，所以不需要再转换 # 但是，我们可能需要将坐标轴标签改为原始值？或者直接使用当前值。 # 另外，标记一些重要的基因 geom_text(data = subset(plot_data, p_val_adj < 0.01 & abs(avg_log2FC) > 1), aes(label = gene), size = 3, vjust = -0.5) + labs(x = "Low group (log1p(CP10K))", y = "High group (log1p(CP10K))", title = "Scatter plot of average expression") 但是，注意：AverageExpression(slot="data")返回的是log1p(CP10K)的平均值，所以每个点代表一个基因在组内细胞的平均log1p(CP10K)表达量。 另一种改进：使用scale_x_continuous(trans='log10')和scale_y_continuous(trans='log10')？但这里已经是log值（且可能包含0，取log10会出问题），所以不建议。 我们也可以使用scale_x_log10()和scale_y_log10()，但注意，我们的数据是log1p后的，实际上已经是连续非负值，但取log10会改变坐标。通常，我们直接使用线性坐标。 但是，由于我们取的是组内平均，所以表达量范围可能比较集中，线性坐标即可。 修正：由于数据是log1p(CP10K)，所以实际上已经是log转换后的值（自然对数），我们不需要再取对数坐标。因此，直接使用线性坐标。 注意：我们也可以将数据转换回原始尺度（通过expm1），但这样会使得点分布不均，且差异分析是基于log空间的，所以建议保留log空间。 因此，最终代码： ```r # 假设pbmc是Seurat对象，且已经进行了分组（group列包含High和Low） # 如果还没有进行差异分析，运行FindMarkers markers <- FindMarkers(pbmc, ident.1 = "High", ident.2 = "Low", group.by = "group") # 计算平均表达量（使用归一化数据，即log1p(CP10K)） avg_expr <- AverageExpression(pbmc, group.by = "group", assays = "RNA", slot = "data", return.seurat = FALSE) avg_expr <- as.data.frame(avg_expr$RNA) # 准备绘图数据 plot_data <- data.frame( gene = rownames(avg_expr), High = avg_expr[,"High"], Low = avg_expr[,"Low"], stringsAsFactors = FALSE ) # 合并差异分析结果 plot_data <- merge(plot_data, markers, by.x = "gene", by.y = "row.names") # 设置显著性标记 plot_data$significant <- "No" plot_data$significant[plot_data$p_val_adj < 0.05 & abs(plot_data$avg_log2FC) > 0.5] <- "Yes" # 绘制散点图 library(ggplot2) ggplot(plot_data, aes(x = Low, y = High)) + geom_point(aes(color = significant), alpha = 0.6, size = 1) + scale_color_manual(values = c("No" = "grey", "Yes" = "red")) + geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "blue") + # 标记Top10基因（按p值排序） ggrepel::geom_text_repel( data = subset(plot_data, significant == "Yes") %>% arrange(p_val_adj) %>% head(10), aes(label = gene), size = 3, max.overlaps = 20 ) + labs( x = "Low group (log1p(CP10K))", y = "High group (log1p(CP10K))", color = "Significant DEG", title = "DEGs between High and Low groups" ) + theme_bw() ``` 注意：这里使用了ggrepel包来避免标签重叠，如果没有安装需要先安装。 另外，阈值设置：avg_log2FC的绝对值大于0.5（即log2FC>0.5，相当于FC>1.4倍），并且调整后p值小于0.05。 如果数据点过多，我们可以使用点透明度（alpha）来避免重叠。 引用：在单细胞分析中，这种散点图是展示差异表达基因的常用方式[^1]。</think>### 单细胞DEGs散点图可视化方法 在单细胞RNA测序分析中，对High/Low组差异表达基因（DEGs）进行散点图可视化的核心步骤如下： #### 1. **准备数据** ```r # 计算两组平均表达量（使用对数归一化数据） avg_expr <- as.data.frame( Seurat::AverageExpression(pbmc, group.by = "group", assays = "RNA", slot = "data")$RNA ) # 合并差异分析结果 deg_results <- Seurat::FindMarkers(pbmc, ident.1 = "High", ident.2 = "Low") plot_data <- merge(avg_expr, deg_results, by = "row.names", suffixes = c("", "")) colnames(plot_data)[1] <- "gene" ``` #### 2. **创建基础散点图** ```r library(ggplot2) library(ggrepel) p <- ggplot(plot_data, aes(x = Low, y = High)) + geom_point(aes(color = p_val_adj < 0.05 & abs(avg_log2FC) > 1), alpha = 0.7, size = 2.5) + scale_color_manual(values = c("grey70", "red3"), name = "Significant DEG", labels = c("No", "Yes")) + geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "blue3") + scale_x_continuous(trans = "log10") + scale_y_continuous(trans = "log10") + labs(x = "Low组平均表达量 (log10)", y = "High组平均表达量 (log10)", title = "High vs Low组差异表达基因") + theme_classic(base_size = 14) + theme(legend.position = "bottom") ``` #### 3. **添加关键基因标签** ```r # 筛选Top10显著基因 top_genes <- plot_data %>% filter(p_val_adj < 0.01 & abs(avg_log2FC) > 1.5) %>% arrange(p_val_adj) %>% head(10) # 添加标签 p + ggrepel::geom_text_repel( data = top_genes, aes(label = gene), size = 4, box.padding = 0.5, max.overlaps = 20, segment.color = "grey50" ) ``` #### 4. **增强可视化效果** ```r p + stat_density_2d(aes(fill = ..level..), geom = "polygon", alpha = 0.2) + scale_fill_gradient(low = "lightblue", high = "darkblue") + annotate("text", x = max(plot_data$Low)*0.1, y = max(plot_data$High)*0.9, label = paste("DEGs:", sum(plot_data$p_val_adj < 0.05)), color = "red3", size = 5) ``` ### 关键参数说明 | 参数 | 推荐值 | 作用 | |------|--------|------| | `slot` | "data" | 使用对数归一化数据 | | `alpha` | 0.5-0.8 | 点透明度避免重叠 | | `logFC阈值` | >1 | 显著差异倍数 | | `p_val_adj` | <0.05 | 调整后p值阈值 | | `max.overlaps` | 20-30 | 标签最大重叠数 | ### 生物学解读技巧 1. **对角线分布**：靠近对角线的基因表示组间无差异，远离对角线的点表示显著DEGs 2. **象限分析**： - 右上象限：高表达基因（如管家基因） - 左上象限：High组特异性表达基因 - 右下象限：Low组特异性表达基因 3. **关注关键基因**：如T细胞相关基因CD3D、CD8A在High组的表达模式[^1] > **注意事项**：当数据存在批次效应时，需先用`harmony`或`Seurat::IntegrateData()`进行整合处理[^1]。 ---