RNA-seq：KEGG富集分析

 代码大部分来自简书：猪莎爱学习，部分修改过

1.数据准备

输入数据为差异分析后得到的Rdata，我使用的是三大R包差异分析，随机选取了一个R包的分析结果作为输入数据

rm(list = ls())
options(stringsAsFactors = F)

library(clusterProfiler)
library(org.Hs.eg.db)

# 读取3个R包的差异分析结果
limma_DEG <- load(file = "limma_voom_DEG.Rdata")
DESeq2_DEG <- load(file = "DESeq2_DEG.Rdata")
edgeR_DEG <- load(file = "edgeR_DEG.Rdata")

# 根据需要更改DEG的值
DEG <- DEG3  # 使用limma产生的结果

 2.设置差异倍数阈值，用于划分上下调基因

# 设置log2倍变化的截止值
# 这里我们计算DEG（差异表达基因）数据框中logFC绝对值的平均数和标准差
# 然后使用平均数加上两倍的标准差作为截止值
logFC_cutoff <- with(DEG, mean( abs(logFC) ) + 2 * sd(abs(logFC)))

# 将截止值四舍五入到两位小数，使得结果更加可读
logFC_cutoff <- round(logFC_cutoff,2)

# 创建一个新列'change'来标记基因表达的变化趋势
# 我们遍历DEG数据框的每一行，并根据p值和log2倍变化值设置分类标签
# 如果p值小于0.05并且log2倍变化的绝对值大于logFC_cutoff，那么我们标记基因表达为显著变化
# 如果log2倍变化大于logFC_cutoff，我们标记为'UP'；反之，如果小于-logFC_cutoff，我们标记为'DOWN'
# 如果p值不小于0.05或者log2倍变化的绝对值不大于logFC_cutoff，我们标记为'STABLE'
DEG$change = as.factor( ifelse(DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff,
                                  ifelse( DEG$logFC > logFC_cutoff , 'UP', 'DOWN' ), 'STABLE' ) )
table(DEG$change)

3.从org.Hs.eg.db提取Ensembl ID和GI号对应关系

这一步是我遇到比较麻烦的报错的一步。首先差异分析得到的数据框DEG没有ensembl id这一列，需要先用bitr函数把基因符号转换为ensembl id；转换完之后得到一个带有基因符号和ensembl id的数据框，此时需要单独将ensembl id这一列提取出来作为一个向量，添加到DEG数据框中，用于下一步的富集分析，由于不是所有的基因符号都能成功转换，所以在向DEG数据框添加ensembl id这一列的时候还需要额外填充NA值；最后再将ensembl id转换为GI号，得到df数据框。

library(org.Hs.eg.db)

# 查看org.Hs.eg.db数据库中可用的键类型，这些类型可用于标识基因
keytypes(org.Hs.eg.db)

# 将未处理的行名赋值给DEG数据框的一个新列，命名为'gene_symbols'
DEG$gene_symbols <- rownames(DEG)

# 对gene_symbols列的基因标签进行处理
gene_symbols <- DEG$gene_symbol
# 使用bitr函数将基因符号转换为Ensembl基因ID
ensembl_ids <- bitr(gene_symbols, fromType = "SYMBOL", toType = "ENSEMBL", OrgDb = org.Hs.eg.db)

# 将Ensembl id从数据框中取出，作为一个向量
gene_list <- ensembl_ids$ENSEMBL


# 将Ensembl id添加到DEG数据框中
# 匹配基因符号到其在`ensembl_ids`中的位置
positions <- match(DEG$gene_symbols, ensembl_ids$SYMBOL)
# 创建与`DEG`行数相等的`ENSEMBL` ID向量
ensembl_ids_full <- rep(NA_character_, nrow(DEG))
# 填充找到的`ENSEMBL` ID
ensembl_ids_full[!is.na(positions)] <- ensembl_ids$ENSEMBL[positions[!is.na(positions)]]
# 将`ENSEMBL` ID向量添加到`DEG`数据框
DEG$ENSEMBL <- ensembl_ids_full
# 现在DEG数据框中应该有一个名为ENSEMBL的新列
head(DEG)


# 将Ensembl id转换为GI号
df <- bitr(gene_list, fromType = "ENSEMBL", toType = c( "ENTREZID" ), OrgDb = org.Hs.eg.db )
head(df)

4.将GI号合并到DEG数据框里

# 定义文件路径
f <- 'DEG_ENTREZID_for_GSEA_GSVA.Rda'

# 检查文件是否存在
if (!file.exists(f)) {  
  # 如果文件不存在
  # 将DEG数据框的行名（假设是基因符号）赋值给新列'SYMBOL'
  DEG$SYMBOL <- rownames(DEG)
  
  # 使用'merge'函数将DEG数据框与df数据框通过'ENSEMBL'列进行合并
  # 结果将保存在DEG数据框中，'by'参数指定了用于匹配的列名
  DEG <- merge(DEG, df, by = 'ENSEMBL')
  
  # 显示合并后的DEG数据框的前几行数据
  head(DEG)
  
  # 输出合并后的DEG数据框的维度（行数和列数）
  dim(DEG)
  
  # 将合并后的DEG数据框保存到磁盘，以便后续使用
  save(DEG, file = "DEG_ENTREZID_for_GSEA_GSVA.Rda")
}

# 如果文件存在，直接从磁盘加载DEG数据框
load("DEG_ENTREZID_for_GSEA_GSVA.Rda")

5.提取上下调的GI集

# 分别提取上下调基因的ENTREZID
gene_up <- DEG[ DEG$change == 'UP', 'ENTREZID' ] 
gene_down <- DEG[ DEG$change == 'DOWN', 'ENTREZID' ]

# 合并上下调基因的ENTREZID
gene_diff <- c( gene_up, gene_down )

# 获取所有基因的ENTREZID
gene_all <- as.character(DEG[ ,'ENTREZID'] )

# 将logFC值与对应的ENTREZID关联起来
geneList <- DEG$logFC
names( geneList ) <- DEG$ENTREZID

# 对logFC值进行降序排序
geneList <- sort( geneList, decreasing = T )   #从大到小排序

# KEGG通路富集分析
f  <- 'enrichKEGG.Rdata'

# 如果文件不存在
if(!file.exists(f)){
  # 对上调基因进行KEGG通路富集分析
  kk.up <- enrichKEGG(  gene          =  gene_up  ,  # 输入基因列表
                        organism      =  'hsa'    ,  # 生物物种，这里是人（Homo sapiens）
                        universe      =  gene_all ,  # 所有可能的基因背景

                        pvalueCutoff  =  0.8      ,  # 设置p-value阈值，此处可能设置过高，一般设置更低如0.05
                        qvalueCutoff  =  0.8      )  # 设置FDR校正后的q-value阈值，同样可能设置过高
  # 对下调基因进行KEGG通路富集分析
  kk.down <- enrichKEGG(gene          =  gene_down,
                        organism      =  'hsa'    ,
                        universe      =  gene_all ,
                        pvalueCutoff  =  0.05     )  # 这里的p-value阈值设置更为合理
  # 保存分析结果
  save(kk.up,kk.down,file = "enrichKEGG.Rdata")
}
# 加载分析结果
load(file = "enrichKEGG.Rdata")
head(kk.up)[ ,1:6 ]
head(kk.down)[ ,1:6 ]

6.ggplot2画图

library(ggplot2)

kegg_down_dt <- as.data.frame( kk.down )
kegg_up_dt <- as.data.frame( kk.up )
down_kegg <- kegg_down_dt[ kegg_down_dt$pvalue < 0.05, ]
down_kegg$group = -1
up_kegg <- kegg_up_dt[ kegg_up_dt$pvalue < 0.05, ]
up_kegg$group = 1
  
dat = rbind( up_kegg, down_kegg )
dat$pvalue = -log10( dat$pvalue )
dat$pvalue = dat$pvalue * dat$group
  
dat = dat[ order( dat$pvalue, decreasing = F ), ]
  
g_kegg <- ggplot( dat, aes(x = reorder( Description, order( pvalue, decreasing=F ) ), y = pvalue, fill = group)) + 
  geom_bar( stat = "identity" ) + 
  scale_fill_gradient( low = "blue", high = "red", guide = FALSE ) + 
  scale_x_discrete( name = "Pathway names" ) +
  scale_y_continuous( name = "log10P-value" ) +
  coord_flip() + theme_bw() + theme( plot.title = element_text( hjust = 0.5 ) ) +
  ggtitle( "Pathway Enrichment" ) 
print( g_kegg )
ggsave( g_kegg, width = 10, height = 6, filename = 'kegg_up_down_anno.png' )