Bioconductorノート

■ ケーススタディ：遺伝子発現データに基づいた分類

つぎのclass.Rを参考にして、遺伝子発現データに基づいた分類についてのケーススタディを説明します。

• http://www.bioconductor.org/mogr/chapter-code/class.R

▼のついてソースコードをRコンソールにコピー＆ペーストしながらすすめていきます。

ストーリー

1. データ：腎細胞ガンの遺伝子発現データ
2. 目的：遺伝子発現データから細胞の表現形（ガンの種類）を分類予測する
3. 目的２：うまく分類できる分類法を模索する
4. データの読み込み
5. 変数選択で利用する遺伝子を選択
6. 複数の分類方法（DLDA、kNN、SVM）を評価
7. 複数回のランダム分割によるパラメータ最適化
8. 細胞の分類予測をおこなう
9. 出口：適切な分類方法とそれによる細胞の分類

腎細胞ガンデータの読み込み

GeneChip による腎臓ガン細胞遺伝子発現データは、kidpack データパッケージの eset データを使用します。

▼データの読み込み。

set.seed(32)
library("kidpack")
data(eset)

▼あとで利用するためにデータセットから一部取り除いておきます。

x        <- t(exprs(eset)) 
y        <- pData(phenoData(eset))$type
y        <- as.factor(y)
y.train  <- y[-c(1:10) ]
x.train  <- x[-c(1:10),]  

▼eset の中身

eset

△実行結果

ExpressionSet (storageMode: lockedEnvironment)
assayData: 4224 features, 74 samples 
  element names: exprs 
phenoData
  sampleNames: 1, 2, ..., 87 (74 total)
  varLabels and varMetadata:
    lfdnr: NULL
    patientid: NULL
    ...: ...
    22q: NULL
    (64 total)
featureData
  featureNames: 1, 2, ..., 4224 (4224 total)
  varLabels and varMetadata: none
experimentData: use 'experimentData(object)'
Annotation [1] ""

▼eset は Biobase パッケージに含まれる exprSet クラスのデータです。

class(eset)

△実行結果

[1] "ExpressionSet"
attr(,"package")
[1] "Biobase"

▼細胞型は $type でアクセスできます。

eset$type

△実行結果

 [1] "ccRCC" "ccRCC" "ccRCC" "ccRCC" "ccRCC" "pRCC"  "ccRCC"
 [8] "pRCC"  "chRCC" "ccRCC" "ccRCC" "ccRCC" "chRCC" "chRCC"
[15] "ccRCC" "ccRCC" "ccRCC" "ccRCC" "ccRCC" "chRCC" "pRCC" 
[22] "ccRCC" "ccRCC" "ccRCC" "ccRCC" "ccRCC" "ccRCC" "ccRCC"
[29] "pRCC"  "ccRCC" "ccRCC" "ccRCC" "pRCC"  "ccRCC" "chRCC"
[36] "chRCC" "ccRCC" "ccRCC" "ccRCC" "chRCC" "pRCC"  "pRCC" 
[43] "ccRCC" "pRCC"  "pRCC"  "ccRCC" "pRCC"  "ccRCC" "ccRCC"
[50] "pRCC"  "ccRCC" "ccRCC" "ccRCC" "ccRCC" "ccRCC" "ccRCC"
[57] "ccRCC" "ccRCC" "ccRCC" "ccRCC" "ccRCC" "ccRCC" "pRCC" 
[64] "ccRCC" "chRCC" "ccRCC" "ccRCC" "ccRCC" "ccRCC" "ccRCC"
[71] "chRCC" "pRCC"  "ccRCC" "ccRCC"

▼三種類の細胞型のそれぞれの数の分布をみます。

table(eset$type)

△実行結果

ccRCC chRCC  pRCC 
   52     9    13

分類器の訓練と評価の準備をする

分類には、訓練データ、学習データと評価データの三種類のデータを用意する必要があります。それぞれを学習データと評価データは訓練データを 2:1 に分割して作成されます。

▼訓練データサイズ（t.size）、学習データサイズ（l.size）、評価データサイズ（v.size）

t.size   <- length(y.train)
l.size   <- round((2/3) * t.size) # 2/3 が l.size へ
v.size   <- t.size - l.size       # 1/3 が v.size へ

▼訓練データ（t.samp）、学習データ（l.samp）、評価データ（v.samp）

l.samp   <- NULL
K        <- nlevels(y.train)
props    <- round(l.size / t.size * table(y.train))
props[1] <- l.size - sum(props[2:K])
for (k in 1:K)
  {
    y.num  <- as.numeric(y.train)
    l.samp <- c(l.samp, sample(which(y.num == k))[1:props[k]])
  }
v.samp     <- (1:t.size)[-l.samp]

データ集合をみる

▼学習集合（l.samp）を表示してみる。

table(y.train[l.samp])

△実行結果

ccRCC chRCC  pRCC 
   31     5     7 

▼評価集合（v.samp）を表示してみる。

table(y.train[v.samp])

△実行結果

ccRCC chRCC  pRCC 
   14     3     4

このようにデータの比率は保存されています。

変数選択で利用する遺伝子を選択する

▼効果的な分類のために分類予測に使用する遺伝子を選別します。ここでは単純なF統計量（群間平方和と群内平方和の比）を用います。F統計量の上位200遺伝子を選択します。

library("multtest")
yl.num     <- as.numeric(y.train[l.samp]) - 1
xl.mtt     <- t(x.train[l.samp,])
f.stat     <- mt.teststat(xl.mtt, yl.num, test = "f")
best.genes <- rev(order(f.stat))[1:200]
x.learn    <- x.train[l.samp, best.genes]
x.valid    <- x.train[v.samp, best.genes]
y.learn    <- y.train[l.samp]
y.valid    <- y.train[v.samp]

遺伝子発現のF統計量が大きいとは、群内での分散が小さくかつ群間での分散が大きい、すなわち、群（ここでは分類対象：細胞タイプ）に特徴的な遺伝子であるとしています。

分類を実行

分類には対角線形判別（DLDA）、k最近傍分類（kNN）、サポートベクターマシン（SVM）の三種類のアルゴリズムを試してみます。

▼対角線形判別（diagonal linear disriminant analysis: DLDA）の実行

library("sma")
yl.num      <- as.numeric(y.learn)
yv.num      <- as.numeric(y.valid) - 1
dlda.predic <- stat.diag.da(x.learn, yl.num, x.valid)
dlda.error  <- mean(dlda.predic$pred - 1 != yv.num)

▼k最近傍分類（kNN）の実行。kNNではパラメータ k を 1, 3, 5 に振って最適な分類を探索します。

library("class")
knn.error <- numeric(3)
for(k in c(1,3,5))
  {
    i              <- ((k - 1) / 2) + 1
    knn.predic     <- knn(x.learn, x.valid, y.learn, k = k)
    knn.error[i]   <- mean(knn.predic != y.valid)
  }

▼サポートベクターマシン（SVM）の実行。SVM では二つのパラメータ（γと c）を振って最適なパラメータを探索します。パラメータγ(gamma）は -1, 0, 1に、パラメータc（cost） は 0, 1, 2 に振っています。実際のパラメータは，cost = 2^cost なので、{2^0, 2^1, 2^2}です。また、gamma = 2^gamma/ncol(x.learn) なので、{2^(-1)/200, 2^0/200, 2^1/200} です。

library("e1071")
svm.error <- matrix(0, nrow = 3, ncol = 3)
for(cost in 0:2)
  {
    for(gamma in (-1):1)
      {
        i              <- cost + 1
        j              <- gamma + 2
        svm.fit        <- svm(x.learn, y.learn, cost = 2^cost,    
                              gamma = 2^gamma/ncol(x.learn))
        svm.predic     <- predict(svm.fit, newdata = x.valid)
        svm.error[i,j] <- mean(svm.predic != y.valid)
      }
  }

分類器 DLDA、kNN、SVM の誤分類率を表示

▼DLDAの誤分類率を表示

dlda.error

△実行結果

[1] 0

▼kNNの誤分類率を表示

knn.error

△実行結果

[1] 0.04761905 0.00000000 0.04761905

パラメータ k = 3 の時がもっとも正確な分類をおこなっています。

▼SVMの誤分類率を表示

svm.error

△実行結果

           [,1]       [,2]      [,3]
[1,] 0.04761905 0.04761905 0.0952381
[2,] 0.04761905 0.04761905 0.0952381
[3,] 0.04761905 0.04761905 0.0952381

行がパラメータc {0,1,2} で列がパラメータγ{-1,0,1}です。

c が {0,1} かつ γ{-1,0} の組み合わせが、SVM のなかでもっとも正確な分類をおこなっています。

ここで三つの分類法を比較しました。その結果、DLDA と k = 3 の kNN がもっとも正確な分類をおこなうことがわかりました。今後の解析（診断）ではどちらかの方法をつかうとよいといえます。

しかし、今回の訓練データの分割がたまたま DLDA や k = 3 の kNN に有利ダッタのかもしれません。このような偶然性への依存を減らすには訓練データの分割を複数回おこなうことが有効です。

複数回のランダム分割

上記のようなパラメータの最適化をおこなう場合は、データの分割を複数回ためす必要があります。ここでは、訓練データの50 回のランダム分割によるパラメータ最適化をみてみます。

▼

randiv.repr <- function(x, y)
  {
    t.size   <- length(y)
    l.size   <- round((2 / 3) * t.size)
    v.size   <- t.size - l.size
    l.samp   <- NULL
    K        <- nlevels(y)
    props    <- round(l.size / t.size * table(y))
    props[1] <- l.size - sum(props[2:K])
    for (k in 1:K)
      {
        y.num   <- as.numeric(y)
        l.samp  <- c(l.samp, sample(which(y.num == k))[1:props[k]])
      }
    v.samp      <- (1:t.size)[-l.samp]
    yl.num      <- as.numeric(y[l.samp]) - 1
    f.stat      <- mt.teststat(t(x[l.samp,]), yl.num, test = "f")
    best.genes  <- rev(order(f.stat))[1:200]
    x.learn     <- x[l.samp, best.genes]
    x.valid     <- x[v.samp, best.genes]
    y.learn     <- y[l.samp]
    y.valid     <- y[v.samp]
    yl.num      <- as.numeric(y.learn)
    yv.num      <- as.numeric(y.valid) - 1
    dlda.predic <- stat.diag.da(x.learn, as.numeric(y.learn), x.valid)
    dlda.error  <- mean(dlda.predic$pred-1 != yv.num)
    knn.error   <- numeric(3)
    for(k in c(1,3,5))
      {
        i              <- ((k - 1) / 2) + 1
        knn.predic     <- knn(x.learn, x.valid, y.learn, k = k)
        knn.error[i]   <- mean(knn.predic != y.valid)
      }
    svm.error <- matrix(0, nrow = 3, ncol = 3)
    for(cost in 0:2)
      {
        for(gamma in (-1):1)
          {
            i              <- cost + 1
            j              <- gamma + 2
            svm.fit        <- svm(x.learn, y.learn, cost = 2^cost,    
                                  gamma = 2^gamma / ncol(x.learn))
            svm.predic     <- predict(svm.fit, newdata = x.valid)
            svm.error[i,j] <- mean(svm.predic != y.valid)
          }
      } 
    list(dlda = dlda.error, knn = knn.error, svm = svm.error,
         lsamp = l.samp, vsamp = v.samp)
  }

▼データの分割と分類を50回くり返す

runs         <- 50
dlda.errors  <- numeric(runs)
knn.errors   <- matrix(0, nrow = 3, ncol = runs)
svm.errors   <- array(0, dim = c(3, 3, runs))
vsamp        <- matrix(0, runs, length(v.samp))
lsamp        <- matrix(0, runs, length(l.samp))
for(r in 1:runs)
  {
    results         <- randiv.repr(x.train, y.train)
    dlda.errors[ r] <- results$dlda
    knn.errors[, r] <- results$knn
    svm.errors[,,r] <- results$svm
    vsamp[r,]       <- results$vsamp
    lsamp[r,]       <- results$lsamp
    cat("This was run", r, "of", runs, "\n")
  }

△実行結果

This was run 1 of 50 
This was run 2 of 50 
This was run 3 of 50 
This was run 4 of 50 
This was run 5 of 50 
This was run 6 of 50 
This was run 7 of 50 
This was run 8 of 50 
This was run 9 of 50 
This was run 10 of 50 
This was run 11 of 50 
This was run 12 of 50 
This was run 13 of 50 
This was run 14 of 50 
This was run 15 of 50 
This was run 16 of 50 
This was run 17 of 50 
This was run 18 of 50 
This was run 19 of 50 
This was run 20 of 50 
This was run 21 of 50 
This was run 22 of 50 
This was run 23 of 50 
This was run 24 of 50 
This was run 25 of 50 
This was run 26 of 50 
This was run 27 of 50 
This was run 28 of 50 
This was run 29 of 50 
This was run 30 of 50 
This was run 31 of 50 
This was run 32 of 50 
This was run 33 of 50 
This was run 34 of 50 
This was run 35 of 50 
This was run 36 of 50 
This was run 37 of 50 
This was run 38 of 50 
This was run 39 of 50 
This was run 40 of 50 
This was run 41 of 50 
This was run 42 of 50 
This was run 43 of 50 
This was run 44 of 50 
This was run 45 of 50 
This was run 46 of 50 
This was run 47 of 50 
This was run 48 of 50 
This was run 49 of 50 
This was run 50 of 50 

DLDA、knn、svmにおける誤判別率の平均値を表示します。

▼DLDAの誤判別率の平均値を表示

mean(dlda.errors)

△実行結果

[1] 0.0467

▼kNNの誤判別率の平均値を表示

apply(knn.errors, 1, mean)

△実行結果

[1] 0.0590 0.0571 0.0600

▼SVMの誤判別率の平均値を表示

apply(svm.errors, c(1,2), mean)

△実行結果

       [,1]   [,2]  [,3]
[1,] 0.0657 0.0810 0.138
[2,] 0.0610 0.0705 0.128
[3,] 0.0562 0.0695 0.128

▼誤判別率の箱ヒゲ図と棒グラフを表示するメソッドを定義。

## Loading the required package
library("grid")

## Check the settings
if (!exists("pushViewport")) pushViewport <- push.viewport
if (!exists("popViewport"))  popViewport  <- pop.viewport

## The function for the barplots
grid.barplot <- function(x, xlab = FALSE, main = FALSE, xscale = NULL, yticks, ylabels)
{
  ## Customizing the data
  y      <- table(x)
  y.alt  <- table(x)
  x.alt  <- x
  x      <- sort(unique(x))
  if(!is.null(xscale))
    {
      index <- which(x > xscale[1] & x < xscale[2])
      x <- x[index]
      y <- y[index]
    }

  ## Horizontal bars
  grid.segments(x0 = unit(x, "native"), y0 = unit(0, "native"),
                x1 = unit(x, "native"), y1 = unit(y.alt, "native"))

  ## Drawing the lines
  grid.lines(unit(x, "native"), unit(y, "native"))

  ## Plotting the mean
  brks   <- as.numeric(attributes(table(x))$dimnames$x)
  start  <- max(which(brks < mean(x.alt)))
  ende   <- start + 1
  xinc   <- (mean(x.alt) - brks[start]) / (brks[ende] - brks[start])
  yh     <- ((y.alt[ende] - y.alt[start])*xinc) + (y.alt[start])             
  farbe  <- gpar(col="red")
  grid.segments(x0 = unit(mean(x.alt), "native"), y0 = unit(0, "native"),
                x1 = unit(mean(x.alt), "native"), y1 = unit(yh, "native"), gp = farbe)

  ## Annotation
  if(xlab) grid.xaxis()
  grid.yaxis(main = FALSE, at = yticks, name = ylabels)
  grid.rect()
}


## The function for the boxplots
grid.boxplot <- function(x, xlab = FALSE, main = FALSE, ylab = NULL, xscale = NULL)
{
  ## Generic boxplot
  x    <- boxplot(x, plot = FALSE)
  farb <- gpar(col="red")

  ## Transformation to the grid-system
  grid.lines(unit(x$stats[1],  "native"), unit(c(0.3, 0.7), "npc"))
  grid.lines(unit(x$stats[1:2],"native"), unit(0.5, "npc"), gp = gpar(lty = 2))
  grid.rect(unit(x$stats[2], "native"), y = unit(0.5, "npc"),
            height = unit(0.6, "npc"), width = unit(diff(x$stats[2:3]), "native"),
            just = c("left", "centre"))
  grid.rect(unit(x$stats[3], "native"), y = unit(0.5, "npc"),
            height = unit(0.6, "npc"), width = unit(diff(x$stats[3:4]), "native"),
            just = c("left", "centre"))
  grid.lines(unit(x$stats[4:5],"native"), unit(0.5, "npc"), gp = gpar(lty = 2))
  grid.lines(unit(x$stats[5], "native"), unit(c(0.3, 0.7), "npc"))
  grid.lines(unit(x$stats[3], "native"), unit(c(0.2, 0.8), "npc"), gp = farb)

  ## Outliers
  n <- length(x$out)
  if(n > 0) {
    if(!is.null(xscale)) index <- which(x$out > xscale[1] & x$out < xscale[2])
      else index <- 1:n
    if(length(index) > 0)
      {
        grid.points(unit(x$out[index],"native"),unit(rep(0.5,length(index)),
                    "npc"), size = unit(0.5, "char"))
      }
  }

  ## Annotation
  if(xlab) grid.xaxis()
  if(!is.null(ylab))
    {
      grid.text(ylab, x = unit(0, "npc") - unit(0.5, "lines"),
                gp = gpar(fontface = "bold"), just = c("right", "centre"))
    }
}


## Plotting the results
grid.boxNbar <- function(daten, mname)
  {
    grid.newpage()
    pushViewport(plotViewport(c(4, 6, 4, 4)))
    #pushViewport(viewport(layout = lyt))
    n         <- ncol(daten)

    ## Calibration of the x-axis
    xsc     <- range(daten) + c(-1, 1) * max(daten)/15
    xsc2    <- xsc          

    ## Calibration of the y-axis
    ysc     <- c(0, max(sapply(apply(daten, 2, function(x) table(x)), max)))
    ysc     <- ysc * 1.1
    ylabels <- c("0", "10", "20", "30", "40", "50")
    ynumb   <- floor((ysc[2] + 3) / 10) - 1
    yticks  <- (0:ynumb) * 10
    ylabels <- ylabels[1:(ynumb + 1)]
          
    ## The plot contains 3 pieces
    lyt <- grid.layout(1, 3, widths = unit(c(0.5 - 0.03 / 2, 0.03, 0.5 - 0.03 / 2), "npc"))
    pushViewport(viewport(layout = lyt))

    ## Left side: boxplots
    lyt <- grid.layout(n, 1, heights = unit(rep(1/n, n), "npc"))
    pushViewport(viewport(layout.pos.row = 1, layout.pos.col = 1, xsc = xsc))
    pushViewport(viewport(layout = lyt))
    grid.rect()
    for(i in n:1)
      {
        pushViewport(viewport(layout.pos.r = i, layout.pos.c = 1, xsc = xsc))
        main  <- (i < 2)
        grid.boxplot(daten[,i], m = main, xl = (i > (n - 1)), yl = mname[i], xsc = xsc)
        popViewport(1)
      }
    popViewport(1)
    grid.text("Boxplot", y = unit(-2.5, "lines"))

    ## Middle: an empty frame  
    popViewport(1)
    #schrift <- gpar(fontface = "bold")
    #grid.text(namen[j], y = unit(1, "npc") + unit(2, "lines"), gp = schrift)

    ## Right side: density plots
    lyt <- grid.layout(n, 1, hei = unit(rep(1/n, n), "npc"))
    pushViewport(viewport(layout.pos.row = 1, layout.pos.col = 3, xscale = xsc))
    pushViewport(viewport(layout = lyt))
    grid.rect()
    for(i in n:1)
      {
        pushViewport(viewport(layout.pos.r = i, layout.pos.c = 1, xs = xsc2, ys = ysc))
        main  <- (i < 2)
        grid.barplot(daten[,i], main, xsc2, yticks, ylabels, xlab = (i > (n - 1)))
        popViewport(1)
      }
    popViewport(1)
    grid.text("Density", y = unit(-2.5, "lines"))
    popViewport(3)
  }

▼図を生成。

daten <- cbind(dlda.errors, knn.errors[2,], svm.errors[3,1,], svm.errors[2,1,])
mname <- c("DLDA", "kNN", "SVM1", "SVM2")
grid.boxNbar(daten, mname)

△実行結果

kNN は k = 3 の場合、SVM1 は c = 2^2 かつ γ = 2^1/200、SVM2 は c = 2^2 かつ γ = 2^0/200の場合です。赤いバーは平均値を示しています。左のパネルが箱ヒゲ図で分類精度の分布を表していて、DLDA の分類精度がほかのものよりよいことが読み取れる。右のパネルでは密度分布図で分類精度の分布を表している。こちらでも同様である。

補足：図をファイルとして保存する場合は、pdf() と dev.off() を利用します。

pdf(file = 'box_and_barplot.pdf')
grid.boxNbar(daten, mname)
dev.off()

分類予測

はじめにデータから取り除いていた 10 個のサンプルに対しての細胞型の予測をおこないます。

▼DLDA による予測を実行

test   <- (1:10)
train  <- (1:length(eset$type))[-test]  # needed!
trEset <- eset[,train]
library("MLInterfaces")
yt.num      <- as.numeric(factor(trEset$type)) - 1
xt.mtt      <- exprs(trEset)
f.stat      <- mt.teststat(xt.mtt, yt.num, test = "f")
best.genes  <- rev(order(f.stat))[1:200]
selEset     <- eset[best.genes,]
dlda.predic <- stat.diag.daB(selEset, "type", train)
dlda.predic@RObject$pred - 1

△実行結果

 [1] 0 0 0 0 0 2 0 2 1 0

▼DLDA による予測結果を混合行列で評価。

confuMat(dlda.predic)

△実行結果

       predicted
given   1 2 3
  ccRCC 7 0 0
  chRCC 0 1 0
  pRCC  0 0 2

ここでは10個すべてのサンプルで正しい細胞型を予測しました。正しく予測した場合、対角成分に数字がならびます。対角成分以外に数字がある場合は誤分類（偽陽性もしくは偽陰性）の数を示しています。

まとめ

ここでは、遺伝子発現データから細胞を分類予測する研究をみました。

このような解析はあたらしい診断、分類や予測方法の開発でもちいられます。

Bioconductor では数多くの分類法が利用できます。

Bioconductor には遺伝子発現データのサンプルデータが豊富にあります。

複雑なことをするにはそれなりに長いプログラムが必要です。

もくじにもどる