使用 R 基于数据帧中的因子的迭代减法答案

【问题标题】：Iterative subtraction based on factors in a data frame using R使用 R 基于数据帧中的因子的迭代减法
【发布时间】：2016-10-21 15:17:54
【问题描述】：

我正在努力想出一个可行的解决方案来解决这个看似相当简单的问题。我有一个包含数据和因子的数据框，我想使用这些因子来决定需要从其他数据点中减去哪些数据点以生成具有比较值的新数据框。

数据框是这样的：

str(means)
Classes ‘grouped_df’, ‘tbl_df’, ‘tbl’ and 'data.frame': 32 obs. of  5 variables:
 $ rat          : Factor w/ 8 levels "Rat1","Rat2",..: 1 1 1 1 2 2 2 2 3 3 ...
 $ gene         : Factor w/ 4 levels "gene1","gene2",..: 1 2 3 4 1 2 3 4 1 2 ...
 $ gene_category: Factor w/ 2 levels "control","experimental": 2 2 1 1 2 2 1 1 2 2 ...
 $ timepoint1   : num  23.4 18.3 42.1 40.1 25.3 ...
 $ timepoint2   : num  23.5 18.4 41.5 39.9 22.8 ...
> head(means)
Source: local data frame [6 x 5]
Groups: rat, gene [6]

 rat   gene gene_category timepoint1 timepoint2
(fctr) (fctr)        (fctr)      (dbl)      (dbl)
1   Rat1  gene1  experimental   23.36667   23.49667
2   Rat1  gene2  experimental   18.26000   18.38000
3   Rat1  gene3       control   42.05500   41.45000
4   Rat1  gene4       control   40.08667   39.89500
5   Rat2  gene1  experimental   25.29333   22.83000
6   Rat2  gene2  experimental   19.72667   19.19333

对于每只大鼠（总共 8 只大鼠），我想从“实验”基因值（基因 1 和 2）中减去“对照”基因值（基因 3 和 4）。我需要迭代地执行此操作，因此每个实验基因值必须从中减去每个对照基因值（在每只大鼠内，但不是在大鼠之间）。应该对每个时间点列进行上述操作。

我一直在摆弄使用 dplyr 的解决方案，我已经完成了分组，但我不知道如何做剩下的：

diffs <- means %>% group_by(rat, gene, gene_category) %>% here_is_where_i_don't_know_what_to_do)

There is a solution here to a similar problem here 但我认为它会给我所有可能的配对操作，这不是我想要的。它也只涉及两个因素，而我需要考虑三个因素。

Here's another solution to a similar problem，但同样有一些事情使它不太理想。它只处理一个因素，我不确定如何将其应用于具有三个因素和两个数据向量的数据集。

我知道在进行配对比较以确定统计显着性（多个 t 检验、ANOVA、MANOVA 等）时，此问题已得到解决，但我熟悉的用于执行这些测试的包/基本统计函数将此基本操作保留在引擎盖下。我想要一个简单的解决方案，使用 base R 或 dplyr/plyr/reshape2 等尽可能少的循环。

【问题讨论】：

您能否提供一些数据的输出？
我不太清楚你的计算。您是否需要将gene1 和gene2 的平均值以及gene3 和gene4 的平均值相减，或者它们是配对的（gene1 和gene3，...）？
@user124123 在上面添加了更清晰的数据框表示。
@bVa 从这个意义上说，它们不是配对的。基因 1 和 2 是实验基因，基因 3 和 4 是对照基因。我需要从实验基因值中减去控制基因值的实验基因和对照基因之间的所有可能比较。这些值已经是三次测量的平均值，所以我不想再次平均它们。这有意义吗？
因此您希望每只大鼠有 4 个结果：gene1 - gene3 |基因1 - 基因4 |基因2 - 基因3 |基因2 - 基因4？

标签： r dplyr

【解决方案1】：

我认为解决方案将涉及生成您想要的比较，然后将它们传递给标准评估 mutate_，而不是与 group_by 和 summarize 对抗。

首先，这里是读入的数据（注意，添加了 rat2 的基因 3/4）：

means <-
  read.table(text =
" rat   gene gene_category timepoint1 timepoint2
1   Rat1  gene1  experimental   23.36667   23.49667
2   Rat1  gene2  experimental   18.26000   18.38000
3   Rat1  gene3       control   42.05500   41.45000
4   Rat1  gene4       control   40.08667   39.89500
5   Rat2  gene1  experimental   25.29333   22.83000
6   Rat2  gene2  experimental   19.72667   19.19333
7   Rat2  gene3       control   42.05500   41.45000
8   Rat2  gene4       control   40.08667   39.89500")

接下来，在每个类中生成一组基因：

geneLists <-
  means %>%
  {split(.$gene, .$`gene_category`)} %>%
  lapply(unique) %>%
  lapply(as.character) %>%
  lapply(function(x){paste0("`", x, "`")})

请注意，反引号“`”是为了防止潜在的无效列名（例如，带有空格的东西）。这给出了：

$control
[1] "`gene3`" "`gene4`"

$experimental
[1] "`gene1`" "`gene2`"

然后，将您想要的比较粘贴在一起：

colsToCreate <-
  outer(geneLists[["experimental"]]
        , geneLists[["control"]]
        , paste, sep = " - ") %>%
  as.character()

给予：

[1] "`gene1` - `gene3`" "`gene2` - `gene3`" "`gene1` - `gene4`" "`gene2` - `gene4`"

然后，使用tidyr 传播数据，每只老鼠生成一行。注意，如果你想同时传播timepoint1 和timepoint2，你可能需要先gather（把两个时间放在一个列中），然后创建一个包含时间和基因的id 列，然后使用spread那个单一的 id 列。这还需要更改 colsToCreate 构造。

展开后，传递列的向量生成，你应该有你想要的：

means %>%
  select(rat, gene, timepoint1) %>%
  spread(gene, timepoint1) %>%
  mutate_(.dots = colsToCreate)

瞧：

   rat    gene1    gene2  gene3    gene4 gene1 - gene3 gene2 - gene3 gene1 - gene4 gene2 - gene4
1 Rat1 23.36667 18.26000 42.055 40.08667     -18.68833     -23.79500     -16.72000     -21.82667
2 Rat2 25.29333 19.72667 42.055 40.08667     -16.76167     -22.32833     -14.79334     -20.36000

实际上，获得两个时间点比我想象的要容易：

means %>%
  select(-gene_category) %>%
  gather("timepoint", "value", starts_with("timepoint")) %>%
  spread(gene, value) %>%
  mutate_(.dots = colsToCreate)

给予：

   rat  timepoint    gene1    gene2  gene3    gene4 gene1 - gene3 gene2 - gene3 gene1 - gene4 gene2 - gene4
1 Rat1 timepoint1 23.36667 18.26000 42.055 40.08667     -18.68833     -23.79500     -16.72000     -21.82667
2 Rat1 timepoint2 23.49667 18.38000 41.450 39.89500     -17.95333     -23.07000     -16.39833     -21.51500
3 Rat2 timepoint1 25.29333 19.72667 42.055 40.08667     -16.76167     -22.32833     -14.79334     -20.36000
4 Rat2 timepoint2 22.83000 19.19333 41.450 39.89500     -18.62000     -22.25667     -17.06500     -20.70167

另请注意，您可以命名包含列计算公式的向量，例如：

colsToCreate2 <-
  setNames(colsToCreate
           , c("nameA", "nameB", "nameC", "nameD"))

means %>%
  select(rat, gene, timepoint1) %>%
  spread(gene, timepoint1) %>%
  mutate_(.dots = colsToCreate2)

给予：

   rat    gene1    gene2  gene3    gene4     nameA     nameB     nameC     nameD
1 Rat1 23.36667 18.26000 42.055 40.08667 -18.68833 -23.79500 -16.72000 -21.82667
2 Rat2 25.29333 19.72667 42.055 40.08667 -16.76167 -22.32833 -14.79334 -20.36000

我不知道为什么，但是这个问题让我很兴奋，我想完成这个想法。在这里，我将gather 比较返回长格式，然后将时间点mutate 使用parse_number 从readr 和separate 将比较基因输出到单独的列中，以允许有效访问和过滤。请注意，每个基因的重复使用消除了独立性假设，因此不要在没有非常仔细地考虑控制的情况下对这些进行统计。

longForm <-
  means %>%
  select(-gene_category) %>%
  gather("timepoint", "value", starts_with("timepoint")) %>%
  spread(gene, value) %>%
  mutate_(.dots = colsToCreate) %>%
  select_(.dots = paste0("-",unlist(geneLists))) %>%
  gather(Comparison, Difference, -rat, -timepoint) %>%
  mutate(time = parse_number(timepoint)) %>%
  separate(Comparison, c("exp_Gene", "cont_Gene"), " - ")

head(longForm)

给予

   rat  timepoint exp_Gene cont_Gene Difference time
1 Rat1 timepoint1    gene1     gene3  -18.68833    1
2 Rat1 timepoint2    gene1     gene3  -17.95333    2
3 Rat2 timepoint1    gene1     gene3  -16.76167    1
4 Rat2 timepoint2    gene1     gene3  -18.62000    2
5 Rat1 timepoint1    gene2     gene3  -23.79500    1
6 Rat1 timepoint2    gene2     gene3  -23.07000    2

然后，我们可以绘制结果：

longForm %>%
  ggplot(aes(x = time
             , y = Difference
             , col = rat)) +
  geom_line() +
  facet_grid(exp_Gene ~ cont_Gene)

【讨论】：

您能否更详细地解释一下在构造 longForm 时“选择”的使用？有没有办法在不硬编码基因名称的情况下进行选择？那肯定会更好。
在构建longForm 时，所有select 正在删除gene_category 行，这允许spread 将所有基因放在同一行中（如果它离开gene_category，它将有一排包含所有 exp 基因，另一排包含所有对照基因，NA 填补空白）。在这种方法中，没有一个基因名称是硬编码的。它们都提取在geneLists 中。目前，colsToCreate 使用硬编码的类别名称，但您可以将其更改为仅使用 geneLists 的第一个和第二个名称（在最后使用单独命名列时需要这样做
感谢您提供详细信息，但实际上我指的是 select 的第二次使用，它确实按名称指定基因（这一行：select(-(gene1:gene4)) %>% ）。
我的错误：我错过了那个。 selection 在某种程度上是可选的。省略它只会在输出中留下单个基因值列，并使事情变得混乱。我现在修改了该调用以使用select_ 从geneLists 变量中删除所有基因。我还修改了geneLists 的创建，以在基因名称中添加反引号，以防将来的实现包含具有无效名称的基因（例如，带有空格）

【解决方案2】：

这是使用data.table 的latest devel version (1.9.7+) 的解决方案：

library(data.table)
setDT(means)

# join on rat being same and gene categories not being same, discard unmatched rows
# then extract interesting columns
means[means, on = .(rat, gene_category > gene_category), nomatch = 0,
      .(rat, gene.exp = gene, gene.ctrl = i.gene,
        timediff1 = timepoint1 - i.timepoint1, timediff2 = timepoint2 - i.timepoint2)]
#    rat gene.exp gene.ctrl timediff1 timediff2
#1: Rat1    gene1     gene3 -18.68833 -17.95333
#2: Rat1    gene2     gene3 -23.79500 -23.07000
#3: Rat1    gene1     gene4 -16.72000 -16.39833
#4: Rat1    gene2     gene4 -21.82667 -21.51500
#5: Rat2    gene1     gene3 -16.76167 -18.62000
#6: Rat2    gene2     gene3 -22.32833 -22.25667
#7: Rat2    gene1     gene4 -14.79334 -17.06500
#8: Rat2    gene2     gene4 -20.36000 -20.70167

如果你想推广到任意数量的“时间点”列：

nm = grep("timepoint", names(means), value = T)

means[means, on = .(rat, gene_category > gene_category), nomatch = 0,
      c(.(rat = rat, gene.exp = gene, gene.ctrl = i.gene),
        setDT(mget(nm)) - mget(paste0('i.', nm)))]

【讨论】：

这也是一个很好的解决方案，所以谢谢你。我接受了上面的“hadleyverse”答案，因为我正在教一些非编码人员在完成与他们的工作后使用此代码，我认为他们可能更容易从 dplyr 和其他包中解析代码，即使还有更多的事情要做。