待ち人来る

bioconductor

各所で話題になっている The R Book がウチにも届いた．パラパラと斜め読みしていたときに，最近ハマッたことに関わる部分があった．p. 372 なんだけど

Calls <- mas5calls(affybatch.example)

ってのがある．通常はコレで問題なく動くし，bioconductor 1.3 以前でも動く(動いていた)んだけど，bioconductor 1.4 にしたときに動かなくなる場合があった．具体的には，cell ファイルを一個だけ読み込んだ AffyBatch オブジェクトに対しては mas5calls が動かない．例えば，

# cell ファイルを一つだけ読み込む．複数読み込んだ場合はちゃんと動く
dataSet <- ReadAffy(filenames = c("hoge.cel")) # 落ちる例
# dataSet <- ReadAffy(filenames = c("hoge.cel", "moga.cel")) # 上手くいく例
Calls <- mas5calls(dataSet) # こいつがコケる

なんてスクリプトがコケるのだ．ソースコード(mas5.R ver.1.4.30)を見れば分かるんだけど，mas5calls の 27行目付近

pms <- pm(object)
mms <- mm(object)
.... 中略 ....
pvals<-sapply(1:length(pms[1,]),function(x) {

pm(object) の返す型が dataSet 無いの cell ファイル数で異なる(ザックリ言えば，一つなら配列で複数なら行列) のが原因だった．配列に対して pms[1,] のような行列みたいにアクセスして，落ちていた．

解決方法

いくつか考えられる．一つは

dataSet <- ReadAffy(filenames=c("hoge.cel"))
dummy <- merge.AffyBatch(dataSet, dataSet) # わざと複数にしちゃう
dummyCalls <- mas5calls(dummy)
write.table(cbind(exprs(dummyCalls)[1,], se.exprs(dummyCalls)[1,]), "hoge.calls")

みたいに，わざと複数あるように見せちゃう方法．もう一つは，mas5calls の中身をパクッて

singleMas5calls <- function(object){
# 全部書くのしんどいからパス．sapply のところを
#.C("DetectionPValue",as.double(pms),as.double(mms),as.character(pns),as.integer(length(mms)),
#as.double(tau),as.double(sat),dpval=double(length(unique.pns)),length(unique.pns), #PACKAGE="affy")$dpval;
# みたいにすれば OK
}

ってやる方法．最後に pms/mms って関数を書き換える方法が挙げられる．さいごの奴が出来ればパッチを投げるんだけどなぁ．そこまでソースコードを読み込んでないからなぁ．