筛选表以只保留非冗余组。

Question

我在R中做了一段时间的系统发育分析，使用了类人猿、穿山甲和系统工具等库。

在解决问题时，我遇到了一个存在/缺席的data.frame，它指定感兴趣的基因是否属于(或不属于)某个群体。

这方面的一个例子是：

             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_2、group_5和group_8是与生物相关的基因群，而事先并不知道是相关的组E 210。其他较小的组是由于这些相关组中所显示的关系而产生的: group_1与gene11和gene25相关，但它是嵌套在更广泛(且相关的) group_2中的组。在其他情况下也是如此: group_8描述了gene65、gene73和gene88之间的关系；与这些基因有关的其他组(group_9、group_10和group_11)只是描述属于更广泛的group_8的基因之间存在的成对关系的子组。

预先知道的是基因组成了不相交的群群，每个簇是由其他(逐渐小的)簇组成的。我有兴趣捕捉那些最大的不相交的团体。

问题的明确定义是由另一个用户(@Shree)完成的：

找到最少数量的组，以便所有其他组都是其中至少一个组的子组。另外，一个群体必须至少有2个基因，即连续两个1s。也假设1,01,0是1,1,1,0的子群，而0,1,1,1不是1,1,1,0的子群。

感谢所有提前！

Answer 1

这是一种使用混合整数规划方法的方法。我使用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")
)

Answer 2

这里有一个选项，我们用split创建一个分组索引(在这里，1和2对应于前3列和下2列)，然后循环遍历list，使用filter_all提取包含所有1s和mutate的行，将NA替换为0。

library(dplyr)
library(purrr)
library(tibble)
split.default(df, rep(1:2, c(3, 2))) %>%
     map_dfr(~ .x %>%
                  rownames_to_column('rn') %>% 
                  filter_at(-1, all_vars(.==1))) %>% 
        mutate_all(replace_na, 0) %>%
        column_to_rownames('rn')
#.       gene1 gene2 gene3 gene4 gene5
#group_1     1     1     1     0     0
#group_7     0     0     0     1     1

数据

df <- structure(list(gene1 = c(1L, 0L, 1L, 1L, 0L, 0L, 0L), gene2 = c(1L, 
1L, 1L, 0L, 0L, 0L, 0L), gene3 = c(1L, 1L, 0L, 1L, 0L, 0L, 0L
), gene4 = c(0L, 0L, 0L, 0L, 1L, 0L, 1L), gene5 = c(0L, 0L, 0L, 
0L, 0L, 1L, 1L)), .Names = c("gene1", "gene2", "gene3", "gene4", 
"gene5"), class = "data.frame", row.names = c("group_1", "group_2", 
"group_3", "group_4", "group_5", "group_6", "group_7"))