专栏名称: 生信媛
生信媛，从1人分享，到8人同行。坚持分享生信入门方法与课程，持续记录生信相关的分析pipeline, python和R在生物信息学中的利用。内容涵盖服务器使用、基因组转录组分析以及群体遗传。
目录
相关文章推荐
51好读  ›  专栏  ›  生信媛

数据框筛选和GO富集的小操作

生信媛  · 公众号  · 生物  · 2020-01-07 07:45

正文

请到「今天看啥」查看全文


最 近我得到了一个数据框，里面有几列，然后我要在这几列里面根据一些过滤条件得到基因，然后进行GO分析。 比如下面这个数据框 

  1. # gene是我随机抽取的基因
  2. library(tidyverse)
  3. library(clusterProfiler)
  4. library(org.At.tair.db)

  6. data.frame(gene = sample(gene,10),
  7.            rho = sample(seq(-1,1,length.out = 1000),10),
  8.            log2Fold = sample(seq(-10,10,length.out = 100),10),
  9.            stringsAsFactors = F) -> test_data

  11. > test_data




    

  12.         gene         rho    log2Fold
  13. 1  AT5G01640 -0.96996997   3.3333333
  14. 2        CH5 -0.67367367   1.5151515
  15. 3  AT5G07470 -0.99399399 -10.0000000
  16. 4  AT5G19140 -0.29129129  -0.1010101
  17. 5      EMB79  0.02902903  -1.5151515
  18. 6  AT5G54070 -0.65565566  -4.1414141
  19. 7  AT5G46315 -0.77177177  -5.7575758
  20. 8  AT3G55040  0.01701702  -8.1818182
  21. 9  AT1G56410  0.43743744   3.7373737
  22. 10 AT3G05360 -0.74774775  -9.7979798

比如我想得到分别得到rho > 0和 < 0的基因，那么我就会用 

  1. # 取出来基因
  2. test_data %>%




    
 
  3.   filter(rho > 0) %>% 
  4.   pull(gene) %>% 
  5.   unique() -> gene_rho_1

  7. test_data %>% 
  8.   filter(rho < 0) %>% 
  9.   pull(gene) %>% 
  10.   unique() -> gene_rho_2

  12. # 做GO
  13. ego_1  enrichGO(gene = gene_rho_1,
  14.                 OrgDb = org.At.tair.db,
  15.                 keyType = "TAIR",
  16.                 ont = 




    
"BP")

  18. ego_2  enrichGO(gene = gene_rho_2,
  19.                 OrgDb = org.At.tair.db,
  20.                 keyType = "TAIR",
  21.                 ont = "BP")

但实际上Y叔的clusterProfiler包里面是有一个叫做 compareCluster 的函数，有了这个函数，我们就可以对 rho >0 和 rho <0一下子都做GO了。 

  1. # compareCluster要求传入的是一个列表
  2. # 所以我们这里做一个列表

  4. gene_rho[["p"]]  test_data %>% 
  5.   filter(rho > 0) %>% 
  6.   pull(gene) %>% 
  7.   unique()
  8. gene_rho[["n"]]  test_data %>% 




    

  9.   filter(rho < 0) %>% 
  10.   pull(gene) %>% 
  11.   unique()

  13. # 然后传入compareCluster
  14. compareCluster(geneClusters = gene_rho,
  15.                       fun = "enrichGO",
  16.                       OrgDb = org.At.tair.db,
  17.                       keyType = "TAIR",
  18.                       ont = "BP") -> ego

但实际上，我不仅是想对rho进行筛选，我还想对log2Fold进行筛选。那我就需要4次filter，才能构造出能够传入compareCluster的列表了。 

  1. # 这样做四次筛选
  2. gene_rho[["p_p"]]  test_data %>% 
  3.   filter(rho > 0




    
 & log2Fold > 0) %>% 
  4.   pull(gene) %>% 
  5.   unique()

本来想构造一个函数，但这种筛选条件比较复杂，不知道怎么给filter传入参数（大家知道的也可以在下面留言）。然后我就突然想到了一个方法 

  1. # 加两列，作为我们后续group的条件
  2. test_data %>% 
  3.   mutate(r = rho > 0,
  4.          l = log2Fold > 0)

  6.         gene         rho    log2Fold     r     l
  7. 1  AT5G01640 -0.96996997   3.3333333 FALSE  TRUE
  8. 2        CH5 -0.67367367   1.5151515 FALSE  TRUE
  9. 3  AT5G07470 -0.99399399 -10.0000000 FALSE FALSE
  10. 4  AT5G19140 -0.29129129  -0.1010101 FALSE FALSE
  11. 5




    
      EMB79  0.02902903  -1.5151515  TRUE FALSE
  12. 6  AT5G54070 -0.65565566  -4.1414141 FALSE FALSE
  13. 7  AT5G46315 -0.77177177  -5.7575758 FALSE FALSE
  14. 8  AT3G55040  0.01701702  -8.1818182  TRUE FALSE
  15. 9  AT1G56410  0.43743744   3.7373737  TRUE  TRUE
  16. 10 AT3G05360 -0.74774775  -9.7979798 FALSE FALSE


  19. # 这样我们的数据就顺利地分成了4组
  20. test_data %>% 
  21.   mutate(r = rho > 0,
  22.          l = log2Fold > 0) %>% 
  23.   group_by(r,l)

  25. # A tibble: 10 x 5
  26. # Groups:   r, l [4]
  27.    gene          rho log2Fold r     l




    

  28.                 
  29.  1 AT5G01640 -0.970     3.33  FALSE TRUE
  30.  2 CH5       -0.674     1.52  FALSE TRUE
  31.  3 AT5G07470 -0.994   -10     FALSE FALSE
  32.  4 AT5G19140 -0.291    -0.101 FALSE FALSE
  33.  5 EMB79      0.0290   -1.52  TRUE  FALSE
  34.  6 AT5G54070 -0.656    -4.14  FALSE FALSE
  35.  7 AT5G46315 -0.772    -5.76  FALSE FALSE
  36.  8 AT3G55040  0.0170   -8.18  TRUE  FALSE
  37.  9




    
 AT1G56410  0.437     3.74  TRUE  TRUE
  38. 10 AT3G05360 -0.748    -9.80  FALSE FALSE

在分成了4组之后，考虑到 compareCluster 需要传入的是list，所以我们就需要把group结果直接切割成含有4组的list。这时候就要用到group_split了 

  1. test_data %>% 
  2.   mutate(r = rho > 0,
  3.          l = log2Fold > 0) %>% 
  4.   group_split(r,l)

  6. [[1]]
  7. # A tibble: 5 x 5
  8.   gene         rho log2Fold r     l
  9.               
  10. 1 AT5G07470 -0.994  -10     FALSE FALSE
  11. 2 AT5G19140 -0.291   -0.101




    
 FALSE FALSE
  12. 3 AT5G54070 -0.656   -4.14  FALSE FALSE
  13. 4 AT5G46315 -0.772   -5.76  FALSE FALSE
  14. 5 AT3G05360 -0.748   -9.80  FALSE FALSE

  16. [[2]]
  17. # A tibble: 2 x 5
  18.   gene         rho log2Fold r     l
  19.             
  20. 1 AT5G01640 -0.970     3.33 FALSE TRUE
  21. 2 CH5       -0.674     1.52 FALSE TRUE

  23. [[3]]
  24. # A tibble: 2 x 5
  25.   gene         rho log2Fold r     l
  26.               
  27. 1




    
 EMB79     0.0290    -1.52 TRUE  FALSE
  28. 2 AT3G55040 0.0170    -8.18 TRUE  FALSE

  30. [[4]]
  31. # A tibble: 1 x 5
  32.   gene        rho log2Fold r     l
  33.              
  34. 1 AT1G56410 0.437     3.74 TRUE  TRUE

每个列表的名字可以由group_keys来查看，然后你手动names一下，更改下列表元素的名字就行了。 

  1. # 可以看到分别是
  2. > test_data %>% 
  3. +   mutate(r = rho > 0,
  4. +          l = log2Fold > 0) %>% 
  5. +   group_keys(r,l)
  6. # A tibble: 4 x 2
  7.   r     l
  8.    
  9. 1 FALSE FALSE
  10. 2 FALSE TRUE
  11. 3 TRUE  FALSE
  12. 4 TRUE  TRUE

但还有个问题就是我们得到了列表内元素还是tibble，并不是我们想要的基因，所以我们还要进行一步操作lappy的操作 

  1. # 这样我们得到的就是一个含有4组元素，且元素里面都是基因的列表了
  2. # 然后就可以传入compareCluster批量做GO了。
  3. test_data %>% 
  4.   mutate(r = rho > 0,
  5.          l = log2Fold > 0) %>% 
  6.   group_split(r,l) %>% 
  7.   lapply(.,function(x) unique(pull(x,gene)))








请到「今天看啥」查看全文


推荐文章
罗永浩  ·  36氪看锤子科技：权衡之下，意外之喜
8 年前
携程每日特惠  ·  浙江居然藏着一个小“北海道”！ 6万平雪场，还有竹林飘雪温泉！（内含重磅福利）
8 年前
儿童摄影引导技巧  ·  拍婚纱照，得罪谁也不要得罪婚纱摄影师啊！
7 年前
哲学园  ·  你是一个脑还是两个脑？ | 改变心理学的40项研究
7 年前
科创通  ·  巴适丨在成都市搞创新创业活动还有补助给你拿！
7 年前
Sov5搜索   ·   小百科   ·   今天看啥   ·   移动版
51好读 - 好文章就要读起来!