最近我得到了一个数据框，里面有几列，然后我要在这几列里面根据一些过滤条件得到基因，然后进行GO分析。比如下面这个数据框
# gene是我随机抽取的基因
library(tidyverse)
library(clusterProfiler)
library(org.At.tair.db)
data.frame(gene = sample(gene,10),
rho = sample(seq(-1,1,length.out = 1000),10),
log2Fold = sample(seq(-10,10,length.out = 100),10),
stringsAsFactors = F) -> test_data
> test_data
gene rho log2Fold
1 AT5G01640 -0.96996997 3.3333333
2 CH5 -0.67367367 1.5151515
3 AT5G07470 -0.99399399 -10.0000000
4 AT5G19140 -0.29129129 -0.1010101
5 EMB79 0.02902903 -1.5151515
6 AT5G54070 -0.65565566 -4.1414141
7 AT5G46315 -0.77177177 -5.7575758
8 AT3G55040 0.01701702 -8.1818182
9 AT1G56410 0.43743744 3.7373737
10 AT3G05360 -0.74774775 -9.7979798
比如我想得到分别得到rho > 0和 < 0的基因，那么我就会用
# 取出来基因
test_data %>%
filter(rho > 0) %>%
pull(gene) %>%
unique() -> gene_rho_1
test_data %>%
filter(rho < 0) %>%
pull(gene) %>%
unique() -> gene_rho_2
# 做GO
ego_1 <- enrichGO(gene = gene_rho_1,
OrgDb = org.At.tair.db,
keyType = "TAIR",
ont = "BP")
ego_2 <- enrichGO(gene = gene_rho_2,
OrgDb = org.At.tair.db,
keyType = "TAIR",
ont = "BP")
但实际上Y叔的clusterProfiler包里面是有一个叫做compareCluster的函数，有了这个函数，我们就可以对 rho >0 和 rho <0一下子都做GO了。
# compareCluster要求传入的是一个列表
# 所以我们这里做一个列表
gene_rho[["p"]] <- test_data %>%
filter(rho > 0) %>%
pull(gene) %>%
unique()
gene_rho[["n"]] <- test_data %>%
filter(rho < 0) %>%
pull(gene) %>%
unique()
# 然后传入compareCluster
compareCluster(geneClusters = gene_rho,
fun = "enrichGO",
OrgDb = org.At.tair.db,
keyType = "TAIR",
ont = "BP") -> ego
但实际上，我不仅是想对rho进行筛选，我还想对log2Fold进行筛选。那我就需要4次filter，才能构造出能够传入compareCluster的列表了。
# 这样做四次筛选
gene_rho[["p_p"]] <- test_data %>%
filter(rho > 0 & log2Fold > 0) %>%
pull(gene) %>%
unique()
本来想构造一个函数，但这种筛选条件比较复杂，不知道怎么给filter传入参数（大家知道的也可以在下面留言）。然后我就突然想到了一个方法
# 加两列，作为我们后续group的条件
test_data %>%
mutate(r = rho > 0,
l = log2Fold > 0)
gene rho log2Fold r l
1 AT5G01640 -0.96996997 3.3333333 FALSE TRUE
2 CH5 -0.67367367 1.5151515 FALSE TRUE
3 AT5G07470 -0.99399399 -10.0000000 FALSE FALSE
4 AT5G19140 -0.29129129 -0.1010101 FALSE FALSE
5 EMB79 0.02902903 -1.5151515 TRUE FALSE
6 AT5G54070 -0.65565566 -4.1414141 FALSE FALSE
7 AT5G46315 -0.77177177 -5.7575758 FALSE FALSE
8 AT3G55040 0.01701702 -8.1818182 TRUE FALSE
9 AT1G56410 0.43743744 3.7373737 TRUE TRUE
10 AT3G05360 -0.74774775 -9.7979798 FALSE FALSE
# 这样我们的数据就顺利地分成了4组
test_data %>%
mutate(r = rho > 0,
l = log2Fold > 0) %>%
group_by(r,l)
# A tibble: 10 x 5
# Groups: r, l [4]
gene rho log2Fold r l
<fct> <dbl> <dbl> <lgl> <lgl>
1 AT5G01640 -0.970 3.33 FALSE TRUE
2 CH5 -0.674 1.52 FALSE TRUE
3 AT5G07470 -0.994 -10 FALSE FALSE
4 AT5G19140 -0.291 -0.101 FALSE FALSE
5 EMB79 0.0290 -1.52 TRUE FALSE
6 AT5G54070 -0.656 -4.14 FALSE FALSE
7 AT5G46315 -0.772 -5.76 FALSE FALSE
8 AT3G55040 0.0170 -8.18 TRUE FALSE
9 AT1G56410 0.437 3.74 TRUE TRUE
10 AT3G05360 -0.748 -9.80 FALSE FALSE
在分成了4组之后，考虑到compareCluster需要传入的是list，所以我们就需要把group结果直接切割成含有4组的list。这时候就要用到group_split了
test_data %>%
mutate(r = rho > 0,
l = log2Fold > 0) %>%
group_split(r,l)
[[1]]
# A tibble: 5 x 5
gene rho log2Fold r l
<fct> <dbl> <dbl> <lgl> <lgl>
1 AT5G07470 -0.994 -10 FALSE FALSE
2 AT5G19140 -0.291 -0.101 FALSE FALSE
3 AT5G54070 -0.656 -4.14 FALSE FALSE
4 AT5G46315 -0.772 -5.76 FALSE FALSE
5 AT3G05360 -0.748 -9.80 FALSE FALSE
[[2]]
# A tibble: 2 x 5
gene rho log2Fold r l
<fct> <dbl> <dbl> <lgl> <lgl>
1 AT5G01640 -0.970 3.33 FALSE TRUE
2 CH5 -0.674 1.52 FALSE TRUE
[[3]]
# A tibble: 2 x 5
gene rho log2Fold r l
<fct> <dbl> <dbl> <lgl> <lgl>
1 EMB79 0.0290 -1.52 TRUE FALSE
2 AT3G55040 0.0170 -8.18 TRUE FALSE
[[4]]
# A tibble: 1 x 5
gene rho log2Fold r l
<fct> <dbl> <dbl> <lgl> <lgl>
1 AT1G56410 0.437 3.74 TRUE TRUE
每个列表的名字可以由group_keys来查看，然后你手动names一下，更改下列表元素的名字就行了。
# 可以看到分别是
> test_data %>%
+ mutate(r = rho > 0,
+ l = log2Fold > 0) %>%
+ group_keys(r,l)
# A tibble: 4 x 2
r l
<lgl> <lgl>
1 FALSE FALSE
2 FALSE TRUE
3 TRUE FALSE
4 TRUE TRUE
但还有个问题就是我们得到了列表内元素还是tibble，并不是我们想要的基因，所以我们还要进行一步操作lappy的操作
# 这样我们得到的就是一个含有4组元素，且元素里面都是基因的列表了
# 然后就可以传入compareCluster批量做GO了。
test_data %>%
mutate(r = rho > 0,
l = log2Fold > 0) %>%
group_split(r,l) %>%
lapply(.,function(x) unique(pull(x,gene)))
