Skip to content

R (프로그래밍 언어) Bioconductor GOstats #
Find similar titles

Structured data

Category
Software

개요 #

Bioconductor에는 많은 생물정보학 관련 R package가 있다. 그 중에서 GO 분석과 관련된 Package인 GO stats를 소개한다. 많은 기능을 가지고 있지만 기본적으로 원하는 Entrez Gene IDs가 갖는 GO term을 그래프로 표현하는 기능을 가지고 있다.

해석 #

GO의 경우 유전자들의 functional analysis에 널리 사용되고 있다 하지만 GO의 경우 그 Level의 Depth가 깊고 그 관계가 일방적이지 않아 Text를 가지고는 확인 하기 어려운 의미가 많다. 따라서 그래프로 표현을 하여 어떠한 Term이 있고 그 관계를 확인 하는 것은 분석에 많은 영향을 줄 수 있다.

Codes #

아래 코드는 연습의 용이함을 위해 Bioconductor에서 제공되는 Manual의 예제를 사용하였다

### R code from vignette source 'GOvis.Rnw'

###################################################
### code chunk number 1: Setup
###################################################
library("Biobase")
library("annotate")
library("GOstats")
library("xtable")
library("multtest")
library("Rgraphviz")
library("hgu95av2.db")
library("GO.db")
library("genefilter")
library("ALL")


###################################################
### code chunk number 2: Example
###################################################
## subset of interest: 37+42 samples
data(ALL)
eset <- ALL[, intersect(grep("^B", as.character(ALL$BT)),
                        which(as.character(ALL$mol) %in% c("BCR/ABL","NEG", "ALL1/AF4")))]

MB = factor(eset$mol.bio)

## intensities above 100 in at least 7 of the samples
f1 <- kOverA(7, log2(100))
f2 = function(x) {gps = split(2^x, MB); mns=sapply(gps, mean);
                  if(max(mns) - min(mns) < 100) FALSE else TRUE}
ff <- filterfun(f1, f2)
selected <- genefilter(eset, ff)
sum(selected)
esetSub <- eset[selected,]



###################################################
### code chunk number 3: mtFstats
###################################################

pvCutOff = 0.05

Fstats = mt.maxT(exprs(esetSub), as.numeric(MB)-1,
                 B=1000, test="f")

eSet = esetSub[Fstats$index[Fstats$adjp<pvCutOff],]

gN = featureNames(eSet)
lls = unlist(mget(gN, hgu95av2ENTREZID, ifnotfound=NA))

eS = eSet[!duplicated(lls),]
gN = featureNames(eS)
lls = unlist(mget(gN, hgu95av2ENTREZID, ifnotfound=NA))
lls = as.character(lls[!is.na(lls)])
syms = unlist(mget(gN, hgu95av2SYMBOL, ifnotfound=NA))



###################################################
### code chunk number 4: inducedGO
###################################################
gMF = makeGOGraph(lls, "MF", chip="hgu95av2")
gBP = makeGOGraph(lls, "BP", chip="hgu95av2")
gCC = makeGOGraph(lls, "CC", chip="hgu95av2")


nMF = list()
lbs = rep("", length(nodes(gMF)))
names(lbs) = nodes(gMF)
nMF$label = lbs

nBP = list()
lbs = rep("", length(nodes(gBP)))
names(lbs) = nodes(gBP)
nBP$label = lbs

nCC = list()
lbs = rep("", length(nodes(gCC)))
names(lbs) = nodes(gCC)
nCC$label = lbs

if( require("Rgraphviz") ) {
  gMFlo = agopen(gMF, "gMF")
  gBPlo = agopen(gBP, "gBP")
  gCClo = agopen(gCC, "gCC")
}

plot(gMFlo, nodeAttrs=nMF, main="Molecular Function")
plot(gBPlo, nodeAttrs=nBP, main="Biological Process")
plot(gCClo, nodeAttrs=nCC, main="Cellular Component")
  • 결과

Image

Image

Image

Image

Suggested Pages #

0.0.1_20140628_0