유전자 정보로 암을 예측하는 법: Classifying Microarray Sample
암은 유전적인 이유가 큰 질병이다. 만약 암 환자의 증상과 유전자의 관계를 찾을 수 있다면 특정 유전자를 가진 사람의 암 발생 여부를 예측할 수 있을 것이다.
그러나 문제가 단순치 않다. 유전자 수가 너무 많다는 것(조사 가능한 암 환자대비) 때문이다. 암 환자수를 데이터수, 유전자를 설명변수라 한다면 많은 설명변수들 간에는 분명 공분산을 가지는 설명변수가 존재할 확률이 높을 것이다. 다시 말해 공분산을 최소화하기 위해 설명변수를 제거하는 수고 또한 필요하다.
DMWR의 마지막 케이스인 “Microarray 샘플 분류하기”는 이런 문제점을 포함해 다음과 같이 요약된다.
- 데이터보다 설명변수의 수가 많다. (128개 vs 12,625개)
- 지도 분류(Supervised Classification) 문제이다.
- 반응변수가 2개 이상의 명목값을 가지는 Multi-class 분류 문제이다. (총 6개)
- Data는 행이 설명변수이고 열은 데이터(sample)로 구성되어 있어 유의해야 한다.
- 성(sex)과 나이의 관계처럼 공분산을 가지는 변수가 존재할 확률이 높다.
Data에 대하여
본 예측 분석의 목적은 특정 유전자를 가질 경우(설명변수), 유전자 변형의 유형(반응변수)을 예측하는 것이다. 이 유전자 변형은 암세포를 의미한다.
최초 데이터는 급성 림프구성 백혈병 진단을 받은 128명의 환자의 자료이다. 환자가 가진 종양의 종류(반응변수)는 6개인데 크게 T-cell ALL(33명)와 B-cell ALL(95명)로 구분된다.
table(ALL$BT)##
## B B1 B2 B3 B4 T T1 T2 T3 T4
## 5 19 36 23 12 5 1 15 10 2table(ALL$BT,ALL$mol.bio)##
## ALL1/AF4 BCR/ABL E2A/PBX1 NEG NUP-98 p15/p16
## B 0 2 1 2 0 0
## B1 10 1 0 8 0 0
## B2 0 19 0 16 0 1
## B3 0 8 1 14 0 0
## B4 0 7 3 2 0 0
## T 0 0 0 5 0 0
## T1 0 0 0 1 0 0
## T2 0 0 0 15 0 0
## T3 0 0 0 9 1 0
## T4 0 0 0 2 0 0본 케이스는 B-cell에 중점을 두고자하니 T-cell 데이터를 제거한다. 개인적으로는 T-cell도 NEG가 발생하는 유전자를 예측할 수 있을 것으로 보인다. 우선 빈도가 작은 NUP-98과 p15/p16를 제거하다. 새로운 데이터를 저장한다.
tgt.cases <- which(ALL$BT %in% levels(ALL$BT)[1:5] &
ALL$mol.bio %in% levels(ALL$mol.bio)[1:4])
ALLb <- ALL[,tgt.cases]ALLb$BT <- factor(ALLb$BT)
ALLb$mol.bio <- factor(ALLb$mol.bio)
save(ALLb, file = "myALL.Rdata")자. 이제 새로운 데이터는 94개의 데이터를 가진다.
es <- exprs(ALLb)
dim(es)## [1] 12625 94각 설명변수에 대한 값은 “측정된 express level”로 되어 있다. 다시 말해 모든 값은 동일 단위이다. 전체 데이터를 한 눈에 보면 다음과 같다.
summary(as.vector(es))## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.985 4.122 5.469 5.624 6.829 14.040library(genefilter)
library(ggplot2)어떻게 변수를 제거할까?
변수 선택 방법 중 하나인 wrapper는 시간이 많이 걸리는 단점을 가진다. 따라서 데이터가 많은 본 문제의 경우, wrapper대신 filter 방법을 사용한다.
각 데이터의 중앙값과 IQR을 살펴보자. 본 예제는 Biobase패키지에서 제공하는 여러 함수를 사용하는데 rowQ나 rowMedians도 이에 해당한다. apply 함수 대신 굳이 또 다른 패키지내 함수를 사용하는 이유는 계산 시간이 엄청 빠르기 때문이다. 실제 system.time()으로 median의 계산 시간을 비교한 결과 전자는 0.08초, 후자는 1.46초 걸린다.
rowIQRs <- function(em)
rowQ(em,ceiling(0.75*ncol(em))) - rowQ(em,floor(0.25*ncol(em)))
library(ggplot2)
dg <- data.frame(rowMed=rowMedians(es), rowIQR=rowIQRs(es))
ggplot(dg,aes(x=rowMed, y=rowIQR)) + geom_point() +
xlab("Median expression level") + ylab("IQR expression level") +
ggtitle("Main Characteristics of Genes Expression Levels")
IQR은 75%에서 25%의 값을 뺀 값을 의미하는데 위 그래프에서는 IQR이 0에 가까운 데이터가 많이 관측된다. 만약 특정 유전자의 특성값이 환자와 관계없이 변동이 거의 없다면 그 유전자(설명변수)는 결과 값(B-Type의 변형)을 구별하는데 유용하지 못한 변수라 할 수 있다.
그러나 쉽게 단정치 못하는 문제가 하나 있다. 설명변수간 상호작용 때문이다. 어떤 유전자의 변동성이 낮은 이유가 다른 유전자와 함께 있을 경우에 더욱 두드러질 수 있다는 말이다.
이처럼 설명변수의 특성을 파악하기 위해서는 각 설명변수를 독립적으로 볼게 아니라 변수간 의존성을 고려한 각 변수의 중요도를 계량화해야한다. RELIEF 방법은 그것을 가능하게 한다.(DMWR 84p) 그러나 이 문제에선 RELIEF 법보다는 다른 방법을 통해 변수 제거를 시도할 것이다.
변수 선택법 1 : 값의 변동성이 극도로 낮은 변수 제거
상호작용을 고려함에도 불구하고 만약 IQR이 극도로 작다면 그 설명변수는 변수간 상호보완을 고려할 필요 없이 의미없는 변수라고 할 수 있다. 보수적인 기준을 적용키 위해 모든 데이터의 IQR을 기준으로 이보다 1/5 작은 변수는 제거하도록 하자. 12,625의 유전자가 3,942으로 줄었다. 새로운 데이터를 다시 es에 저장한다.
ALLb <- resFilter$eset
es <- exprs(ALLb)
dim(es)## [1] 3942 94변수 선택법 2: ANOVA통한 변수 제거
만약 변수간에 평균값 차이가 서로 유의미하지 않다면(귀무가설 채택) 두 변수중 하나는 제거해도 된다. ANOVA분석은 2개 이상의 변수의 평균 값이 유의미한 가를 볼 수 있는 Tool이다. 같은 반응 변수일 때 p값 0.01을 기준으로 그 보다 큰 변수를 제거한다. ANOVA결과 746개의 유전자가 남게 된다.
f <- Anova(ALLb$mol.bio, p = 0.01)
ff <- filterfun(f)
selGenes <- genefilter(exprs(ALLb), ff)
ALLb <- ALLb[selGenes, ]
ALLb## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 746 features, 94 samples
## element names: exprs
## protocolData: none
## phenoData
## sampleNames: 01005 01010 ... LAL5 (94 total)
## varLabels: cod diagnosis ... mol.bio (22 total)
## varMetadata: labelDescription
## featureData: none
## experimentData: use 'experimentData(object)'
## pubMedIds: 14684422 16243790
## Annotation: hgu95av2es <- exprs(ALLb)
dim(es)## [1] 746 94상기 변수 선택법 1, 2를 사용할 수 있는 이유는 모든 설명변수내 값이 동일 단위(무차원)이기 때문이다. 이전에 풀었던 문제와는 이 점에서 차이가 있음을 상기하자. 결과에 대한 IQR vs Median 그래프는 다음과 같다.
dg <- data.frame(rowMed=rowMedians(es), rowIQR=rowIQRs(es))
ggplot(dg,aes(x=rowMed, y=rowIQR)) + geom_point() +
xlab("Median expression level") + ylab("IQR expression level") +
ggtitle("Distribution Properties of the Selected Genes")
상기 그래프는 각각의 변수에 따른 샘플이 서로 다른 스케일을 가지고 있음을 보여준다. Modeling 적용시 군집분석을 적용할 예정인데 군집분석의 경우 단위는 같으나 스케일이 크게 차이 난다면 각 샘플간 거리에 있어 평균 값이 큰 변수의 영향력이 커지게 된다. 따라서 표준화가 필요하다. 이 절차는 뒤에 Modleing Stage에서 하기로 한다.
변수 선택법 3: Random Forest의 변수 중요도 사용
Random Forest를 이용해 746개 설명변수의 중요도의 순위를 찾아 보자. 기준은 각각의 설명변수 제거시 분류 정확도의 평균 감소량이 가장 적은 순으로 상위 30개를 찾는다.
dt <- data.frame(t(es), Mut = ALLb$mol.bio)
dt$Mut <- droplevels(dt$Mut)
set.seed(1234)
rf <- randomForest(Mut ~ ., dt, importance = TRUE)
imp <- importance(rf)
rf.genes <- rownames(imp)[order(imp[,"MeanDecreaseAccuracy"],
decreasing = TRUE)[1:30]]상위 30개 변수와 각각의 반응변수에 대응하는 값의 중앙값은 구하면 아래와 같다.
sapply(rf.genes, function(g) tapply(dt[, g], dt$Mut, median))## X1635_at X1674_at X37015_at X1467_at X40504_at X40202_at
## ALL1/AF4 7.302814 3.745752 3.752649 3.708985 3.218079 8.550639
## BCR/ABL 8.693082 5.833510 4.857105 4.239306 4.924310 9.767293
## E2A/PBX1 7.562676 3.808258 6.579530 3.411696 3.455316 7.414635
## NEG 7.324691 4.244791 3.765741 3.515020 3.541651 7.655605
## X34210_at X39631_at X41470_at X40795_at X32116_at X36638_at
## ALL1/AF4 5.641130 4.412177 9.616743 3.867134 7.115400 9.811828
## BCR/ABL 9.204237 5.021266 5.205797 4.544239 7.966959 8.486446
## E2A/PBX1 8.198781 4.901999 3.931191 4.151637 7.359097 6.259730
## NEG 8.791774 4.537882 4.157748 3.909532 7.636213 5.856580
## X37225_at X35162_s_at X39424_at X40480_s_at X34699_at X40454_at
## ALL1/AF4 5.220668 4.398885 6.896571 6.414368 4.253504 4.007171
## BCR/ABL 3.460902 4.924553 8.054597 8.208263 6.315966 3.910912
## E2A/PBX1 7.445655 4.380962 7.689701 6.722296 6.102031 7.390283
## NEG 3.387552 4.236335 8.044947 7.362318 6.092511 3.807652
## X37579_at X1616_at X39837_s_at X39929_at X41237_at X1307_at
## ALL1/AF4 7.614200 4.059859 6.633188 6.694605 10.94079 3.368915
## BCR/ABL 8.231081 3.643486 7.374046 7.664980 12.11895 4.945270
## E2A/PBX1 9.494368 6.033233 6.708400 9.313867 11.35610 4.678577
## NEG 8.455750 3.747773 6.878846 7.372117 11.93624 4.863930
## X37039_at X37027_at X37951_at X36873_at X37981_at X38994_at
## ALL1/AF4 11.12880 9.118515 3.418433 7.040593 6.170311 8.939007
## BCR/ABL 12.22497 9.421987 3.881780 3.490262 6.882755 8.760193
## E2A/PBX1 11.49101 6.688977 3.461861 3.634471 8.080002 4.741548
## NEG 12.00980 7.408175 3.419113 3.824670 7.423079 8.046647d <- gather(dt[,c(rf.genes,"Mut")],Gene,ExprValue,1:length(rf.genes))
dat <- group_by(d,Mut,Gene) %>%
summarise(med=median(ExprValue), iqr=IQR(ExprValue))## Warning: package 'bindrcpp' was built under R version 3.3.3ggplot(dat, aes(x=med,y=iqr,color=Mut)) +
geom_point(size=4) + facet_wrap(~ Gene) +
labs(x="MEDIAN expression level",y="IQR expression level",color="Mutation")
상기 그래프를 보면 하나의 유전자내에 각각의 변형에 따른 중앙값이 차이가 큰 경우를 볼 수 있다. X41470_at 의 ALL1/AF4 또는 X40454_at의 E2A/PBX1이 그렇다.
변수 선택법 4: 군집 앙상블(Clustering Ensembles)
다른 변수 선택법으로 군집 앙상블이 있다. 앙상블에는 데이터를 추출하는 방법과 설명변수를 추출하는 방법이 있다. 일반적으론 전자를 사용하나 설명변수가 충분히 있으므로 여기선 후자를 선택한다. 즉, 변수간 거리를 군집화 하는 것이다. 앙상블이라 불리는 이유는 군집 분석을 통해 얻은 각 군집에서 하나의 변수만을 추출해 데이터 셋을 만들기 때문이다. 변수간의 군집화를 통해 94개의 데이터가 서로 유사한지 아닌지를 알 수 있다. 변수 군집화는 varclus()함수로 구현 가능하다. cutree()함수는 군집 갯수를 만들어 준다. 각각의 군집에 포함된 변수 갯수는 다음과 같다.
vc <- varclus(t(es))
clus30 <- cutree(vc$hclust, 30)
table(clus30)## clus30
## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
## 18 34 30 22 34 35 19 16 40 52 19 22 17 24 30 26 20 17 18 21 43 30 32 14 23
## 26 27 28 29 30
## 28 18 17 11 16각각의 군집에서(이 예에서는 30개) 하나의 변수를 무작위로 추출하자. 아래는 각기 다른 군집에서 1개씩 뽑은 30개 변수 set의 결과이다. 매 시도 때마다 set.seed가 없으므로 변수명은 달라진다.
getVarsSet <- function(cluster,nvars=30,seed=NULL,verb=FALSE) {
if (!is.null(seed)) set.seed(seed)
cls <- cutree(cluster,nvars)
tots <- table(cls)
vars <- c()
vars <- sapply(1:nvars,function(clID)
{
if (!length(tots[clID])) stop('Empty cluster! (',clID,')')
x <- sample(1:tots[clID],1)
names(cls[cls==clID])[x]
})
if (verb) structure(vars,clusMemb=cls,clusTots=tots)
else vars
}
getVarsSet(vc$hclust)## [1] "266_s_at" "33598_r_at" "41138_at" "41610_at" "40468_at"
## [6] "1914_at" "41742_s_at" "39044_s_at" "430_at" "32961_at"
## [11] "38651_at" "1065_at" "40635_at" "555_at" "34740_at"
## [16] "1453_at" "34892_at" "38062_at" "35670_at" "34789_at"
## [21] "317_at" "39358_at" "34785_at" "39058_at" "38413_at"
## [26] "32562_at" "947_at" "37326_at" "1616_at" "32725_at"Metric 선정
예측을 위해서는 Metric 선정이 필요하다. 앞서 언급했듯이 이 문제는 ALL1/AF4, BCR/ABL, E2A/PBX1, NEG의 총 4개의 반응변수를 가진 Multi-Class 분류이다. 분류 문제는 일반적으로 분류 오류율 또는 ROC로 측정하는데 multi-class 분류(2개 이상 반응변수)의 경우 ROC는 정확치 않다고 알려져 있다. 따라서 여기서는 표준화된 정확도를 기준으로 한다. 이 값은 모든 데이터를 정확히 예측하면 1 모두 틀리면 0으로 각각이 오류가 발생할 때마다 1에서 점수가 차감된다.
검증법 선정
데이터 수가 94개 뿐이므로 bootstrap을 사용하자.
Modeling 적용 및 비교 그리고 경우의 수
위에서 언급한 변수 선택 방법을 여러가지 모델링에 적용하여 최상의 모델링을 찾을 시간이다. 변수 선택법의 경우, 1번 방법을 제외한 나머지 3가지 방법은 모두 반응 변수와 관련이 있다. 따라서 1번 방법 이후의 방법을 서로 비교하고자 한다. 3가지 모델은 - 1) ANOVA 방법(상기 변수선택법 2) - 2) ANOVA + RandomForest(상기 변수선택법 2+3) - 3) ANOVA + Clustering Ensemble (상기 변수선택법 2+4) 이다.
각각의 변수 선택법을 k-근접이웃(knn), 서포트 벡터머신(svm), 랜덤포리스트의 3가지 모델을 통해 예측하고 예측률을 비교하자.
각각의 모델은 최적의 값을 찾기 위해 여러 가지 변수를 사용하는데 예를 들어 랜덤포리스트의 경우 트리 갯수를 500, 700, 1000개로 하고 변수 갯수를 5개와 15개를 한다. 다시 말해 랜덤포리트의 모델 결과는 6개가 되고 이것을 3개의 변수 선택법에 각각 적용하면 총 18개의 결과가 나온다. 유사한 방법으로 KNN도 18개 , SVM도 18개의 결과를 얻는다. 총 54개의 결과는 검증이 필요한데 검증은 100번 반복 복원추출하는 bootstrap으로 한다.
생각만해도 어마어마한 숫자들이다. 따라서 내 컴퓨터로는 계산되지 않을 것임을 알 수 있다. 다행히 이 문제에서는 각각의 모델에 대한 결과값을 제공한다.
load("knn.Rdata")
load("svm.Rdata")
load("randomForest.Rdata")rankWorkflows(svm, maxs = TRUE) ## $ALL
## $ALL$acc
## Workflow Estimate
## 1 svm.v8 0.8126319
## 2 svm.v12 0.8112365
## 3 svm.v10 0.8064391
## 4 svm.v6 0.8046412
## 5 svm.v7 0.7978988rankWorkflows(randomForest, maxs = TRUE) ## $ALL
## $ALL$acc
## Workflow Estimate
## 1 randomForest.v7 0.7831709
## 2 randomForest.v8 0.7827674
## 3 randomForest.v9 0.7817705
## 4 randomForest.v11 0.7795233
## 5 randomForest.v12 0.7788864rankWorkflows(knn, maxs = TRUE) ## $ALL
## $ALL$acc
## Workflow Estimate
## 1 knn.v7 0.808435
## 2 knn.v8 0.808435
## 3 knn.v9 0.808435
## 4 knn.v10 0.808435
## 5 knn.v11 0.8084353가지 모델 결과 서포트벡터머신과 k-근접이웃의 정확도가 높다. 세가지를 하나로 묶어 순위를 매겨 보자.
all.trials <- mergeEstimationRes(svm, knn, randomForest, by ="workflows")Top 10 결과를 보면 서포트벡터 머신의 8번째 모델이 최고 모델이다.
rankWorkflows(all.trials, top=10, maxs = TRUE)## $ALL
## $ALL$acc
## Workflow Estimate
## 1 svm.v8 0.8126319
## 2 svm.v12 0.8112365
## 3 knn.v7 0.8084350
## 4 knn.v8 0.8084350
## 5 knn.v9 0.8084350
## 6 knn.v10 0.8084350
## 7 knn.v11 0.8084350
## 8 knn.v12 0.8084350
## 9 svm.v10 0.8064391
## 10 svm.v6 0.8046412서포트벡터머신의 8번째 값을 그래프로 나타내보자.
getWorkflow("svm.v8", all.trials)## Workflow Object:
## Workflow ID :: svm.v8
## Workflow Function :: ALLb.wf
## Parameter values:
## learner.pars -> cost=100 gamma=0.01
## learner -> svm
## featSel.meth -> s28번째 모델은 cost가 100, gamma가 0.01이다. s2는 랜덤포리스트를 통해 중요도 높은 변수 30개를 선택했다는 의미이다.
top10WFnames <- rankWorkflows(all.trials, top=10,
maxs = TRUE)[["ALL"]][["acc"]][,"Workflow"]
sapply(top10WFnames, function(WFid) getWorkflow(WFid,all.trials)@pars$featSel.meth)## svm.v8 svm.v12 knn.v7 knn.v8 knn.v9 knn.v10 knn.v11 knn.v12 svm.v10
## "s2" "s2" "s2" "s2" "s2" "s2" "s2" "s2" "s2"
## svm.v6
## "s1"plot(subset(all.trials,workflows=top10WFnames))
마지막으로 이 결과의 차이가 통계적 유의성이 있는지 확인하자. pairedComparisons 함수는 Wilcoxon signed rank test로 baseline이 되는 결과와 나머지 결과간의 1:1 비교를 통해 통계적 유의성이 있는지를 볼 수 있다. 결과를 보면 10개 모두 결과 차이가 유의성이 없음을 알 수 있다.
ps <- pairedComparisons(subset(all.trials,workflows=top10WFnames),baseline="svm.v8")## Warning in pairedComparisons(subset(all.trials, workflows =
## top10WFnames), : With less 2 tasks the Friedman, Nemenyi and Bonferroni-
## Dunn tests are not calculated.ps$acc$WilcoxonSignedRank.test## , , ALL
##
## MedScore DiffMedScores p.value
## svm.v8 0.8235294 NA NA
## svm.v12 0.8235294 0.000000000 0.7000463
## knn.v7 0.8169856 0.006543766 0.7941389
## knn.v8 0.8169856 0.006543766 0.7941389
## knn.v9 0.8169856 0.006543766 0.7941389
## knn.v10 0.8169856 0.006543766 0.7941389
## knn.v11 0.8169856 0.006543766 0.7941389
## knn.v12 0.8169856 0.006543766 0.7941389
## svm.v10 0.8086312 0.014898200 0.2240278
## svm.v6 0.8055556 0.017973856 0.4055620confusion matrix를 통해 예측 결과를 살펴보자. 100번의 bootstrap 검 중 한개에 대한 예측 값을 보자
iteration <- 1 # any number between 1 and 100 in this case
itInfo <- getIterationsInfo(all.trials,workflow="svm.v8",it=iteration)
table(itInfo$trues, itInfo$preds)##
## ALL1/AF4 BCR/ABL E2A/PBX1 NEG
## ALL1/AF4 3 0 0 0
## BCR/ABL 0 12 0 3
## E2A/PBX1 0 0 0 1
## NEG 0 1 1 14한가지 의문은 Bootstrap은 100개의 평균일 경우에 장점을 가지는데 위와 같이 1번의 결과를 보는게 무슨 의미가 있을까 싶다.
결론
Algae 문제에서 Torgo가 제공하는 궁극의 패키지 performanceEstimation을 잘 이용하면 내가 짠 함수와 함께 유용하게 쓸 수 있다고 생각했었다. Torgo 패키지에서 제공하지 않는 정규화 선형회귀(Ridge, Lasso)를 제외하고는 나머지 모델의 경우 상당한 이점이 있다고 생각했다. (굳이 내가 함수를 만들어야 하는 수고를 피할 수 있으니까!!) 그러나 본 예제에서 보듯, 다른 사람이 만든 함수(또는 패키지)를 능수능란하고 자유롭게 쓰려면 상당한 Technical 숙련도가 상당해야 함을 알 수 있다. 특히 각 모델별 18개씩의 결과를 저장하는 함수 작성은 이해가 쉽지 않다.
그럼에도 본 예제를 통해 다음의 3가지의 중요한 사항을 배울 수 있다.
- 데이터보다 변수가 많을 경우, 그리고 단위가 같을 경우, IQR의 Low Variability, ANOVA평균, 랜덤포리스트의 중요도를 통해 변수 제거가 가능하다.
- k-평균 군집을 할 때 데이터간의 거리가 아닌 변수간 거리를 해도 된다. 그리고 이것을 변수 선택의 방법으로 사용한다.
- 랜덤포리스트, 부스팅 외 또다른 앙상블인 군집 앙상블을 배웠다. 랜덤포리스트가 추출을 먼저(bootstrap)하고 모델(tree)을 생성하는 반면 군집앙상블은 군집(clustering)을 형성한 후 sample 추출을 하는점이 다르다.
마지막으로 교훈은 Torgo 패키지를 손에 익도록 자주 써야 겠다.
'R-PAGO 노트' 카테고리의 다른 글
| Entropy Based Discretization: 어떻게 연속변수를 구분 짓는가? (0) | 2017.06.13 |
|---|---|
| 능형(Ridge) 회귀, Lasso 이론 및 예제 (0) | 2017.04.15 |
| 부스팅(Boosting) 이론과 예제 (0) | 2017.04.09 |
| 랜덤포리스트(Random Forest) 총정리(배깅 포함) (0) | 2017.04.01 |
| 검정오차(test error)가 중요한 이유 (0) | 2017.03.30 |