GSEA分析中的gmt格式文件如何自定义

DoubleHelix

发布于 2021-09-29 15:33:08

4.6K1

发布于 2021-09-29 15:33:08

在我前面的文章：clusterProfiler包进行KEGG,GO,GSEA富集分析，有介绍在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=gcSample
write.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 Essential
1     Proteins     A2RU49         HYKK Lysine catabolism R-HSA-71064.2       Yes
2     Proteins     Q9BQT8     SLC25A21 Lysine catabolism R-HSA-71064.2       Yes
3     Proteins     Q92947         GCDH Lysine catabolism R-HSA-71064.2       Yes
4     Proteins     Q8N5Z0        AADAT Lysine catabolism R-HSA-71064.2       Yes
5     Proteins     Q8IUZ5       PHYKPL Lysine catabolism R-HSA-71064.2       Yes
6     Proteins     Q9P0Z9        PIPOX Lysine catabolism R-HSA-71064.2       Yes

MoleculeName和 catabolism.Type这2列是我们要的。

可以自己构建类似上面gcSample的list，然后自己写一个函数输入就行。

name <- unique(geneset$catabolism.Type)
description <- rep(NA,length(name))
names(description) <- name
genes <- lapply(name, function(name){
   as.vector(geneset[geneset$catabolism.Type == name,"MoleculeName"])
})
names(genes) <- name

gmtinput <- 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)

参考：

上次说的gmt函数(学徒作业)

https://blog.csdn.net/coding_Joash/article/details/120422166