单基因富集分析

前面给大家介绍了这么多的富集分析，其实主要就是两种：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"))

什么？你还不会富集分析？赶紧翻看万字长文，带你彻底了解富集分析：富集分析常见类型