github-actions[bot]
commented
4 years ago 使用fgsea包在60秒内获得基因调控信号通路 by 医知圈
点击标题右下方“医知圈”关注本公号
在这里，我将向您展示如何从DESeq获取结果，并使用fgsea软件包快速评估通路的富集，以进行快速预先基因集富集分析（GSEA）。
关于DESeq结果的FGSEA（人类）
首先，让我们从人类数据集开始。
从DESeq2加载一些结果
从DESeq2加载一些结果（使用生成results(..., tidy=TRUE)）。有关如何生成此数据，请参见最后的附录部分。
library(tidyverse)res <- read_csv("data/deseq-results-tidy-human-airway.csv")res## # A tibble: 64,102 x 7## row baseMean log2FoldChange lfcSE stat pvalue padj#### 1 ENSG000000… 709. 0.381 0.101 3.79 1.52e-4 1.28e-3## 2 ENSG000000… 0 NA NA NA NA NA ## 3 ENSG000000… 520. -0.207 0.112 -1.84 6.53e-2 1.97e-1## 4 ENSG000000… 237. -0.0379 0.143 -0.264 7.92e-1 9.11e-1## 5 ENSG000000… 57.9 0.0882 0.287 0.307 7.59e-1 8.95e-1## 6 ENSG000000… 0.318 1.38 3.50 0.394 6.94e-1 NA ## 7 ENSG000000… 5817. -0.426 0.0883 -4.83 1.38e-6 1.82e-5## 8 ENSG000000… 1282. 0.241 0.0887 2.72 6.58e-3 3.28e-2## 9 ENSG000000… 610. 0.0476 0.167 0.286 7.75e-1 9.03e-1## 10 ENSG000000… 369. 0.500 0.121 4.14 3.48e-5 3.42e-4## # ... with 64,092 more rows
将Ensembl基因ID映射到符号。首先创建一个映射表。
library(org.Hs.eg.db)ens2symbol <- AnnotationDbi::select(org.Hs.eg.db, key=res$row, columns="SYMBOL", keytype="ENSEMBL")ens2symbol <- as_tibble(ens2symbol)ens2symbol
## # A tibble: 64,381 x 2## ENSEMBL SYMBOL #### 1 ENSG00000000003 TSPAN6 ## 2 ENSG00000000005 TNMD ## 3 ENSG00000000419 DPM1 ## 4 ENSG00000000457 SCYL3 ## 5 ENSG00000000460 C1orf112## 6 ENSG00000000938 FGR ## 7 ENSG00000000971 CFH ## 8 ENSG00000001036 FUCA2 ## 9 ENSG00000001084 GCLC ## 10 ENSG00000001167 NFYA ## # ... with 64,371 more rows
现在加入他们。
res <- inner_join(res, ens2symbol, by=c("row"="ENSEMBL"))res## # A tibble: 64,381 x 8## row baseMean log2FoldChange lfcSE stat pvalue padj SYMBOL#### 1 ENSG0… 709. 0.381 0.101 3.79 1.52e-4 1.28e-3 TSPAN6## 2 ENSG0… 0 NA NA NA NA NA TNMD ## 3 ENSG0… 520. -0.207 0.112 -1.84 6.53e-2 1.97e-1 DPM1 ## 4 ENSG0… 237. -0.0379 0.143 -0.264 7.92e-1 9.11e-1 SCYL3 ## 5 ENSG0… 57.9 0.0882 0.287 0.307 7.59e-1 8.95e-1 C1orf…## 6 ENSG0… 0.318 1.38 3.50 0.394 6.94e-1 NA FGR ## 7 ENSG0… 5817. -0.426 0.0883 -4.83 1.38e-6 1.82e-5 CFH ## 8 ENSG0… 1282. 0.241 0.0887 2.72 6.58e-3 3.28e-2 FUCA2 ## 9 ENSG0… 610. 0.0476 0.167 0.286 7.75e-1 9.03e-1 GCLC ## 10 ENSG0… 369. 0.500 0.121 4.14 3.48e-5 3.42e-4 NFYA ## # ... with 64,371 more rows
此外，您以后所关心的只是基因符号和检验统计量。获取那些，并删除NAs。最后，如果您有相同符号的多个测试统计信息，则需要以某种方式处理该问题。在这里，我只是将它们平均化。
res2 <- res %>% dplyr::select(SYMBOL, stat) %>% na.omit() %>% distinct() %>% group_by(SYMBOL) %>% summarize(stat=mean(stat))res2
## # A tibble: 19,662 x 2## SYMBOL stat#### 1 A1BG 0.680 ## 2 A1BG-AS1 -1.79 ## 3 A1CF -0.126 ## 4 A2M -1.26 ## 5 A2M-AS1 0.875 ## 6 A2ML1 0.673 ## 7 A2MP1 0.973 ## 8 A3GALT2 -0.0491 ## 9 A4GALT -4.14 ## 10 A4GNT 0.00777## # ... with 19,652 more rows
使用fgsea包
我们将使用fgsea软件包进行快速预分解基因集富集分析（GSEA）。
基因集富集分析是用于分析基因表达数据的广泛使用的工具。然而，由于用于分析的大量所需样本具有良好的统计功率，当前的实施方案是缓慢的。在本文中，我们提出了一种新算法，可以有效地多次重复使用一个样本，从而加快分析速度。我们表明，可以在几分钟内产生数十万个排列，从而产生非常准确的p值。反过来，这允许应用标准FDR校正程序，这些程序比当前使用的程序更准确。
library(fgsea)
该fgsea()功能需要检查基因集列表，以及基因水平统计的命名向量，其中名称应与通路列表中的基因名称相同。首先，让我们创建测试统计的命名向量。请参阅?tibble::deframe此处的帮助 - deframe()将两列数据帧转换为命名向量或列表，使用第一列作为名称，第二列作为值。
ranks <- deframe(res2)head(ranks, 20)
## A1BG A1BG-AS1 A1CF A2M A2M-AS1 ## 0.679946437 -1.793291412 -0.126192849 -1.259539478 0.875346116 ## A2ML1 A2MP1 A3GALT2 A4GALT A4GNT ## 0.672710814 0.973319532 -0.049131427 -4.144839902 0.007772497 ## AAAS AACS AADAC AADACL2 AADACL4 ## 0.163986128 1.416071728 -0.306541616 0.489315094 -1.876311694 ## AADACP1 AADAT AAED1 AAGAB AAK1 ## 0.161119304 3.079128034 -7.492301790 1.554279946 1.141522348
让我们使用MSigDB中的Hallmark基因集。Hallmark基因集总结并代表特定的明确定义的生物状态或过程并显示连贯表达。这些基因集是通过基于识别其他MSigDB集合中基因集之间的重叠和保留显示坐标表达的基因的计算方法产生的。该gmtPathways()函数将获取您从MSigDB下载的GMT文件并将其转换为列表。列表中的每个元素是该途径中基因的特征向量。
# Load the pathways into a named listpathways.hallmark <- gmtPathways("data/msigdb/h.all.v6.2.symbols.gmt")# Look at them all if you want (uncomment)# pathways.hallmark# Show the first few pathways, and within those, show only the first few genes. pathways.hallmark %>% head() %>% lapply(head)## $HALLMARK_TNFA_SIGNALING_VIA_NFKB## [1] "JUNB" "CXCL2" "ATF3" "NFKBIA" "TNFAIP3" "PTGS2" ## ## $HALLMARK_HYPOXIA## [1] "PGK1" "PDK1" "GBE1" "PFKL" "ALDOA" "ENO2" ## ## $HALLMARK_CHOLESTEROL_HOMEOSTASIS## [1] "FDPS" "CYP51A1" "IDI1" "FDFT1" "DHCR7" "SQLE" ## ## $HALLMARK_MITOTIC_SPINDLE## [1] "ARHGEF2" "CLASP1" "KIF11" "KIF23" "ALS2" "ARF6" ## ## $HALLMARK_WNT_BETA_CATENIN_SIGNALING## [1] "MYC" "CTNNB1" "JAG2" "NOTCH1" "DLL1" "AXIN2" ## ## $HALLMARK_TGF_BETA_SIGNALING## [1] "TGFBR1" "SMAD7" "TGFB1" "SMURF2" "SMURF1" "BMPR2"
现在，使用1000个排列运行fgsea算法：
fgseaRes <- fgsea(pathways=pathways.hallmark, stats=ranks, nperm=1000)
整理结果：
fgseaResTidy <- fgseaRes %>% as_tibble() %>% arrange(desc(NES))# Show in a nice table:fgseaResTidy %>% dplyr::select(-leadingEdge, -ES, -nMoreExtreme) %>% arrange(padj) %>% DT::datatable()路PVALpadjNES尺寸
| 1 | HALLMARK_E2F_TARGETS | 0.003125 | 0.016 | 1.50255883270908 | 196 |
| 2 | HALLMARK_P53_PATHWAY | 0.00307692307692308 | 0.016 | 1.46976609146818 | 193 |
| 3 | HALLMARK_MTORC1_SIGNALING | 0.0031055900621118 | 0.016 | 1.40459809605177 | 195 |
| 4 | HALLMARK_COMPLEMENT | 0.0029806259314456 | 0.016 | -1.500920162306 | 177 |
| 五 | HALLMARK_IL2_STAT5_SIGNALING | 0.00294985250737463 | 0.016 | -1.50488334166041 | 186 |
| 6 | HALLMARK_XENOBIOTIC_METABOLISM | 0.00294117647058824 | 0.016 | -1.52285404669765 | 191 |
| 7 | HALLMARK_OXIDATIVE_PHOSPHORYLATION | 0.0029673590504451 | 0.016 | -1.61712192519081 | 183 |
| 8 | HALLMARK_ANDROGEN_RESPONSE | 0.0032 | 0.016 | -1.69156156907097 | 95 |
| 9 | HALLMARK_TNFA_SIGNALING_VIA_NFKB | 0.00147275405007364 | 0.016 | -1.77543957704369 | 192 |
| 10 | HALLMARK_ADIPOGENESIS | 0.00147492625368732 | 0.016 | -1.87910665606575 | 186 |
显示50个参赛作品中的1到10个
上一页12345下一页
绘制标准化的浓缩分数。为指示通路是否重要的条添加颜色：
ggplot(fgseaResTidy, aes(reorder(pathway, NES), NES)) + geom_col(aes(fill=padj<0.05)) + coord_flip() + labs(x="Pathway", y="Normalized Enrichment Score", title="Hallmark pathways NES from GSEA") + theme_minimal()
该α干扰素响应基因组被显著上调，与调查结果一致最初的研究。
每种途径中有哪些基因？首先，获得所有途径及其中的基因。继续将其加入到原始数据中以从中获取基因。（可选）筛选列表以仅包含重要的列表等。
pathways.hallmark %>% enframe("pathway", "SYMBOL") %>% unnest() %>% inner_join(res, by="SYMBOL")## # A tibble: 7,689 x 9## pathway SYMBOL row baseMean log2FoldChange lfcSE stat pvalue#### 1 HALLMARK… JUNB ENSG… 682. 0.0640 0.144 0.445 6.56e- 1## 2 HALLMARK… CXCL2 ENSG… 1.71 1.39 1.55 0.899 3.69e- 1## 3 HALLMARK… ATF3 ENSG… 374. -1.92 0.159 -12.1 1.56e-33## 4 HALLMARK… NFKBIA ENSG… 562. -0.799 0.177 -4.51 6.37e- 6## 5 HALLMARK… TNFAI… ENSG… 356. -0.562 0.149 -3.79 1.53e- 4## 6 HALLMARK… PTGS2 ENSG… 71.6 1.37 0.365 3.76 1.70e- 4## 7 HALLMARK… CXCL1 ENSG… 37.3 -0.841 0.438 -1.92 5.47e- 2## 8 HALLMARK… IER3 ENSG… 136. 1.30 0.217 5.99 2.08e- 9## 9 HALLMARK… IER3 ENSG… 0 NA NA NA NA ## 10 HALLMARK… IER3 ENSG… 0 NA NA NA NA ## # ... with 7,679 more rows, and 1 more variable: padj
让我们尝试一组不同的途径。让我们来看看KEGG的途径：
fgsea(pathways=gmtPathways("data/msigdb/c2.cp.kegg.v6.2.symbols.gmt"), ranks, nperm=1000) %>% as_tibble() %>% arrange(padj)## # A tibble: 186 x 8## pathway pval padj ES NES nMoreExtreme size leadingEdge#### 1 KEGG_STARCH… 0.00176 0.0427 -0.654 -1.77 0 33
## 2 KEGG_GLYCOS… 0.00183 0.0427 -0.751 -1.87 0 21 ## 3 KEGG_FOCAL_… 0.00141 0.0427 -0.458 -1.59 0 191 ## 4 KEGG_NEUROT… 0.00152 0.0427 -0.506 -1.66 0 120 ## 5 KEGG_REGULA… 0.00141 0.0427 -0.483 -1.68 0 192 ## 6 KEGG_INSULI… 0.00151 0.0427 -0.543 -1.80 0 126 ## 7 KEGG_ADIPOC… 0.00159 0.0427 -0.646 -1.94 0 62 ## 8 KEGG_ACUTE_… 0.00162 0.0427 -0.582 -1.70 0 52 ## 9 KEGG_GLYCER… 0.00316 0.0456 -0.551 -1.67 1 65 ## 10 KEGG_AMINOA… 0.00245 0.0456 0.617 1.86 0 41 ## # ... with 176 more rows
或miR目标：
fgsea(pathways=gmtPathways("data/msigdb/c3.mir.v6.2.symbols.gmt"), ranks, nperm=1000) %>% as_tibble() %>% arrange(padj)## # A tibble: 221 x 8## pathway pval padj ES NES nMoreExtreme size leadingEdge#### 1 TGAATGT_MIR… 0.00124 0.0346 -0.433 -1.63 0 446
## 2 AATGTGA_MIR… 0.00128 0.0346 -0.405 -1.51 0 370 ## 3 CAGTGTT_MIR… 0.00136 0.0346 -0.419 -1.52 0 286 ## 4 GTGCCTT_MIR… 0.00115 0.0346 -0.399 -1.54 0 664 ## 5 CTATGCA_MIR… 0.00141 0.0346 -0.486 -1.71 0 196 ## 6 TTGCCAA_MIR… 0.00133 0.0346 -0.413 -1.51 0 291 ## 7 TTTGTAG_MIR… 0.00132 0.0346 -0.411 -1.50 0 306 ## 8 TGCCTTA_MIR… 0.00121 0.0346 -0.388 -1.47 0 519 ## 9 GTGCCAA_MIR… 0.00135 0.0346 -0.445 -1.61 0 274 ## 10 TGTTTAC_MIR… 0.00242 0.0478 -0.361 -1.37 1 532 ## # ... with 211 more rows
或GO注释：
fgsea(pathways=gmtPathways("data/msigdb/c5.all.v6.2.symbols.gmt"), ranks, nperm=1000) %>% as_tibble() %>% arrange(padj)## # A tibble: 5,917 x 8## pathway pval padj ES NES nMoreExtreme size leadingEdge#### 1 GO_CELLULAR… 0.00108 0.0870 -0.376 -1.49 0 1200
## 2 GO_CHEMICAL… 0.00115 0.0870 -0.347 -1.35 0 742 ## 3 GO_RESPONSE… 0.00128 0.0870 -0.418 -1.54 0 355 ## 4 GO_FORMATIO… 0.00154 0.0870 -0.515 -1.67 0 98 ## 5 GO_ACTIN_FI… 0.00124 0.0870 -0.434 -1.62 0 409 ## 6 GO_SINGLE_O… 0.00109 0.0870 -0.332 -1.31 0 1169 ## 7 GO_SMALL_MO… 0.00106 0.0870 -0.339 -1.35 0 1541 ## 8 GO_ALPHA_AM… 0.00153 0.0870 -0.551 -1.75 0 83 ## 9 GO_ORGANOPH… 0.00114 0.0870 -0.360 -1.40 0 767 ## 10 GO_ESTABLIS… 0.00106 0.0870 -0.328 -1.31 0 1515 ## # ... with 5,907 more rows
非人类生物
如果您有非人类生物，首先要创建一个表格，将fly / mouse / worm / etc基因ID映射到人类符号：
library(biomaRt)mart <- useDataset("mmusculus_gene_ensembl", mart=useMart("ensembl"))bm <- getBM(attributes=c("ensembl_gene_id", "hsapiens_homolog_associated_gene_name"), mart=mart) %>% distinct() %>% as_tibble() %>% na_if("") %>% na.omit()bm## # A tibble: 26,256 x 2## ensembl_gene_id hsapiens_homolog_associated_gene_name#### 1 ENSMUSG00000064370 MT-CYB ## 2 ENSMUSG00000064368 MT-ND6 ## 3 ENSMUSG00000064367 MT-ND5 ## 4 ENSMUSG00000064363 MT-ND4 ## 5 ENSMUSG00000065947 MT-ND4L ## 6 ENSMUSG00000064360 MT-ND3 ## 7 ENSMUSG00000064358 MT-CO3 ## 8 ENSMUSG00000064357 MT-ATP6 ## 9 ENSMUSG00000064356 MT-ATP8 ## 10 ENSMUSG00000064354 MT-CO2 ## # ... with 26,246 more rows
然后简单地将您的非人类结果加入到上面的表中，并使用人类符号和相关的测试统计数据继续执行之前的操作。
附录
此代码用于生成人体呼吸道结果：
library(DESeq2)library(airway)ddsSE <- DESeqDataSet(airway, design = ~ cell + dex)ddsSE <- DESeq(ddsSE)res <- results(ddsSE, tidy = TRUE)readr::write_csv(res, path="data/deseq-results-tidy-human-airway.csv")
如果您觉得有价值，请把此文放到您朋友圈，大家都会感谢你
ixxmu
https://mp.weixin.qq.com/s/P3m2YySdMFuEUrh0zkGPuw