内容

1理论

看到幻灯片．

最近的一次推特提供了一个很好的总结的努力基准基因集富集分析方法使用GSEABenchmarkR包中。

2实践

数据输入和按摩

库(气道)数据(气道)气道$dex <-重新level(气道$dex， "untrt")

微分表达式分析

library(DESeq2) des <- DESeqDataSet(气道，design = ~ cell + dex) des <- DESeq(des) res <- results(des)

过渡到整齐的数据

library(dplyr) library(tibble) tbl <- res %>% as.data.frame() %>% as_tibble(rownames = "ENSEMBL") tbl . dplyr
## # A tibble: 64102 x 7 ## ENSEMBL baseMean log2FoldChange lfcSE stat pvalue padj ##        ## 1 ENSG0000000…709。-0.381 0.101 -3.79 1.52e-4 1.28e-3 ## 2 ENSG0000000…0 NA NA NA NA NA ## 3 ENSG0000000…520。0.207 0.112 1.84 6.53e-2 1.97e-1 ## 4 ENSG0000000…0.0379 0.143 0.264 7.92e-1 9.11e-1 ## 5 ENSG0000000…57.9 -0.0882 0.287 -0.307 7.59e-1 8.95e-1 ## 6 ENSG0000000…0.426 0.0883 4.83 1.38e-6 1.82e-5 ## 8 ENSG0000000…1282。-0.241 0.0887 -2.72 6.58e-3 3.28e-2 ## 9 ENSG0000000…-0.0476 0.167 -0.286 7.75e-1 9.053 e-1 ## 10 ENSG0000000…-0.500 0.121 -4.14 3.48e-5 3.42e-4 ## #

2.1示例:使用超几何测试limma:: goana ()

需要ENTREZ标识符

library(org. hs . e.g. .db) tbl <- tbl %>% mutate(ENTREZID = mapIds(org. hs . e.g. .db, ENSEMBL， "ENTREZID"， "ENSEMBL")) %>% dplyr::select(ENSEMBL, ENTREZID, everything()) tbl
## # A tibble: 64102 x 8 ## ENSEMBL ENTREZID baseMean log2FoldChange lfcSE stat pvalue ##        ## 1 ENSG00…7105 709。-0.381 0.101 -3.79 1.52e-4 ## 2 ENSG00…64102 0 NA NA NA NA ## 3 ENSG00…8813 520。0.207 0.112 1.84 6.53e-2 ## 4 ENSG00…57147 237。0.0379 0.143 0.264 7.92e-1 ## 5 ENSG00…55732 57.9 -0.0882 0.287 -0.307 7.59e-1 ## 6 ENSG00…2268 0.318 -1.38 3.50 -0.394 6.94e-1 ## 7 ENSG00…3075 5817。0.426 0.0883 4.83 1.38e-6 ## 8 ENSG00…2519 1282。-0.241 0.0887 -2.72 6.58e-3 ## 9 ENSG00…2729 610。-0.0476 0.167 -0.286 7.75e-1 ## 10 ENSG00…4800 369。-0.500 0.121 -4.14 3.48e-5 ## #…多了64,092行，多了一个变量:padj

宇宙-必须是可测试的DE

tbl <- tbl %>% filter(!is.na(padj)， !is.na(ENTREZID)
## # A tibble: 14,550 x 8 ## ENSEMBL ENTREZID baseMean log2FoldChange lfcSE stat pvalue padj# #         ## 1 ENSG0000…7105 709。-0.381 0.101 -3.79 1.52e-4 1.28e-3 ## 2 ENSG0000…0.207 0.112 1.84 6.53e-2 1.97e-1 ## 3 ENSG0000…57147 237。0.0379 0.143 0.264 7.92e-1 9.11e-1 ## 4 ENSG0000…55732 57.9 -0.0882 0.287 -0.307 7.59e-1 8.95e-1 ## 5 ENSG0000…3075 5817。0.426 0.0883 4.83 1.38e-6 1.82e-5 ## 6 ENSG0000…2519 1282。-0.241 0.0887 -2.72 6.58e-3 3.28e-2 ## 7 ENSG0000…2729 610。-0.0476 0.167 -0.286 7.75e-1 9.032 e-1 ## 8 ENSG0000…4800 369。-0.500 0.121 -4.14 3.48e-5 3.42e-4 ## 9 ENSG0000…90529 183。-0.124 0.180 -0.689 4.91e-1 7.24e-1 ## 10 ENSG0000…57185 2814。-0.0411 0.103 -0.400 6.89e-1 8.57e-1 ## #

limma:: goana ()——超几何

library(limma) go <- goana(tbl$ENTREZID[tbl$padj < .05]， tbl$ENTREZID， "Hs") %>% as_tibble(rownames = "GOID") %>% dplyr::select(GOID, Ont, everything()) go
## #小猫咪:20656 x 6 # # GOID安大略省的术语N DE P.DE # # <空空的> <空空的> <空空的> <双> <双> <双> # # 1:00017……BP细胞激活939 333 6.80 e-12 # # 2: 00022……BP免疫效应处理808 251 2.98 e - 4 # # 3: 00022 BP参与免疫细胞激活…498 148 2.43 e - 2 # # 4: 00022……BP骨髓白细胞激活463 147 2.01 e - 3 # # 5: 00022 BP骨髓细胞的激活参与我……399 118 4.61 e - 2 # # 6: 00022 BP中性粒细胞激活参与……361 108 4.08 e - 2 # # 7: 00023……BP白细胞激活涉及i…494 148 1.86e- 2 ## 8 GO:00023…BP免疫系统过程1958 652 8.14e-16 ## 9 GO:00024…BP白细胞介导的免疫538 165 5.33e- 3 ## 10 GO:00024…BP髓系白细胞介导的免疫408 124 1.89e- 2 ## #…与20,646多行

哪些GO id的表达差异最大?(为什么这些基因集看起来很大，N？)

安排(P.DE)
## #小猫咪:20656 x 6 # # GOID安大略省的术语N DE P.DE # # <空空的> <空空的> <空空的> <双> <双> <双> # # 1:3081 CC细胞外围0071944 1104 8.45 e-45 # # 2: 5799 BP应对刺激0050896 1859 8.77 e-45 # # 3: 0032502 BP发展过程4168 1410 9.73 e-44 # # 4: 0048856 BP解剖结构发展3891 1330 3.45 e-43 # # 5: 0005886 CC质膜3013 1078 3.83 e-43 # # 6: 0032501 BP的多细胞有机体的过程4678 1546 1.83 e-42 # # 7: 0048731 3202 1128 BP系统开发7.52e-42 ## 8 GO:0023052 BP信号4060 1372 7.66e-42 ## 9 GO:0007275 BP多细胞生物发育3579 1231 1.46e-40 ## 10 GO:0007154 BP细胞通信4076 1371 1.57e-40 ## #

检查?

go %>% filter(grepl("糖皮质激素"，术语))%>% arrange(P.DE)
## #小猫咪:22 x 6 # # GOID一项N DE P.DE # # <空空的> <空空的> <空空的> <双> <双> <双> # # 1:0051年…92年英国石油公司对糖皮质激素43 1.11 e-5 # # 2: 0071 BP细胞对糖皮质激素…42 21 6.42军医# # 3:2000 BP积极监管glucocortico…2 2 6.64依照# # 4:2000 BP糖皮质激素调节的感受…8 4 1.25 e 1 # # 5: 0006……BP糖皮质激素生物合成过程8 4 1.25 e 1 # # 6: 0043 BP糖皮质激素介导信号p…3 2 1.65 e 1 # # 7: 0008……BP糖皮质激素代谢过程13 5 2.26e-1 ## 8 GO:0004…MF糖皮质激素受体活性1 1 2.58e-1 ## 9 GO:0031…BP糖皮质激素负调控…4 2 2.75e-1 ## # 10 GO:0031…BP糖皮质激素负调控…4 2 2.75e-1 ## #

是什么基因?

genesets <- AnnotationDbi::select(org.Hs.eg.db, tbl$ENTREZID， "GOALL"， "ENTREZID") %>% as_tibble() %>% dplyr::select(GOID = GOALL, Ont = ONTOLOGYALL, ENTREZID) %>% distinct() %>% arrange(Ont, GOID) genessets
## 1 GO:0000002 BP 3980 ## 2 GO:0000002 BP 1890 ## 3 GO:0000002 BP 50484 ## 4 GO:0000002 BP 4205 ## 5 GO:0000002 BP 9093 ## 6 GO:0000002 BP 56652 ## 7 GO:0000002 BP 10891 ## 9 GO:0000002 BP 4358 ## 10 GO:0000002 BP 10000 ## ##…与1,574,070多行

3.出处

sessionInfo ()
## R version 3.6.1 Patched (2019-07-16 r76845) ##平台:x86_64-apple-darwin17.7.0(64位)##运行在macOS High Sierra 10.13.6 ## ##矩阵产品:default ## BLAS: /Users/ma38727/bin/R-3-6-branch/lib/libRblas。/Users/ma38727/bin/R-3-6-branch/lib/libRlapack。dylib # # # #语言环境:# # [1]en_US.UTF-8 / en_US.UTF-8 en_US.UTF-8 / C / en_US.UTF-8 / en_US。UTF-8 ## ##附加的基本包:## [1]stats4并行统计图形grDevices utils数据集##[8]方法基础## ##其他附加包:# # # # [1] tibble_2.1.3 limma_3.41.14 [3] GO.db_3.8.2 org.Hs.eg.db_3.8.2 # # [5] AnnotationDbi_1.47.0 dplyr_0.8.3 # # [7] airway_1.5.0 DESeq2_1.25.7 # # [9] EnrichmentBrowser_2.15.5 graph_1.63.0 # # [11] SummarizedExperiment_1.15.5 DelayedArray_0.11.4 # # [13] BiocParallel_1.19.0 matrixStats_0.54.0 # # [15] Biobase_2.45.0 GenomicRanges_1.37.14 # # [17] GenomeInfoDb_1.21.1 IRanges_2.19.10 # # [19] S4Vectors_0.23.17 BiocGenerics_0.31.5 # # [21] BiocStyle_2.13.2 # # # #加载(而不是通过一个名称空间附):## [1] bitops_1.0-6 bit64_0.9-7 RColorBrewer_1.1-2 ## [4] tools_3.6.1 backports_1.1.4 utf8_1 - 1.4 ## [7] R6_2.4.0 rpart_4.1-15 Hmisc_4.2-0 ## [10] DBI_1.0.0 lazyeval_0.2.2 colorspace_1.4-1 ## [13] nnet_7.3-12 tidyselect_0.2.5 gridextra2.3 ## [16] bit_1.1-14 compiler_3.6.1 cli_1.1.0 ## [19] htmlTable_1.13.1 bookdown_0.12 KEGGgraph_1.45.1 ## [22] scales_1.0.0 checkmate_1.9.4 genefilter_1.67.1 ## [25] rappdirs_0.3.1 string_1 .4.0 digest_0.6.20 ## [28] foreign_0.8-71 rmarkdown_1.14 XVector_0.25.0 ## #[31] base64enc_0.1-3 pkgconfig_2.0.2 htmltools_0.3.6 ## [34] htmlwidgets_1.3 rlang_0.4.0 rstudioapi_0.10 ## [37] RSQLite_2.1.2 acepack_1.4.1 RCurl_1.95-4.12 ## [40] magrittr_1.5 GenomeInfoDbData_1.2.1 Formula_1.2-3 ## [43] Matrix_1.2-17 fansi_0.4.0 Rcpp_1.0.1 ## [46] munsell_0.5.0 stringi_1.4.3 yaml_2.2.0 ## [49] zlibbioc_1.31.0 grid_3.6.1 blob_1.2.0 ## [52] crayon_1.3.4 lattice_0.20-38 splines_3.6.1 ## [55] annotate_1.63.0 locfit_1.5-9.1 zeallot_0.1.0 ## [58] knitr_1.23 pillar_1.4.2 codetools_0.2-16 ## [61] geneplotter_1.63.0 XML_3.98-1.20 glue_1.3.1 ## [64] evaluate_0.14 latticeExtra_0.6-28 data.table_1.12.2 ## [67] BiocManager_1.30.4 vctrs_0.2.0 gtable_0.3.0 ## [70] purrr_0.3.2 assertthat_0.2.1 ggplot2_3.2.0 ## [73] xfun_0.8 xtable_1.8-4 survival_2.44-1.1 ## [76] memoise_1.1.0 cluster_2.1.0 GSEABase_1.47.0