电脑建立网站平台,福州seo推广公司,成都市建设局网站,建设网站收费使用KEGG通路的基因列表进行单细胞GSEA GSVA分析的过程#xff0c;我们需要遵循以下步骤#xff1a; 获取KEGG通路的基因列表#xff1a;这通常涉及使用专门的R包#xff0c;如KEGGREST或biomaRt#xff0c;来查询KEGG数据库并检索特定通路的基因列表。 准备单细胞表达数…使用KEGG通路的基因列表进行单细胞GSEA GSVA分析的过程我们需要遵循以下步骤 获取KEGG通路的基因列表这通常涉及使用专门的R包如KEGGREST或biomaRt来查询KEGG数据库并检索特定通路的基因列表。 准备单细胞表达数据这包括加载单细胞RNA-seq数据通常使用Seurat或其他单细胞分析包进行预处理。 执行GSVA分析使用GSVA包对单细胞数据执行基因集变异分析GSVA根据KEGG通路的基因列表评估每个单细胞样本的通路活性。 可视化GSVA结果最后基于GSVA分析结果绘制热图或其他类型的图表来展示不同单细胞样本中通路活性的变化。 今天我们主要关注第一步如何获取KEGG通路的基因列表
问题来源
不管是转录组数据还是单细胞数据都可以做gsva分析。gsva需要两个文件作为输入
1. 表达矩阵
2. 基因集 表达矩阵容易获得但是如果我们想做kegg数据库的通路分析怎么办如何获取kegg的通路列表
获取kegg的通路列表代码 方法一使用msigdb library(msigdbr) genesets msigdbr(species Homo sapiens ) #msigdbr提供多个物种的基因集数据 # View(msigdbr_collections()) #查看msigdbr包中所有的基因集 unique(genesets$gs_subcat) # 有多个数据库来源的基因集可选这里选用KEGG genesets - subset(genesets, gs_subcatCP:KEGG, select c(gs_name, gene_symbol)) unique(genesets$gs_name) #查看有多少条通路186个
但是这里面的kegg只有186个基因集合 方法二 使用 KEGGREST #BiocManager::install(KEGGREST) #BiocManager::install(EnrichmentBrowser) library(KEGGREST) library(EnrichmentBrowser) KEGGREST:: listDatabases() KEGGREST::keggList(database kegg) keggList(organism) ## returns the list of KEGG organisms with #step2: check and obtain a list of entry identifiers (in this case: sar) and associated definition for a given database or a given set of database entries. res - keggList(pathway, hsa) ## returns the list of human pathways length(res) resas.data.frame(res) head(res) #step 3: download the pathways of that organism: hsapathway - downloadPathways(hsa) head(hsapathway) idTypes(org hsa ) #step 4: retrieve gene sets for an organism from databases such as GO and KEGG: hsa_kegg_genesets - getGenesets(org hsa, db kegg, gene.id.type SYMBOL, cache TRUE, return.typelist) #step5: Parse and write the gene sets to a flat text file in GMT format for other pathway enrichment analysis programs (e.g., GSEA): writeGMT(hsa_kegg_genesets, gmt.file kegg_hsa_kegg_genesets_gmt) save(hsa_kegg_genesets,file ~/heart_muscle/hsa_kegg_genesets.rds) 我们可以看到这里的kegg数据集合有357个 做gsva分析 library(GSVA);print(getwd()) load(~/heart_muscle/hsa_kegg_genesets.rds) genesets_kegghsa_kegg_genesets print(length(hsa_kegg_genesets)) #kegg--- gssea.res - gsva(expr, genesets_kegg [1:50], methodssgsea,kcdfPoisson,min.sz 3,max.sz200 ,parallel.sz10 ) saveRDS(gssea.res, paste0(new_dir,file_name,_gssea.res_kegg.rds ) ) gssea.df - data.frame(Genesetsrownames(gssea.res),gssea.res, check.names F) write.csv(gssea.df, paste0(new_dir,file_name,gssea_res_kegg.csv ), row.names F) # ssgsea----- #library(GSVA);print(getwd()) #gssea.res - gsva(expr, genesets_GO [1:50], methodssgsea,kcdfPoisson,min.sz 3,max.sz200 ,parallel.sz10 ) #,parallel.sz10 #saveRDS(gssea.res, paste0(new_dir,file_name,_gssea.res_go.rds ) ) #gssea.df - data.frame(Genesetsrownames(gssea.res), gssea.res, check.names F) #write.csv(gssea.df, paste0(new_dir,file_name,gssea_res_go.csv), row.names F) #print(done-------);print(getwd())
这里的expr是行为基因列为样本的 表达矩阵。
参考https://biobeat.wordpress.com/category/r/https://www.researchgate.net/post/How_i_can_get_a_list_of_KEGG_pathways_and_its_list_of_genes 分析完成之后可以使用新得到的通路矩阵进行差异分析。下期见~ 看完记得顺手点个“在看”哦
