💡专注R语言在🩺生物医学中的使用
设为“星标”，精彩不错过
前面给大家介绍了这么多的富集分析，其实主要就是两种：ORA和GSEA。通常都是需要一个基因集才可以做。
单个基因能做富集分析吗？肯定是不行的，所以需要我们用间接的方法实现。
对于单基因，你如果要做富集分析，有两种思路：
这个思路同样也适用于其他分子，比如lncRNA，比如miRNA（miRNA其实应该是找靶基因做，这样更合理）。
下面我们进行演示，我们选择HOPX这个基因，来自一篇文章：https://doi.org/10.1186/s12935-023-02962-2
首先我们从TCGA下载黑色素瘤的转录组数据，使用easyTCGA，1行代码解决，即可得到6种表达矩阵和临床信息，而且是官网最新的数据：
library(easyTCGA)
getmrnaexpr("TCGA-SKCM")
加载数据：
load(file = "G:/easyTCGA_test/output_mRNA_lncRNA_expr/TCGA-SKCM_mrna_expr_tpm.rdata")
这个数据是直接从GDC的官网数据中提取出来的，没有经过任何转化，所以我们先进行log2转换。
expr <- log2(mrna_expr_tpm+1)
dim(expr)
## [1] 19938 473
一共有19938个mRNA和473个样本。
提取下这个HOPX的表达矩阵看看：
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
批量计算HPOX和其他所有mRNA的相关性和P值，你自己写的循环太慢了，所以我这里推荐一种更快的方法，基于WGCNA，不过是借助linkET包实现的，这个方法我们在之前的免疫浸润中也讲过：免疫浸润结果可视化
# 自定义一个函数
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)
}
这个函数接受两个表达矩阵，然后返回相关系数和P值，不过要确认你的两个表达矩阵的样本顺序是一样的：
# 确保两个两个矩阵样本顺序一样
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
这个速度还是很快的，我还没找到比这更快的方法！
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
有些结果是NA，因为有的分子可能表达量在所有样本中都一样！我们把NA去掉即可。
最后筛选P值小于0.05和相关系数大于0.7的mRNA（这个东西没有标准，只要你能解释得通就行！）
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
根据这个结果得到82个mRNA，然后对这82个mRNA进行富集分析即可，不过我们就不演示了，因为富集分析在之前已经详细介绍过了！
我们这里根据HOPX表达量中位数进行分组，把所有样本分为高表达组和低表达组。
然后进行差异分析，这里也是用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的结果
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个是下调的，310个上调的，这个下调的很有意思！我对它很好奇，让我们看看它是谁！
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个差异基因进行富集分析了，我们还是演示下吧，进行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.kegg.jp/link/hsa/pathway"...
## Reading KEGG annotation online: "https://rest.kegg.jp/list/pathway/hsa"...
顺手再画个图：
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"))
什么？你还不会富集分析？赶紧翻看万字长文，带你彻底了解富集分析：富集分析常见类型
联系客服
