代码大部分来自简书:猪莎爱学习,部分修改过
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' )