ほくそ笑む

R言語と統計解析について

samr のバグとその解決法について考えた

SAM (Significance Analysis of Microarrays) は、マイクロアレイ発現解析の一手法で、対照実験において発現の有意な遺伝子を選定するための統計解析手法です。
t検定の改良版みたいなやつだと思ってください。
この SAM は統計言語 R では、samr パッケージで行うことが出来ます。
詳しくは門田先生のサイトを見ていただくとして、以下のような感じで実行できます。

install.packages("samr")
library(samr)
data <- read.table("mydata.txt", header=TRUE, row.names=1, sep="\t", quote="")
data.cl <- c(rep(1, 5), rep(2, 5))
data.tmp = list(x=as.matrix(data), y=data.cl, geneid=rownames(data), genenames=rownames(data), logged2=TRUE)
out <- samr(data.tmp, resp.type="Two class unpaired", nperms=200)

samr のバグ

弊社のマイクロアレイ受託解析では、基本解析プラン(一番安いやつ)で、t検定や他の検定と一緒にこの SAM も行っているのですが、ここ3年ほど全く問題のなかった解析プログラムが突然警告を吐きだしました。

 警告メッセージ: 
1: In init.fit$sd < s0 :
   長いオブジェクトの長さが短いオブジェクトの長さの倍数になっていません 
2: In sd + s0 :
   長いオブジェクトの長さが短いオブジェクトの長さの倍数になっていません 

これは一体何が原因なのでしょうか?
SAM では、統計量を計算する際、発現値の大きいものと小さいもので同じ統計量を使って比較するのは不公平ということで、補正項 s0 を挿入しています。
この s0 は、変動係数が最小になるように選ぶのですが、samr() 関数は、変動係数の最小値が2つ以上ある場合を想定せずにプログラムが作成されています。
そのため、変動係数の最小値が2つ以上算出されるようなデータを入力すると、s0 がベクトルとなり、s0 を用いた計算の部分で上記のような警告が出てしまうのです。
これは samr() 関数のバグです。

解決法

それでは、このバグを回避する方法を考えてみましょう。
変動係数の最小値が2つ算出される原因として、次の2つが考えられます。

  1. 本当は最小値は一つしかないが、s0 の計算精度が悪くて二つ出る
  2. 本当に最小値が二つある

samr() 関数では、s0 の推定に est.s0() という内部関数が使われています。
先程のデータを使って実際に s0 を計算してみます。

init.fit <- samr:::ttest.func(data, data.cl)
s0.hat <- samr:::est.s0(init.fit$tt, init.fit$sd)$s0.hat
s0.hat
        0%         5%        10% 
0.03333333 0.13579545 0.26666667 

確かに s0 の推定値が3つも出てしまっています。
s0 の計算精度は est.s0() 関数の s0.perc パラメータで指定できます。
少し細かく計算するように指定してみましょう。

step <- 0.01 # default:0.05
init.fit <- samr:::ttest.func(data, data.cl)
s0.hat <- samr:::est.s0(init.fit$tt, init.fit$sd, seq(0, 1, by = step))$s0.hat
s0.hat
      13% 
0.3333333 

s0 の推定値が一つになりました。つまり、原因は計算精度にあったわけです。
ところが、もっと計算精度を細かくしてみると、別の問題が出てきます。

step <- 0.001 # default:0.05
init.fit <- samr:::ttest.func(data, data.cl)
s0.hat <- samr:::est.s0(init.fit$tt, init.fit$sd, seq(0, 1, by = step))$s0.hat
s0.hat
    12.5%     12.6%     12.7%     12.8%     12.9%       13%     13.1%     13.2%     13.3%     13.4% 
0.3282953 0.3333333 0.3333333 0.3333333 0.3333333 0.3333333 0.3333333 0.3333333 0.3333333 0.3349959 

計算精度が細かすぎると、推定値がたくさん出てきてしまうのです。
そうですよね。よくよく考えたら、計算精度を無限に細かく出来ない限り、原因が計算精度なのか、本当に最小値が二つあるのかの切り分けは不可能ですよね。
結局、いずれの原因にせよ、適切な計算精度で s0 の推定値を算出する必要があるようです。

まとめ

samr() 関数を使っていて上記の警告が出た場合、データによって適切な計算精度を指定し、s0 を算出しなおすことが必要です。
samr() 関数のパラメータには、s0 を直接指定することができるので、上記で求めた s0.hat を入力してもう一度 samr() を実行すればいいでしょう。

out <- samr(data.tmp, resp.type="Two class unpaired", nperms=200, s0=s0.hat)

以上です。