KEGG, Gene-Ontology Geneset Object 만들기
Introduction
구글에 검색해서 5년 이상 지난 쓰레기만 뒤적뒤적대다가 해결한 내용.
Pathway analysis 에 있어서 가장 중요한것은 아마 Pathway를 어떻게 정의 하느냐 일 것이다.
가장 직관적인 예로는, 실험맨들이 우리가 이런 실험을 했는데 이런 Gene들이 express 되더라,
하는 내용을 바탕으로
Experiment1 -> geneA, geneB, geneC
Experiment2 -> geneD, geneE, geneF
…
와 같이 정해둔 경우가 있다.
이 외에 Gene의 물리적인 Position에 따라서나, 계산을 통해 이 Gene 들은 같이 상호작용 할 것이다.
와 같은 예시들도 있다.
이러한 정보를 Pathway Analysis에 사용 하기 위해서는 경험상 2가지 방법으로 데이터를 다루는데
- dataframe 방식
- list 방식
이다.
위의 예를 다시 사용하면
dataframe 방식은 |Geneset|Gene| |:—:|:—:| |Experiment1|geneA| |Experiment1|geneB| |Experiment1|geneC| |Experiment2|geneD| |Experiment2|geneE| |Experiment2|geneF|
와 같이 말 그대로 data.frame 처럼 (matrix던 뭐던 그냥 data.frame이라고 표기함)
저장하는 방법이고.
list 방식은
$Experiment1
geneA geneB geneC
$Experiment2
geneD geneE geneF
처럼 저장하는 방법이다.
나는 처음 접할때 list형태를 접했기 때문에, list방법을 매우 선호하고,
따라서 dataframe 형태로 저장되어 있는 것도 list형태로 바꿔서 사용하는 편이다.
Method - Human
Human의 경우는 사실 Geneset을 구하기가 매우 편리하다. 왜냐면 연구를 하는 사람이 많으니까,
개인적으로 제일 애용하는 하는 곳은 MSigDB.
GSEA로 유명한 broad에서 친히 geneset을 정리해두었다. c2(Curated, KEGG) 나, c5(GO)를 바로 다운로드 받아서
GSA 라이브러리를 사용하여 로드 후 아래처럼 저장하여 필요할때마다 로드해서 쓰면 된다.
library(GSA)
tmp <- GSA.read.gmt('c2.all.v7.0.symbols.gmt')
genesets <- tmp$genesets
names(genesets) <- tmp$geneset.names
save(genesets, file = 'c2v7.RData')
Method - mouse의 경우
이제 문제가 달라진다, Human은 앞에서도 말했던것 처럼 수요가 워낙 많아서 온갖 Geneset들이 널려있다. 그러나 만약 Human이 아니면? 특히 나처럼 Experimental Background가 전혀 없다면 이러한 geneset을 어디서 구해야 할까?
내가 찾은 한가지 방법은 GSKB이다.
South dakota bioinformatics group에서 만들어둔 geneset인데 꽤 많은 species에 대하여 정리 해두었고
gmt 포맷도 지원하니 위의 human처럼 gsa를 사용해서 읽을 수도 있다.
문제는 우리 교수님은 이걸 인정하지 않으신다는거지.
많은 사람들이 나름 “공신력” 있는 geneset으로 KEGG 와 GO를 선호하기 때문에,
mouse (mus musculus) 에 대해서 KEGG geneset과 GO geneset을 만드는 방법을 정리한다.
우선 작년에 친애하는 소라짱이 ppt로 정리한 biomaRt를 이용하여 ensembl을 통해 geneset 을 만드는 방법을 참조 하라고 받았다.
if (!requireNamespace(“BiocManager”, quietly = TRUE))
install.packages(“BiocManager”)
BiocManager::install(“biomaRt”) ## Install package
library(biomaRt) ## Load package
listMarts() ## Check the list of biomarts
ensembl = useMart(“ENSEMBL_MART_ENSEMBL”) ## Select the marts
dataset = listDatasets(ensembl) ## List of species
head(dataset)
ensembl = useDataset(“hsapiens_gene_ensembl”, mart=ensembl) ## Select the species
filters = listFilters(ensembl) ## List of filters
attributes = listAttributes(ensembl) ## List of attributes (features)
attributes[1:5,]
attributes[grep(‘GO’, attributes$description),] ## Search for the name of attributes related to GO
attributes[grep(‘KEGG’, attributes$description),] ## Search for the name of attributes related to KEGG pathway
## Get the data matrix
data = getBM(attributes = c(“go_id”, “name_1006”, “namespace_1003",“kegg_enzyme”, “external_gene_name”,“ensembl_gene_id”), mart = ensembl)
data ## See it
## Extract genes involved in MAP kinase activity (GO:0004707)using ‘go’ filter
getBM(attributes = c('external_gene_name'), filters = 'go', values = 'GO:0004707', mart = ensembl)
근데 이건 Human에 대한거였고, Mouse쪽으로 바꿔 사용 함을 감안해도 GO쪽은 유용했으나,
KEGG는 kegg pathway 이름을 주지 않고 enzyme이랑 섞어서 알 수 없는 숫자를 줬기 때문에 전혀 사용 할 수가 없었다.
library(biomaRt) # 라이브러리 로드
ensembl <- useMart("ENSEMBL_MART_ENSEMBL") # 마트 설정
dataset <- listDatasets(ensembl) # ensembl로 데이터셋 설정
dataset[grep("Mouse", dataset$description), ] # mouse 데이터
ensembl <- useDataset("mmusculus_gene_ensembl", mart = ensembl)
attributes <- listAttributes(ensembl) # attributes 설정
attributes[grep("GO", attributes$description), ] # GO 검색
dataGO <- getBM(attributes = c("external_gene_name", "name_1006"), mart = ensembl) # Gene - GO name 쌍의 dataframe형태 pathway
dataGO <- dataGO[order(dataGO$name_1006), ] # GO name 기준으로 정렬
dataGO <- dataGO[which(dataGO$name_1006 != ""), ] # GO name 비어있는 거 삭제
dataGO$external_gene_name <- toupper(dataGO$external_gene_name) #
dataGO$name_1006 <- toupper(dataGO$name_1006) # Gene과 GO name 전부 대문자 처리함. ( 필수 아님, 취향 )
# list 형태로 변환
genesets <- list()
terms <- unique(dataGO$name_1006)
for (i in 1:length(terms)) {
genesets[[i]] <- dataGO$external_gene_name[which(dataGO$name_1006 == terms[i])]
}
names(genesets) <- terms
lgs <- sapply(1:length(genesets), function(i) {
length(genesets[[i]])
})
genesets <- genesets[intersect(which(lgs >= 15), which(lgs <= 500))] # geneset 중에서 gene이 15 ~ 500개 가 아닌것 제외
save(genesets, file = "mouseGO.RData") # 저장
mouse KEGG geneset
google에 그대로 치면 나오는게 bioconductor의 pathview, CHRONOS 등이 나온다.
해봤는데 전부 구리다, 내가 원하는 걸 찾을 수 없었다. 쓰지 말자
대신 KEGGREST를 쓰면 해결 된다.
BiocManager::install("KEGGREST") # package install
library(KEGGREST)
k <- names(keggList("pathway", "mmu"))
kv <- sapply(unname(keggList("pathway", "mmu")), function(i) {
strsplit(i, " - ")[[1]][1]
}, USE.NAMES = FALSE)
k <- lapply(k, function(i) {
g <- keggGet(i)[[1]]$GENE
g <- g[which(1:length(g) %% 2 == 0)]
g <- sapply(g, function(j) {
strsplit(j, ";")[[1]][1]
}, USE.NAMES = FALSE)
return(unlist(g))
})
names(k) <- kv # geneset 이름 설정
genesets <- k # geneset으로 변경
save(genesets, file='mouse_kegg.Rdata") 저장
```R
------
번외, <br>
예전에는
```R
source('http://bioconductor.org/biocLite.R') # 하도 쳐서 poweroverwhelming처럼 손이 코드를 외움
biocLite(~~)
로 bioconductor 패키지를 설치 했던 것 같은데 이제 나도 늙었나보다.
- github blog 에는 detail summary를
Comments