在R中计算特征部分匹配的两个矩阵的成对相关性方法

问题背景 #

分析高通量生物数据时，需要计算同一ProteinID对应的每个磷酸化位点与该ProteinID的成对相关系数，并为每个相关系数匹配对应的p值。现有两个列相同的数值矩阵：

• mat1：行名为磷酸化位点ID（格式如X1_Protein_1，下划线前为PhosphositeID，后为ProteinID，同一ProteinID对应多个PhosphositeID）
• mat2：行名为ProteinID，与mat1的ProteinID部分匹配

解决方案 #

基于已有的数据预处理代码，通过dplyr分组计算结合cor.test()函数，可同时得到分组的相关系数和对应p值：

library(tidyverse)

set.seed(123)
mat1 <- matrix(rnorm(100), nrow=10, ncol=10)
rownames(mat1) <- c("X1_Protein_1", "X2_Protein_1", "X3_Protein_2", "X4_Protein_3", 
"X5_Protein_3", "X6_Protein_4", "X7_Protein_5", "X8_Protein_6", "X9_Protein_7", "X10_Protein_8")
colnames(mat1) <- paste0("Sample", 1:10)

mat2 <- matrix(rnorm(70), nrow=7, ncol=10)
rownames(mat2) <- c("Protein_1", "Protein_2", "Protein_3", "Protein_4", "Protein_5", 
"Protein_6", "Protein_7")
colnames(mat2) <- paste0("Sample", 1:10)

# 提取mat1行名中的磷酸化位点和蛋白ID
mat1_df <- as.data.frame(mat1) %>%
  rownames_to_column(var="Phosphosite") %>%
  separate(Phosphosite, into=c("PhosphositeID", "ProteinID"), sep="_", extra = "merge")

# 筛选两者共有的ProteinID
common_proteins <- intersect(mat1_df$ProteinID, rownames(mat2))
mat1_subset <- mat1_df %>% filter(ProteinID %in% common_proteins)
mat2_subset <- mat2[rownames(mat2) %in% common_proteins, ]

# 转换为长格式并合并数据
mat1_df_final <- mat1_subset %>% 
  pivot_longer(cols = -c(ProteinID, PhosphositeID),
               names_to = "Sample", values_to = "phos_value")

mat2_df_final <- mat2_subset %>% 
  as.data.frame() %>% 
  rownames_to_column(var="ProteinID") %>% 
  pivot_longer(cols = -c(ProteinID),names_to = "Sample", values_to = "prot_value") 

to_cor <- mat1_df_final %>%
  left_join(mat2_df_final) %>%
  drop_na()

# 按磷酸化位点和蛋白ID分组计算相关系数及p值
cor_results <- to_cor %>%
  group_by(PhosphositeID, ProteinID) %>%
  summarize(
    pearson_correlation = cor(phos_value, prot_value, use = "pairwise.complete.obs"),
    p_value = cor.test(phos_value, prot_value)$p.value,
    .groups = "drop"  # 取消分组状态
  )

# 查看结果
print(cor_results)

代码说明 #

1. 分组逻辑：group_by(PhosphositeID, ProteinID)确保仅计算同一ProteinID下每个磷酸化位点与对应蛋白的相关性，避免全局无关组合的计算。
2. 统计量计算：
   • cor()默认计算皮尔逊相关系数，可通过method = "spearman"或"kendall"替换为非参数相关系数
   • cor.test()返回的结果中提取p.value，得到每个成对组合的显著性检验p值
3. 输出结果：最终的cor_results数据框包含四列，清晰对应每个磷酸化位点-蛋白组合的相关系数和p值，便于后续分析或可视化。

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