首页 > 解决方案 > 过滤表以仅保留非冗余组

问题描述

我在 R 中进行系统发育分析已经有一段时间了,使用了 ape、phagorn 和 phytools 等库。

在解决问题时,我遇到了一个存在/不存在数据框架,它指定感兴趣的基因是否属于(或不属于)某个组。

这方面的一个例子是:

             gene11    gene25  gene33  gene54  gene55 gene65 gene73   gene88
group_1         1        1        0      0       0     0      0       0      
group_2         1        1        1      0       0     0      0       0      
group_3         1        0        1      0       0     0      0       0      
group_4         0        1        1      0       0     0      0       0      
group_5         0        0        0      1       1     0      0       0      
group_6         0        0        0      1       0     0      0       0      
group_7         0        0        0      0       1     0      0       0      
group_8         0        0        0      0       0     1      1       1      
group_9         0        0        0      0       0     1      1       0      
group_10        0        0        0      0       0     1      0       1      
group_11        0        0        0      0       0     0      1       1  

正如预期的那样,在处理生物实体组时,这些实体有多种关联方式:基因 11、25 和 33 形成一个组,它们的关系也可以描述为更小的组,描述成对关系。

所以这里很重要group_2group_5group_8是基因的生物学相关组,它们事先并不被称为相关组。其他较小的组是由于这些相关组中显示的关系而出现的:group_1 与gene11 和gene25 相关,但它是嵌套在更广泛(和相关)group_2 中的组。其他情况也是如此:group_8 描述了gene65、gene73 和gene88 之间的关系;与这些基因有关的其他组(group_9、group_10 和 group_11)只是描述作为更广泛组 group_8 的一部分的基因之间存在的成对关系的子组。

事先知道的是,基因形成了不相交的群,每个群由其他(逐渐变小的)群组成。我有兴趣捕捉最大的不相交的群体。

另一个用户(@Shree)对问题进行了明确定义:

找到最小数量的组,以使所有其他组是这些组中至少一个的子组。此外,一个组必须具有至少 2 个基因,即连续两个 1。还假设 1,01,0 是 的子组,1,1,1,00,1,1,1不是 的子组1,1,1,0

提前感谢大家!

标签: ralgorithmcluster-analysismathematical-optimizationhierarchical-clustering

解决方案


这是一种使用混合整数编程方法的方法。我正在使用ompr数学建模和glpk(免费开源)作为求解器。建模逻辑在代码中作为注释提供。

我认为这个问题可以在数学上描述如下 -

过滤数据框以最小化行数,使所有列的总和为 1。选定的行称为主要组,每隔一行应该是主要组的子组。一列(基因)只能属于一个主要组。subgroup <= primary group当处于所有位置(列)时,任何未选择的行都是主要组的子组。因此,(0,0,1,1)是 的子群,(0,1,1,1)(1,0,1,1)不是 的子群(0,1,1,1)

library(dplyr)
library(ROI)
library(ROI.plugin.glpk)
library(ompr)
library(ompr.roi)

gene_mat <- as.matrix(df)
nr <- nrow(gene_mat)
nc <- ncol(gene_mat)

model <- MIPModel() %>% 
  # binary variable x[i] is 1 if row i is selected else 0
  add_variable(x[i], i = 1:nr, type = "binary") %>% 
  # minimize total rows selected
  set_objective(sum_expr(x[i], i = 1:nr), "min") %>% 
  # sum of columns of selected rows must be = 1
  add_constraint(sum_expr(gene_mat[i,j]*x[i], i = 1:nr) == 1, j = 1:nc) %>% 
  solve_model(with_ROI(solver = "glpk"))

# get rows selected
group_rows <- model %>% 
  get_solution(x[i]) %>% 
  filter(value > 0) %>% 
  pull(i) %>% 
  print()

result <- df[group_rows, ]

        gene11 gene25 gene33 gene54 gene55 gene65 gene73 gene88
group_2      1      1      1      0      0      0      0      0
group_5      0      0      0      1      1      0      0      0
group_8      0      0      0      0      0      1      1      1

重要的提示 -

上述公式并未解决subgroup <= primary group,而是依赖于 OP 提到“事先已知的是基因形成不相交组的簇”这一事实。这意味着数据中不存在如下所示的情况,因为第 1、3、4 行不会形成不相交的组,即第 3 列将属于不允许的 2 个主要组。

1 1 0 0 0
0 1 0 0 0
1 0 1 0 0 <- this row is not a subgroup of any row
0 0 1 1 1

无论如何,这里的代码进行安全检查以确保所有未选择的行都是仅一个主要组的子组 -

test <- lapply(group_rows, function(x) {
  sweep(df, 2, as.numeric(df[x, ]), "<=") %>% 
    {which(rowSums(.) == ncol(df))}
})

# all is okay if below returns TRUE
length(Reduce(intersect, test)) == 0

数据 -

df <- structure(list(
  gene11 = c(1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L,0L, 0L), 
  gene25 = c(1L, 1L, 0L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), 
  gene33 = c(0L, 1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 0L), 
  gene54 = c(0L, 0L, 0L, 0L, 1L, 1L, 0L, 0L, 0L, 0L, 0L), 
  gene55 = c(0L, 0L, 0L, 0L, 1L, 0L, 1L, 0L, 0L, 0L, 0L), 
  gene65 = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 1L, 0L), 
  gene73 = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 1L, 0L, 1L), 
  gene88 = c(0L, 0L, 0L, 0L, 0L, 0L, 0L, 1L, 0L, 1L, 1L)), 
  class = "data.frame", 
  row.names = c("group_1", "group_2", "group_3", "group_4", 
                "group_5", "group_6", "group_7", "group_8",
                "group_9", "group_10", "group_11")
)

推荐阅读