28a05a492e5be075caa75cfbcef6d1e1.gif

在GSEA分析中,在MSigDB(Molecular Signatures Database)数据库中定义了很多基因集,下载的基因集是gmt格式文件。下载的gmt格式文件,打开后可以看见是下面这个样子的:

f2898697e33c8c3846d13923960f82a4.png

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

ad3df8b9c1973000e9c38961451bc877.png

72c9b6266847f2a34bdede64fb77eb91.png

ccd8d7ff1d4e0433756d695436af8d1f.png

fe6aa3e2ecf956f2c4541bdb28a2a621.png

293bab34382db6170bb5a054a43a165c.png

929b72284073a82ddaa1b60170f5e08b.png

9123bdb1d0541cca7716d8d08ad4f8a8.png

a0001ffc3b797a5aed03d87c9ec27fa0.png

8a870544a5331720894f6b0eaa6fd7e5.png

Logo

有“AI”的1024 = 2048,欢迎大家加入2048 AI社区

更多推荐