如何用group_by和summarise手动计算log2倍数变化

以下为用户提供的数据集与样本信息：

example_mat <-structure(c(229.705708222381, 1434.77719289672, 303.918321648074, 
235.006609181359, 365.762166169484, 263.278080962576, 394.552456179431, 
948.547343280687, 279.699857634049, 299.967963259704, 553.994887101256, 
431.035046305611, 357.731674441394, 1420.98970680887, 395.672912639724, 
359.538400069886, 437.227602095038, 386.639284497265, 149.584664657699, 
1883.65874013399, 82.1792293489829, 645.430127134147, 1108.03455301999, 
222.530272731516, 128.261702714215, 1440.4126745603, 138.387626612706, 
512.20298386532, 672.530112258089, 151.88885847736, 168.886143824754, 
1484.04207233241, 73.6631052852648, 648.594658731234, 774.36093604754, 
140.139566152455, 144.335247151679, 934.169794065035, 85.7992858068316, 
534.842276945389, 376.073505078542, 221.314045632575, 268.50050365722, 
219.769095520145, 95.5517806609325, 513.113062149208, 414.694728068447, 
280.922235143142, 118.580483523586, 810.061857323772, 281.450063302968, 
664.336443836956, 350.026728473235, 188.585829218233, 206.981437576129, 
1080.15429283915, 150.181880287796, 622.869721450026, 329.244891400168, 
280.146968998389, 130.819453108194, 942.081756063866, 297.977643190885, 
415.170069933642, 698.612218334728, 192.595305964841, 166.172073025442, 
264.552554070355, 456.353155771362, 398.482284568472, 753.974779100512, 
235.61711846891, 96.1731056297284, 1136.19893381788, 382.842939718342, 
491.345930685215, 843.980651327552, 200.360636728601, 53.1514867212154, 
656.04120753043, 259.682977980795, 410.025754706519, 923.317255042828, 
148.824162819403, 42.8797353525399, 163.541316228292, 431.788962968599, 
257.278412115239, 644.193233435832, 141.602846978155, 48.5032485391316, 
82.6351641777798, 627.847606089871, 195.809410769087, 795.812559364271, 
200.300452300488, 86.6629185441187, 461.553071796767, 185.984465639625, 
454.736887192173, 687.460904406155, 126.586285513881, 157.662405351439, 
100.635577883897, 155.146515904342, 410.089979876881, 635.681400299951, 
144.244328300253, 104.411017961498, 446.447111283645, 189.019946309608, 
531.056039631756, 459.048441037619, 214.222605817556, 113.574061770839, 
209.955719695634, 210.476701630362, 443.355626454055, 496.495783796374, 
140.665122376727), dim = c(6L, 20L), dimnames = list(c("STX2", 
"IKZF3", "CTLA4", "CHAF1A", "ST3GAL5", "PCIF1"), c("AMLA1_0d", 
"AMLA2_0d", "AMLA3_0d", "AMLB1_0d", "AMLB2_0d", "AMLB3_0d", "AML10b_0d", 
"AML11b_0d", "AML8_0d", "AML9b_0d", "AMLA1_14d", "AMLA2_14d", 
"AMLA3_14d", "AMLB1_14d", "AMLB2_14d", "AMLB3_14d", "AML10b_14d", 
"AML11b_14d", "AML8_14d", "AML9b_14d")))

pd_des_example <- as.factor(colnames(example_mat)) %>%
  as.data.frame() %>%
  setNames(c("SampleName")) %>%
  mutate(Type = ifelse(grepl("AMLA", SampleName), "AMLA",
                       ifelse(grepl("AMLB", SampleName), "AMLB", "AML"))) %>%
  mutate(Time = ifelse(grepl("_0d", SampleName), "d1", "d14"),
         Group =  sub("^([^_]*_[^_]*).*", "\\1", SampleName))

目标计算逻辑 #

需基于以下公式计算log2倍数变化（log2Fold Change）：

basal_mean <- mean(log2(row[subseted_pd$Time == "d1"] + 1))
stimulated_mean <- mean(log2(row[subseted_pd$Time == "d14"]+ 1))
log2FC <- basal_mean - stimulated_mean

尝试的代码（存在问题） #

g_first_mat <- example_mat %>%  as.data.frame %>%
  tibble::rownames_to_column(var = "gene")%>%
  gather(SampleName, measurement, -gene) %>%
  left_join(pd.des) %>%
  spread(Time, measurement) %>%
  dplyr::group_by(Type) %>%
  summarize(d1_geometric_mean = geometric.mean(d1, na.rm=TRUE),
            d14_geometric_mean = geometric.mean(d14, na.rme=TRUE)) %>%
  mutate(geometric_mean_ratio = d1_geometric_mean / d14_geometric_mean)

问题分析与修正 #

核心错误点 #

1. 变量/参数拼写错误：left_join(pd.des) 应为 left_join(pd_des_example)；na.rme=TRUE 应为 na.rm=TRUE。
2. 函数依赖缺失：geometric.mean 不是Base R或dplyr内置函数，需使用psych::geometric.mean或自行定义。
3. 分组逻辑偏差：原目标需按基因+Type分组计算，而非仅按Type聚合。

修正后的代码 #

首先加载必要依赖包：

library(dplyr)
library(tidyr)
library(psych) # 用于几何均值计算

方案1：按原log2算术均值逻辑计算

log2fc_result <- example_mat %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  gather(SampleName, measurement, -gene) %>%
  left_join(pd_des_example, by = "SampleName") %>%
  group_by(gene, Type) %>%
  summarize(
    basal_mean = mean(log2(measurement[Time == "d1"] + 1), na.rm = TRUE),
    stimulated_mean = mean(log2(measurement[Time == "d14"] + 1), na.rm = TRUE),
    log2FC = basal_mean - stimulated_mean
  ) %>%
  ungroup()

方案2：使用几何均值计算

geo_mean_result <- example_mat %>%
  as.data.frame() %>%
  tibble::rownames_to_column("gene") %>%
  gather(SampleName, measurement, -gene) %>%
  left_join(pd_des_example, by = "SampleName") %>%
  group_by(gene, Type) %>%
  summarize(
    d1_geo_mean = psych::geometric.mean(measurement[Time == "d1"] + 1, na.rm = TRUE),
    d14_geo_mean = psych::geometric.mean(measurement[Time == "d14"] + 1, na.rm = TRUE),
    log2FC_geo = log2(d1_geo_mean / d14_geo_mean)
  ) %>%
  ungroup()

额外说明 #

若不想依赖psych包，可自行定义几何均值函数：

geometric_mean <- function(x, na.rm = FALSE) {
  if (na.rm) x <- x[!is.na(x)]
  exp(mean(log(x[x > 0])))
}

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