我有一个由标签分隔的细菌菌株名称和基因的二进制矩阵,列为存在(1(或不存在(0(,由ROARY(泛基因组管道(输出。
这是数据的模拟版本:
strain <- rep(letters[1:4], 5)
gene <- c(rep("G1", 4), rep("G2", 4), rep("G3", 4), rep("G4", 4), rep("G5", 4))
pres_abs <- c(1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0, 0, 1, 1, 0, 1, 1)
test <- tibble(strain, gene, pres_abs)
strain gene pres_abs
<chr> <chr> <dbl>
1 a G1 1
2 b G1 1
3 c G1 1
4 d G1 1
5 a G2 1
6 b G2 1
7 c G2 0
8 d G2 0
9 a G3 0
10 b G3 0
11 c G3 1
12 d G3 1
13 a G4 0
14 b G4 0
15 c G4 0
16 d G4 1
17 a G5 1
18 b G5 0
19 c G5 1
20 d G5 1
顺便说一句,当我使用read_tsv()
:将它读入R时,它的结构是这样的
gene a b c d
<chr> <dbl> <dbl> <dbl> <dbl>
1 G1 1 1 1 1
2 G2 1 1 0 0
3 G3 0 0 1 1
4 G4 0 0 0 1
5 G5 1 0 1 1
我的基质中有几千个基因和大约30个菌株。
我想识别菌株的某个子集中缺失的所有基因(0(,并将其保存为载体(列表?(,用于进一步分析(例如,作为类似数据帧的过滤项(。
对于上面的例子,我只想要菌株a和菌株b中都不存在的基因(因此存在于c和/或d中(。所以我希望得到基因G3和G4。
在搜索了一些解决方案后,我使用pivot_longer
扩展了数据,使其结构类似于我的示例中的test
。我试着这样过滤:
test %>% filter(strain %in% c("a", "b") & pres_abs == 0)
这得到了我想要的G3和G4,但也得到了G5,我没有,因为它存在于基因a中。
strain gene pres_abs
<chr> <chr> <dbl>
1 a G3 0
2 b G3 0
3 a G4 0
4 b G4 0
5 b G5 0
有人能帮我输入正确的过滤条件吗?
以下是我们对您感兴趣的菌株的pres_abs
求和的方法。
test %>%
group_by(gene) %>%
filter(sum(pres_abs[strain %in% c("a", "b")]) == 0)
# # A tibble: 8 x 3
# # Groups: gene [2]
# strain gene pres_abs
# <chr> <chr> <dbl>
# 1 a G3 0
# 2 b G3 0
# 3 c G3 1
# 4 d G3 1
# 5 a G4 0
# 6 b G4 0
# 7 c G4 0
# 8 d G4 1
以上返回了这些菌株的所有观察结果。或者,您可以进行两步过滤:
test %>%
group_by(gene) %>%
filter(strain %in% c("a", "b")) %>%
filter(sum(pres_abs) == 0)
# # A tibble: 4 x 3
# # Groups: gene [2]
# strain gene pres_abs
# <chr> <chr> <dbl>
# 1 a G3 0
# 2 b G3 0
# 3 a G4 0
# 4 b G4 0
我们可以通过操作进行分组,即按"基因"分组,检查all
"a"、"b"是否在"菌株"中找到,其中"pres_abs"值为0,为了避免在pres_abs中得到1值,创建第二个条件,即"pres_abs'"为0
library(dplyr)
test %>%
group_by(gene) %>%
filter(all(c("a", "b") %in% strain[pres_abs == 0]),
pres_abs == 0) %>%
ungroup
-ouptut
# A tibble: 5 x 3
strain gene pres_abs
<chr> <chr> <dbl>
1 a G3 0
2 b G3 0
3 a G4 0
4 b G4 0
5 c G4 0
如果我们也需要1值,
test %>%
group_by(gene) %>%
filter(all(c("a", "b") %in% strain[pres_abs == 0])) %>%
ungroup
-输出
# A tibble: 8 x 3
strain gene pres_abs
<chr> <chr> <dbl>
1 a G3 0
2 b G3 0
3 c G3 1
4 d G3 1
5 a G4 0
6 b G4 0
7 c G4 0
8 d G4 1
不要为使数据变长而烦恼。我把它改回宽,然后从那里过滤。
test %>%
pivot_wider(gene, names_from = strain, values_from = pres_abs, values_fill = 0) %>%
filter(across(c(a,b), ~ .==0)) %>%
pull(gene)
[1] "G3" "G4"
编辑:显示带有跨工作示例的any_of
test %>%
pivot_wider(gene, names_from = strain, values_from = pres_abs, values_fill = 0) %>%
filter(across(any_of(c("a","b","z","q","t")), ~ .==0)) %>%
pull(gene)
[1] "G3" "G4"