前面給大家介紹了這么多的富集分析，其實(shí)主要就是兩種：ORA和GSEA。通常都是需要一個(gè)基因集才可以做,。
單個(gè)基因能做富集分析嗎,？肯定是不行的，所以需要我們用間接的方法實(shí)現(xiàn),。
對(duì)于單基因,，你如果要做富集分析，有兩種思路：
- 批量計(jì)算和這個(gè)基因相關(guān)的其他基因,，把其他基因進(jìn)行富集分析,，這個(gè)富集分析結(jié)果就可以近似的看做是單基因的結(jié)果
- 根據(jù)這個(gè)基因的表達(dá)量進(jìn)行分組，然后做差異分析,，用差異基因做富集分析,，這個(gè)富集結(jié)果，也是基于單基因的富集
這個(gè)思路同樣也適用于其他分子,，比如lncRNA,，比如miRNA（miRNA其實(shí)應(yīng)該是找靶基因做,，這樣更合理）。
下面我們進(jìn)行演示,，我們選擇HOPX這個(gè)基因,，來自一篇文章：https:///10.1186/s12935-023-02962-2
數(shù)據(jù)準(zhǔn)備
首先我們從TCGA下載黑色素瘤的轉(zhuǎn)錄組數(shù)據(jù)，使用easyTCGA,，1行代碼解決,，即可得到6種表達(dá)矩陣和臨床信息，而且是官網(wǎng)最新的數(shù)據(jù)：
library(easyTCGA)
getmrnaexpr("TCGA-SKCM")
加載數(shù)據(jù)：
load(file = "G:/easyTCGA_test/output_mRNA_lncRNA_expr/TCGA-SKCM_mrna_expr_tpm.rdata")
這個(gè)數(shù)據(jù)是直接從GDC的官網(wǎng)數(shù)據(jù)中提取出來的,，沒有經(jīng)過任何轉(zhuǎn)化,，所以我們先進(jìn)行l(wèi)og2轉(zhuǎn)換。
expr <- log2(mrna_expr_tpm+1)
dim(expr)
## [1] 19938 473
一共有19938個(gè)mRNA和473個(gè)樣本,。
提取下這個(gè)HOPX的表達(dá)矩陣看看：
HOPX_expr <- expr["HOPX",]
HOPX_expr[,1:4]
## TCGA-EB-A3Y6-01A-21R-A239-07 TCGA-D9-A4Z6-06A-12R-A266-07
## HOPX 1.198746 0.4733194
## TCGA-FW-A5DY-06A-11R-A311-07 TCGA-EE-A2GH-06A-11R-A18T-07
## HOPX 1.760987 2.010207
相關(guān)性分析
批量計(jì)算HPOX和其他所有mRNA的相關(guān)性和P值,，你自己寫的循環(huán)太慢了，所以我這里推薦一種更快的方法,，基于WGCNA,，不過是借助linkET包實(shí)現(xiàn)的，這個(gè)方法我們?cè)谥暗拿庖呓?rùn)中也講過：免疫浸潤(rùn)結(jié)果可視化
# 自定義一個(gè)函數(shù)
cormatrixes <- function(x,y){
tem <- linkET::correlate(t(x),t(y),engine = "WGCNA")
tem1 <- as.data.frame(tem[[1]])
tem1 <- cbind(rownames(tem1),tem1)
tem1_long <- reshape2::melt(tem1,value.name = "correlation")
tem2 <- as.data.frame(tem[[2]])
tem2 <- cbind(rownames(tem2),tem2)
tem2_long <- reshape2::melt(tem2,value.name = "pvalue")
result <- cbind(tem1_long,tem2_long$pvalue)
names(result) <- c("v1","v2","correlation","pvalue")
return(result)
}
這個(gè)函數(shù)接受兩個(gè)表達(dá)矩陣,，然后返回相關(guān)系數(shù)和P值,，不過要確認(rèn)你的兩個(gè)表達(dá)矩陣的樣本順序是一樣的：
# 確保兩個(gè)兩個(gè)矩陣樣本順序一樣
identical(colnames(expr),colnames(HOPX_expr))
## [1] TRUE
cor_res <- cormatrixes(HOPX_expr,expr)
##
## Using rownames(tem1) as id variables
## Using rownames(tem2) as id variables
這個(gè)速度還是很快的，我還沒找到比這更快的方法,！
head(cor_res)
## v1 v2 correlation pvalue
## 1 HOPX MT-CO2 -0.08073411 0.07941864
## 2 HOPX MT-CO3 -0.05334072 0.24692973
## 3 HOPX MT-ND4 -0.08936978 0.05208752
## 4 HOPX MT-CO1 -0.06033840 0.19019751
## 5 HOPX MT-CYB -0.05946605 0.19669609
## 6 HOPX MT-ATP6 -0.05107928 0.26756587
有些結(jié)果是NA,，因?yàn)橛械姆肿涌赡鼙磉_(dá)量在所有樣本中都一樣！我們把NA去掉即可,。
最后篩選P值小于0.05和相關(guān)系數(shù)大于0.7的mRNA（這個(gè)東西沒有標(biāo)準(zhǔn),，只要你能解釋得通就行！）
cor_res <- na.omit(cor_res)
suppressMessages(library(dplyr))
hopx_related_mrna <- cor_res %>%
filter(correlation > 0.7, pvalue < 0.05) %>%
distinct(v2) %>%
pull(v2)
length(hopx_related_mrna)
## [1] 82
head(hopx_related_mrna)
## [1] CSTA AQP3 TACSTD2 TUBA4A SLURP1 ST14
## 19938 Levels: MT-CO2 MT-CO3 MT-ND4 MT-CO1 MT-CYB MT-ATP6 FTL MT-ND2 ... AC006486.3
根據(jù)這個(gè)結(jié)果得到82個(gè)mRNA,，然后對(duì)這82個(gè)mRNA進(jìn)行富集分析即可,，不過我們就不演示了，因?yàn)楦患治鲈谥耙呀?jīng)詳細(xì)介紹過了,！
根據(jù)表達(dá)量分組
我們這里根據(jù)HOPX表達(dá)量中位數(shù)進(jìn)行分組,，把所有樣本分為高表達(dá)組和低表達(dá)組。
然后進(jìn)行差異分析,，這里也是用easyTCGA1行代碼解決：
sample_group <- ifelse(expr["HOPX",] > median(t(expr["HOPX",])), "high","low")
sample_group <- factor(sample_group, levels = c("low","high"))
library(easyTCGA)
deg_res <- diff_analysis(exprset = expr,
group = sample_group,
is_count = F,
logFC_cut = 1,
pvalue_cut = 0.05,
adjpvalue_cut = 0.05,
save = F
)
## => log2 transform not needed
## => Running limma
## => Running wilcoxon test
## => Analysis done.
# 提取limma的結(jié)果
deg_limma <- deg_res$deg_limma
head(deg_limma)
## logFC AveExpr t P.Value adj.P.Val B
## HOPX 1.411491 1.330710 17.92159 3.648490e-55 7.274359e-51 114.22524
## IL18 1.579073 2.882566 15.19126 1.021080e-42 1.017914e-38 86.06859
## TNFSF10 1.627203 3.868731 14.82989 4.077339e-41 2.709799e-37 82.44504
## DAPP1 1.120535 1.551681 13.49540 2.489357e-35 1.240820e-31 69.35296
## ANKRD22 1.459943 1.975013 13.06324 1.661069e-33 6.623677e-30 65.22542
## ZBED2 1.241665 1.613864 12.12598 1.202387e-29 3.995531e-26 56.49488
## genesymbol
## HOPX HOPX
## IL18 IL18
## TNFSF10 TNFSF10
## DAPP1 DAPP1
## ANKRD22 ANKRD22
## ZBED2 ZBED2
只有1個(gè)是下調(diào)的,，310個(gè)上調(diào)的，這個(gè)下調(diào)的很有意思,！我對(duì)它很好奇,，讓我們看看它是誰(shuí)！
deg_limma %>%
dplyr::filter(logFC < 0)
## logFC AveExpr t P.Value adj.P.Val B
## VGF -1.197905 5.679119 -3.760606 0.0001907947 0.0008006871 0.02293662
## genesymbol
## VGF VGF
接下來就可以使用這311個(gè)差異基因進(jìn)行富集分析了,，我們還是演示下吧,，進(jìn)行GO和KEGG的富集分析：
suppressMessages(library(clusterProfiler))
deg_entrezid <- bitr(deg_limma$genesymbol,fromType = "SYMBOL"
,toType = "ENTREZID",OrgDb = "org.Hs.eg.db")
##
## 'select()' returned 1:1 mapping between keys and columns
# GO
go_res <- enrichGO(deg_entrezid$ENTREZID,
OrgDb = "org.Hs.eg.db",
ont = "ALL",
readable = T
)
# KEGG
kegg_res <- enrichKEGG(deg_entrezid$ENTREZID)
## Reading KEGG annotation online: "https://rest./link/hsa/pathway"...
## Reading KEGG annotation online: "https://rest./list/pathway/hsa"...
順手再畫個(gè)圖：
library(enrichplot)
p1 <- dotplot(go_res,label_format=50,showCategory=20)
p2 <- dotplot(kegg_res,label_format=50,showCategory=20)
cowplot::plot_grid(p1,p2,labels = c("A","B"))
plot of chunk unnamed-chunk-12什么,？你還不會(huì)富集分析？趕緊翻看萬字長(zhǎng)文,，帶你徹底了解富集分析：富集分析常見類型
