生信笔记 | 自定义GSEA分析中的gmt格式文件
在GSEA分析中，在MSigDB（Molecular Signatures Database）数据库中定义了很多基因集，下载的基因集是gmt格式文件。下载的gmt格式文件，打开后可以看见是下面这个样子的：
gmt（Gene Matrix Transposed，基因矩阵转置）是多列注释文件，列与列之间都是Tab制表符分割。
第1列：是基因所属基因集的名字，可以是通路名字，也可以是自己定义的任何名字。
第2列 ：一般是描述信息，说明这套基因列表从哪里收集的，也可以为空或者用NA表示。官方提供的格式是URL，也可以是任意字符串。
第3列-第n列：是基因集内所有基因的名字，有几个写几列。
每一行的列数可以不一样，主要是基因集内的基因数量不一样。
gmt文件可用 read.gmt()函数读入，读入的数据是一个数据框。
gmt <- read.gmt("./c5.go.cc.v7.2.symbols.gmt")class(gmt)如何制作自定义的gmt文件？下面是来自生信技能树的案例代码：
library(clusterProfiler)data(gcSample)names(gcSample)file="sink-examp.txt"gs=gcSamplewrite.gmt <- function(gs,file){ sink(file) lapply(names(gs), function(i){ cat( paste(c(i,'tmp',gs[[i]]),collapse='\t') ) cat('\n') }) sink()}write.gmt(gs,file)gcSample数据是来自clusterProfiler包，只是用来练习，自己自定义的可能并不是这样的list，可以处理成类似gcSample数据的list，用上面代码写出gmt文件。
下面是我处理的一个基因集
geneset <- read.table("data/AAsMet.txt",header = T,sep = "\t")head(geneset)> head(geneset) MoleculeType Identifier MoleculeName catabolism.Type ID Essential1 Proteins A2RU49 HYKK Lysine catabolism R-HSA-71064.2 Yes2 Proteins Q9BQT8 SLC25A21 Lysine catabolism R-HSA-71064.2 Yes3 Proteins Q92947 GCDH Lysine catabolism R-HSA-71064.2 Yes4 Proteins Q8N5Z0 AADAT Lysine catabolism R-HSA-71064.2 Yes5 Proteins Q8IUZ5 PHYKPL Lysine catabolism R-HSA-71064.2 Yes6 Proteins Q9P0Z9 PIPOX Lysine catabolism R-HSA-71064.2 YesMoleculeName和 catabolism.Type这2列是我们要的。
可以自己构建类似上面gcSample的list，然后自己写一个函数输入就行。
name <- unique(geneset$catabolism.Type)description <- rep(NA,length(name))names(description) <- namegenes <- lapply(name, function(name){ as.vector(geneset[geneset$catabolism.Type == name,"MoleculeName"])})names(genes) <- namegmtinput <- list(name=name,description=description,genes=genes)get_gmt <- function(gmtinput,filename){ output <- file(filename, open="wt") lapply(gmtinput[["name"]],function(name){ outlines = paste0(c(name, gmtinput[["description"]][[name]], gmtinput[["genes"]][[name]]),collapse='\t') writeLines(outlines, con=output) }) close(output)}get_gmt(gmtinput=gmtinput,filename="data/catabolism.gmt")我自己定义了一个输入对象gmtInfo
setClass("gmtInfo",slots=list(name="vector",description="vector",genes ="list"))gmtInfo <- new("gmtInfo",name=name,description=description,genes = genes)定义用来处理gmtInfo对象的函数：
write.gmt1 <- function(filename,gmtInfo){ if(class(gmtInfo) == "gmtInfo"){ output <- file(filename, open="wt") lapply(gmtInfo@name,function(name){ writeLines(paste(c(name, gmtInfo@description[[name]],gmtInfo@genes[[name]]),collapse='\t'), con=output) }) close(output) }}write.gmt1(filename="data/catabolism1.gmt",gmtInfo = gmtInfo)
write.gmt2 <- function(filename,gmtInfo){ if(class(gmtInfo) == "gmtInfo"){ sink(filename) lapply(gmtInfo@name, function(name){ cat(paste(c(name, gmtInfo@description[[name]],gmtInfo@genes[[name]]),collapse='\t')) cat('\n') }) sink() }}write.gmt2(filename="data/catabolism2.gmt",gmtInfo = gmtInfo)参考：
https://blog.csdn.net/coding_Joash/article/details/120422166
|
|
|
|
|
|
|
|
|
|
|
|
赞 (0)