个生物过程通常是由一组基因共同参与，而不是由单个基因独自完成，常见的富集分析有go，kegg，gsea，gsva等，我们在这里可以尝试进行测试
测试1：go分析
GO(gene ontology)是基因本体联合会(Gene Onotology Consortium)所建立的数据库，旨在建立一个适用于各种物种的、对基因和蛋白质功能进行限定和描述的、并能随着研究不断深入而更新的语言词汇标准。分别描述基因的分子功能（molecular function）、细胞组分（cellular component）、参与的生物过程（biological process）。首先将差异表达转录本向GO数据库（http://www.geneontology.org/）的各term映射，并计算每个term的转录本数，从而得到具有某个GO功能的转录本列表及转录本数目统计。然后应用超几何检验，找出与整个转录本组背景相比，在差异表达转录本中显著富集的GO条目,GO 分析对实验结果有提示的作用，通过差异基因的GO 分析，可以找到富集差异基因的GO分类条目，寻找不同样品的差异基因可能和哪些基因功能的改变有关。
1 加载我们需要用到的包
suppressMessages({
library(Seurat)
#library(corrplot)
library(tidyverse)
library(RColorBrewer)
library(clusterProfiler)
library(org.Hs.eg.db)
#library(org.Mm.eg.db)小鼠库
})
2 读取我们需要的数据,并筛选出差异基因
PRO<-readRDS('lung.diff_PRO.rds')
genes.use <- HVFInfo(object = PRO)
cluster.markers <- FindAllMarkers(object = PRO, genes.use = genes.use)
write.csv(cluster.markers,"lung_markers.csv")
cluster.marker<-read.csv('lung_markers.csv')
3 进行go富集分析，在这里我们可以使用一个cluster和所有cluster 的对比分析，也可以针对2个cluster进行富集分析
do_go <- function(genes){
converted <- select(org.Hs.eg.db,
keys = genes,
columns = c("ENTREZID", "SYMBOL"),
keytype = "SYMBOL")
gene_id <- na.omit(converted$ENTREZID)
ego <- enrichGO(gene = gene_id,
OrgDb = org.Hs.eg.db,
ont = "ALL",
pAdjustMethod = "BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05,
readable = TRUE)
return (ego)
#针对所有cluster
for (cluster in unique(cluster.marker$cluster)){
markers <- as.vector(Cancer_data.markers[cluster.marker$cluster==cluster & cluster.marker$avg_log2FC >0.25,]$gene)
tryCatch({
ego <- do_go(markers)
result_file <- paste("lung",cluster,"_go.csv",sep="")
write.csv(ego@result,result_file)
pdf(file = paste("lung",cluster,"_go_dot.pdf",sep=""),width =15,height = 12)
print (dotplot(ego,showCategory=20))
dev.off()
pdf(file = paste("lung",cluster,"_go_bar.pdf",sep=""),width =15,height = 12)
print (barplot(ego,showCategory=20))
dev.off()
},error=function(e){cat("There is an error:",conditionMessage(e),"\n")})
}
}
#只针对2个cluster
cluster.markers <- FindMarkers(object = PRO, ident.1 = 'CancerCells',ident.2 = 'AT2', logfc.threshold = 0.25)
cluster1<- row.names.data.frame(cluster.markers)
cluster1=bitr(cluster1,fromType = "SYMBOL",toType = c("ENTREZID"),OrgDb = "org.Hs.eg.db")
ego<-enrichGO(gene=cluster1[,"ENTREZID"],keyType = "ENTREZID",OrgDb=org.Hs.eg.db,ont = "ALL",pAdjustMethod = "BH",pvalueCutoff = 0.01,qvalueCutoff = 0.05,readable = TRUE)
write.csv(ego@result,'luncancer.csv')
pdf(file = 'lung_cancer_dot.pdf',width =15,height = 12)
print (dotplot(ego,showCategory=20,split="ONTOLOGY")+ facet_grid(ONTOLOGY~., scale="free"))
dev.off()
pdf(file = 'lung_cancer_bar.pdf',width =15,height = 12)
print (barplot(ego,showCategory=20))
dev.off()
#展示
p1<-dotplot(ego,showCategory=10,split="ONTOLOGY")+ facet_grid(ONTOLOGY~., scale="free")
p2<-barplot(ego,showCategory=10,split="ONTOLOGY")+ facet_grid(ONTOLOGY~., scale="free")
p<-p1+p2
p
image.png
测试2：KEGG：京都基因和基因组百科全书：KEGG项目于1995年5月在日本教育，科学，体育和文化部的人类基因组计划下启动。KEGG由三个数据库组成：用于表示相互作用分子网络的高阶函数的pathway，用于收集所有完全测序的基因组和一些部分基因组的基因目录的GENES，以及用于化学品收集的LIGAND
#针对单一cluster和所有cluster对比
do_kegg <- function(genes){
converted <- select(org.Hs.eg.db,
keys = genes,
columns = c("ENTREZID", "SYMBOL"),
keytype = "SYMBOL")
gene_id <- na.omit(converted$ENTREZID)
ekegg <- enrichKEGG(gene = genelist, organism = 'hsa')
return (ekegg)
}
for (cluster in unique(cluster.marker$cluster)){
markers <- as.vector(cluster.marker[cluster.marker$cluster==cluster & cluster.marker$avg_log2FC >0.25,]$gene)
tryCatch({
ekegg <- do_kegg(markers)
#result_file <- paste("lung",cluster,"_go.csv",sep="")
#write.csv(ego@result,result_file)
#pdf(file = paste("lung",cluster,"_go_dot.pdf",sep=""),width =15,height = 12)
#print (dotplot(ego,showCategory=20))
#dev.off()
pdf(file = paste("lung",cluster,"_kegg_bar.pdf",sep=""),width =15,height = 12)
print (barplot(ego,showCategory=20))
dev.off()
},error=function(e){cat("There is an error:",conditionMessage(e),"\n")})
}
#针对单一cluster和单一cluster对比
cluster.markers <- FindMarkers(object = PRO, ident.1 = 'CancerCells',ident.2 = 'AT2', logfc.threshold = 0.25)
genelist <- bitr(row.names(cluster.markers), fromType="SYMBOL",
toType="ENTREZID", OrgDb='org.Hs.eg.db')
genelist <- pull(genelist,ENTREZID)
ekegg <- enrichKEGG(gene = genelist, organism = 'hsa')
p1 <- barplot(ekegg, showCategory=20)
p2 <- dotplot(ekegg, showCategory=20)
plotc = p1/p2
plotc
#ggsave("enrich/enrichKEGG.png", plot = plotc, width = 12, height = 10)
image.png
测试3：GSEA
GSEA（Gene Set Enrichment Analysis）：基因集富集分析，由Broad Institute研究所提出的一种富集方法，同时还提供对应的分析软件GSEA和一个基因集数据库MSigdb，它的基本思想是使用预定义的基因集（通常来自功能注释或先前实验的结果），将基因按照在两类样本中的差异表达程度排序，然后检验预先设定的基因集合是否在这个排序表的顶端或者底端富集。基因集合富集分析检测基因集合而不是单个基因的表达变化，因此可以包含这些细微的表达变化，预期得到更为理想的结果。
1 加载我们需要的包
ibrary(fgsea)
library(enrichplot)
library(msigdbr)
2 我们继续使用上边使用的数据，先得的我们要对比的cluster 的差异基因
cluster.markers <- FindMarkers(object = PRO, ident.1 = 'CancerCells',ident.2 = 'AT2', logfc.threshold = 0.25)
3 把数据整理成我们需要的模式,并使用KEGG的数据库进行分析
nrDEG = cluster.markers[,c('avg_log2FC', 'p_val')]
colnames(nrDEG)=c('log2FoldChange','pvalue') ##更改列名
#library(org.Hs.eg.db)
#library(clusterProfiler)
## 把SYMBOL转换为ENTREZID，可能有部分丢失
gene <- bitr(rownames(nrDEG),fromType = "SYMBOL",toType = "ENTREZID",OrgDb = org.Hs.eg.db)
## 基因名、ENTREZID、logFC一一对应起来
gene$logFC <- nrDEG$log2FoldChange[match(gene$SYMBOL,rownames(nrDEG))]
## 构建genelist
geneList=gene$logFC
names(geneList)=gene$ENTREZID
geneList=sort(geneList,decreasing = T) # 降序，按照logFC的值来排序
## GSEA分析
kk_gse <- gseKEGG(geneList = geneList,
organism = 'hsa',
nPerm = 1000,
minGSSize = 10,
pvalueCutoff = 0.9,
verbose = FALSE)
kk_gse=DOSE::setReadable(kk_gse, OrgDb='org.Hs.eg.db',keyType='ENTREZID')
sortkk<-kk_gse[order(kk_gse$enrichmentScore, decreasing = T),]
gseaplot2(kk_gse,
"hsa04510",
color = "firebrick",
rel_heights=c(1, .2, .6))
#gseaplot2(kk_gse, row.names(sortkk)[1:4]) 也可以看一下前4条通路
image.png
4 使用msigdbr数据库进行分析
#选择自己用到的数据库的数据集
genesets <- msigdbr(species = "Homo sapiens", category = "C2") %>% dplyr::select("gs_name","gene_symbol" )%>% as.data.frame()
#https://www.gsea-msigdb.org/gsea/msigdb 不同数据集包含的内容，在这里我们选择c2数据集：（专家）校验基因集合，基于通路、文献等：
genesets <- split(genesets$gene_symbol, genesets$gs_name)
#整理自己的数据格式
cluster.markers$gene<-rownames(cluster.markers)
gsea_genes<-cluster.markers %>%arrange(avg_log2FC) %>%dplyr::select(gene,avg_log2FC)
ranks <- deframe(gsea_genes)
#分析与画图
fgseaRes <- fgsea(pathways = genesets,stats = ranks,minSize=5,maxSize=500,nperm=10000)
p<-ggplot(fgseaRes %>% as_tibble() %>% arrange(desc(NES)) %>% filter(pval < 0.05) %>% head(n= 20), aes(reorder(pathway, NES), NES)) +
geom_col(aes(fill= NES)) +
coord_flip() #+
#labs(x="KEGG", y="Normalized Enrichment Score",title="KEGG gene sets NES from GSEA") ##输出差异排秩前20的条目
pdf('fgsea.pdf',width=20,height=10)
print(p)
dev.off()
#可以查看一下自己感兴趣的通路
plotEnrichment(genesets[["RHODES_UNDIFFERENTIATED_CANCER"]],ranks) + labs(title="RHODES_UNDIFFERENTIATED_CANCER")
image.png
image.png
测试4：GSVA
基因富集分析（Gene set enrichment analysis，GSEA）需要预先有一个有生物学意义的基因集（如同一个通路中的基因），再将基因集内（具有相同意义/功能）的基因一并计算并归纳为单个富集分数。这个分析思路增加了可解释性，多用于评估基因集所在通路/功能的活性变化。基因集变异分析 (Gene Set Variation Analysis，GSVA)是其中一个GSEA的算法。也可以一开始从基因表达量和多个通路信息出发，根据通路活性变化情况，对样本进行无监督分类。
1 加载相关的R包
suppressMessages({
library(tidyverse)
library(GSVA)
library(limma)
library(GSEABase)
library(pheatmap)
library(Seurat)
library(msigdbr)
})
2 读入相关的数据,读入要使用的数据
PRO<-readRDS('lung.diff_PRO.rds')
expr <- AverageExpression(PRO, assays = "RNA", slot = "data")[[1]]
expr <- expr[rowSums(expr)>0,]#选取非零基因
expr <- as.matrix(expr)
3 读入使用的数据库
msgdC2 = msigdbr(species = "Homo sapiens", category = "C2",subcategory = "KEGG")
keggSet = msgdC2 %>% split(x = .$gene_symbol, f = .$gs_description)
4 进行GSVA分析
expr_gsca<- gsva(expr, gset.idx.list = keggSet, kcdf="Gaussian",method = "zscore",parallel.sz=1)
saveRDS(expr_gsca,'gsva.rds')
gsva.df <- data.frame(Genesets=rownames(expr_gsca), expr_gsca, check.names = F)
write.csv(gsva.df, "gsva_res.csv", row.names = F)
image.png
网友评论