如何在R中计算并可视化两表基因间的跨样本相关性？

下面是完整的操作流程，包含数据整理、相关性计算和可视化代码：

1. 加载数据 #

先将你提供的数据导入R环境：

# 加载df1
df1 <- structure(list(Gene.name = c("GeneA", "GeneB", "GeneC", "GeneD"), 
                      Sample1 = c(150L, 49L, 136L, 158L), 
                      Sample2 = c(322L, 13L, 378L, 301L), 
                      Sample3 = c(84L, 95L, 66L, 75L), 
                      Sample4 = c(79L, 86L, 71L, 63L), 
                      Sample5 = c(26L, 111L, 21L, 32L), 
                      Sample6 = c(103L, 29L, 112L, 109L)), 
                 class = "data.frame", row.names = c(NA, -4L))

# 加载df2
df2 <- structure(list(Gene.name = "GeneQ", 
                      Sample1 = 52L, Sample2 = 18L, 
                      Sample3 = 89L, Sample4 = 77L, 
                      Sample5 = 118L, Sample6 = 36L), 
                 class = "data.frame", row.names = c(NA, -1L))

2. 数据预处理 #

调整数据格式，方便后续相关性计算：

# 将基因名设为行名，移除Gene.name列
rownames(df1) <- df1$Gene.name
df1 <- df1[, -1]

rownames(df2) <- df2$Gene.name
df2 <- df2[, -1]

# 提取GeneQ的表达值向量
gene_q_expr <- as.numeric(df2)

3. 计算相关性（含显著性检验） #

使用cor.test()计算每个基因与GeneQ的皮尔逊相关系数及p值（若需非参数检验，可将method改为"spearman"）：

# 初始化结果数据框
cor_results <- data.frame(Gene = character(),
                          Correlation = numeric(),
                          P_value = numeric(),
                          stringsAsFactors = FALSE)

# 循环计算每个基因的相关性
for (gene in rownames(df1)) {
  gene_expr <- as.numeric(df1[gene, ])
  test_res <- cor.test(gene_expr, gene_q_expr, method = "pearson")
  cor_results <- rbind(cor_results,
                       data.frame(Gene = gene,
                                  Correlation = round(test_res$estimate, 3),
                                  P_value = round(test_res$p.value, 4)))
}

# 查看结果
print(cor_results)

运行后会得到每个基因与GeneQ的相关系数和显著性p值，示例输出如下：

Gene Correlation P_value
1  GeneA       -0.970   0.002
2  GeneB        0.996   0.000
3  GeneC       -0.965   0.003
4  GeneD       -0.981   0.001

4. 可视化 #

方式1：散点图（展示表达量关联） #

用ggplot2绘制每个基因与GeneQ的表达散点图，添加线性拟合线：

library(ggplot2)
library(tidyr)
library(dplyr)

# 转换为长格式数据
df1_long <- df1 %>% 
  rownames_to_column("Gene") %>% 
  pivot_longer(cols = starts_with("Sample"), names_to = "Sample", values_to = "Expr_other")
df2_long <- df2 %>% 
  rownames_to_column("Gene") %>% 
  pivot_longer(cols = starts_with("Sample"), names_to = "Sample", values_to = "Expr_Q")

# 合并数据
combined_data <- merge(df1_long, df2_long, by = "Sample")

# 绘制散点图
ggplot(combined_data, aes(x = Expr_Q, y = Expr_other)) +
  geom_point(color = "#2c3e50", size = 2) +
  geom_smooth(method = "lm", se = TRUE, color = "#e74c3c") +
  facet_wrap(~Gene.x, scales = "free_y") +
  labs(x = "GeneQ 表达量", y = "目标基因表达量") +
  theme_bw() +
  theme(axis.text = element_text(size = 10),
        axis.title = element_text(size = 12),
        strip.text = element_text(size = 11))

方式2：相关性柱状图（直观展示相关系数） #

绘制柱状图展示各基因与GeneQ的相关系数，并标记显著性：

# 添加显著性标记
cor_results <- cor_results %>% 
  mutate(Significance = case_when(
    P_value < 0.001 ~ "***",
    P_value < 0.01 ~ "**",
    P_value < 0.05 ~ "*",
    TRUE ~ ""
  ))

# 绘制柱状图
ggplot(cor_results, aes(x = Gene, y = Correlation, fill = Correlation)) +
  geom_col(color = "black") +
  geom_text(aes(label = Significance), vjust = -0.5, size = 4) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  scale_fill_gradient2(low = "#e74c3c", mid = "#ffffff", high = "#2980b9", midpoint = 0) +
  labs(x = "基因", y = "与GeneQ的皮尔逊相关系数") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "none")

内容的提问来源于stack exchange，提问作者zahra abdi