リファレンスに基づいたクラスタリングの危険性

電子顕微鏡による単粒子解析界隈でちょっとした騒ぎ(?)になっている事件がある。"Avoiding the pitfalls of single particle cryo-electron microscopy: Einstein from noise" (Richard Henderson, PNAS 2013) と、このページからリンクされている Related Letters だ。

＃　単粒子解析の流れと、今回の問題を簡単に説明しておこうと思ったが、面倒になったまま、この記事も下書きに1ヶ月以上放置してしまったので、説明をせずに、実験部分だけ公開してしまう

これを実体験するため、Rで簡単な実験を行った。まず、ニワトリのアイコンを GIMP で 8bit のグレースケールPNG 画像に変換した。これを参照画像(リファレンス)とする。同サイズのランダム・ノイズ画像を10万枚作成し、参照画像との相関係数を求める。相関係数は 0 を中心とする左右対称な分布となった(理論的には何分布になるんだっけ?)。

ここから、相関係数が 0.005 以上だった画像(35716枚あった)を集めて平均化したものを表示した。

ご覧のとおり、単なるのノイズから参照画像が復元されてしまった!! 実際、この平均画像とリファレンスの相関係数を調べると、0.937 にも達した。

なお、image 命令は勝手にコントラストを調節して色をつけてくれる。実際の値のヒストグラムを描くと

となっていて、S/N比はかなり悪いのが分かる。だが、S/Nが悪かろうがなんだろうが、ノイズから「それらしい」画像が出てきてしまうのが落とし穴なのである。
なお、相関係数 0.005 という閾値には深い意味はない。全体として相関係数が正になるように画像を集めれば、参照画像に似たものが出てきてしまう。

コードは以下のとおり。計算時間は1分ほど。

library('png')
ref <- readPNG('reference.png')
ref.vec <- as.vector(ref)

X <- dim(ref)[1]
Y <- dim(ref)[2]
N <- 100000

thresh <- 0.005
filtered <- rep(0, length(ref.vec))
n.summed <- 0
cc <- rep(NA, N)
for (i in 1:N) {
  if (i %% 5000 == 0) {
    print(i)
  }
  noise <- runif(X * Y)
  cc[i] <- cor(noise, ref.vec)
  if (cc[i] > thresh) {
    n.summed <- n.summed + 1
    filtered <- filtered + noise
  }
}

hist(cc)
n.summed
filtered <- filtered / n.summed

par(mfrow=c(1, 2), pty='s')
image(matrix(filtered, X, Y), col=gray((0:255)/255), 
      useRaster=TRUE, pty='s', main='Reconstructed', xaxt='n', yaxt='n')
image(matrix(ref, X, Y), col=gray((0:255)/255),
      useRaster=TRUE, main='Original', xaxt='n', yaxt='n')

hist(filtered, main='Reconstructed', xlab='intensity')
hist(ref, main='Original', xlab='intensity')

cor(ref.vec, filtered)

なお、アニメ絵のようにコントラストが高い画像の場合は、試行を1万枚まで減らしても十分に見えてくる。もともとは、ゆるゆり公式Webサイトが配っている赤座あかりの Twitter アイコンで実験していた。