組み合わせにより作成する2x2の分割表の周辺度数

[ryamada22]

成分が{0,1}であるnSample x nMarker行列型のデータを考える。これはZDD的にはtransactionデータのこと。 今、あるマーカーを購入した人数はtransactionデータをZDD::lcm()関数を使えば出る。特に、k人以上が購入したマーカーも出る。さらに言えば、複数のマーカーの組み合わせのうちk人以上が購入したものも出る。 例えば p1 = ZDD::lcm("FQ","trans.txt",100,"order.txt") とすれば100人以上が購入したマーカー組み合わせが列挙される。 今、ケースとコントロールとがna,no人ずつ居て、トータルn=na+noとしたとき、p個(p=1,2,...)のマーカーのすべてを購入した(x)か、そうでない(y)かで２ｘ２表を作るとすると、フェノタイプに関する周辺度数はna,no。すべて購入したか否かの周辺度数はZDD::lcm()が教えてくれる。上に挙げたコマンドでは100人以上が購入したマーカーセットを列挙するコマンド。

一方、周辺度数がna+no,x+yの２ｘ２表を考えたとき、その周辺度数を満足する条件でもっとも検定統計量が大きくなるときのp値(得られうる最小p値)というのは、周辺度数のみで決まる。

そうすると、na+noが与えられているときに、ある小さなp値を定めたとき、xの値がいくつ以上であればそれが起き得るかもわかる。そのようにしてxの値を定め、それによってマーカー組み合わせを取り出してみたときに、p値とマーカー組み合わせ数との関係はどうなるだろうか。

[ryamada22]

ZDD::lcm()のLCMはLinear-time Closed itemset Minerの頭文字。LCMについては　http://research.nii.ac.jp/~uno/papers/0409DM.pdf

変数順序付けの影響について：http://www.sofken.com/FIT2009/pdf/D/D_031.pdf

[ryamada22]

さて。マーカー組み合わせを網羅しつつ、周辺度数がまあまあ大きい場合を列挙できることになった。 同様にしてケースのみについても同様に、マーカー組み合わせ該当者数がまあまあ大きい場足を列挙することもできるだろう。

そうすると、マーカーの組み合わせを網羅しつつ、周辺度数がまあまあ大きく、かつ、ケースのそれ(２ｘ２表の[1,1]セルの値)がまあまあ大きい、というマーカー組み合わせ２ｘ２表が列挙できる。

それぞれの列挙は、ZDDによって、組み合わせが保持されている。

ZDD演算をすることで、周辺度数がまあまあ大きく、ケースの値が「小さい」組み合わせも列挙できる。

今、周辺度数の「まあまあ」の値と、[1,1]セルのまあまあの値の閾値として、検定p値がある値より小さくなる/大きくなるという閾値を設定すれば、「p値が十分小さくなるかもしれないマーカー組み合わせだけ」を列挙することができる。

それらについてp値を計算してやれば、「検定p値が十分小さくなる場合を列挙」できるのでは…

[ryamada22]

実験用Rソース n.sample <- 200 n.gene <- 1000 g <- matrix(sample(0:1,n.sample*n.gene,replace=TRUE,prob=c(0.7,0.3)),nrow=n.sample)

n.e <- 10 rs <- rpois(n.e,2) effects <- sample(1:n.e)

cmb <- list() for(i in 1:n.e){ cmb[[i]] <- sample(1:n.gene,rs[i]) }

scores <- rep(0,n.sample) for(i in 1:n.e){ if(length(cmb[[i]])!=0){ tmp <- matrix(g[,cmb[[i]]],nrow=n.sample) tmp2 <- apply(tmp,1,prod) tmp3 <- which(tmp2==1) if(length(tmp3)>0){ scores[tmp3] <- scores[tmp3] + effects[i] } } }

phenotype <- rep(0,n.sample) phenotype[which(scores > quantile(scores,0.5))] <- 1

g. <- list() g.diff <- list() icnt <- 1 for(i in 1:n.sample){ if(phenotype[i]==0){ g.[[i]] <- which(g[i,]==1)

g.diff[[i]] <- which(g[i,]==0)

}else{ g.[[i]] <- which(g[i,]==1) g.diff[[icnt]] <- g.[[i]] icnt <- icnt+1 } }

tab <- c() tab.diff <- c() for(i in 1:length(g.)){ ca= paste(g.[[i]], collapse=" ") tab <- append(tab,ca) } for(i in 1:length(g.diff)){ ca.diff <- paste(g.diff[[i]],collapse=" ") tab.diff <- append(tab.diff,ca.diff) } out <- file("trans.txt", "w") write.table(tab,out,row.names=F,col.names=F,quote=F) close(out) out <- file("trans_diff.txt","w") write.table(tab.diff,out,row.names=F,col.names=F,quote=F)

close(out)

order <- seq(1:n.gene) order <- paste(order,collapse=" ") out2 <- file("order.txt", "w") write.table(order,out2,row.names=F,col.names=F,quote=F)

close(out2)

実験用rubyコード $LOAD_PATH.push("/home/ryamada/.gem/ruby/2.2.0/gems/zdd-1.0.0-linux/lib/nysol/") require 'zdd'

p = ZDD::lcm("FQ","trans.txt",50,"order.txt") pdiff = ZDD::lcm("FQ","trans_diff.txt",30,"order.txt")

puts p.count puts pdiff.count

[skoyama427]

上記ソースは確認しておきます。

要件を満たしているかどうかわからないのですが、今までのところをアップしておきます vsop6_1.R -(1) vsop6.rb -(2) vsop6_2.R -(3) の順で実行すると、ケースで頻出するパターンに基づいて作った2x2表のp値の分布と、頻出パターン中に現れるマーカーのヒストグラムを出力します。バグがないことを祈ります。 未実装 to do ・頻出パターンの出現閾値を周辺度数から自動調整 ・コントロールに頻出するパターンの除去

初歩的な質問なのですが、多重検定の補正は今回のように「そうでない組み合わせも考え得たが、今回はこういう条件でこれだけしか考えなかった」という場合は検定しなかったが考え得た場合について補正はしなくても良いのでしょうか？

[ryamada22]

初歩的な質問なのですが、多重検定の補正は今回のように「そうでない組み合わせも考え得たが、今回はこういう条件でこれだけしか考えなかった」という場合は検定しなかったが考え得た場合について補正はしなくても良いのでしょうか？

これは初歩的な問題ではありません。哲学的ともいえる問題です。GWASでN個のSNPの検定をするつもりで１個ずつ検定をしたときに、最初の４個をやって、そのあと、検定をしなかったときの多重検定はどう考えるのか、とか、将来的にはNGSで全バリアントを確かめるつもりだけれど、今は予算が足りないのでSNPチップでデータを得てｐ値を出したが、多重検定はどうすればよいのか、とか。

[statgenetJimu]

SNPマーカーの他に、フェノタイプも「買い物アイテム」として登録し、ZDD::lcm()で頻出アイテムセットを抽出する。その後で、フェノタイプというアイテムでZDDの除算をすると、フェノタイプ1を持つことと併存する頻出アイテムセットのみが取り出せる。これが、アイテムセットマイニング with ZDD で言うところの、「フェノタイプ併存頻出マーカー探索」だろうと思います。 さらに、SNPマーカーの0/1を反転させると、フェノタイプと併存するアイテムセットに関する「裏～OR的な」情報抽出になるのではないか…(これはまだ考え中)、と思います。 問題は、それでどれくらい、本物らしいものが落とさずに取り出せるのかということなのかと思いますが、それをどう評価するのかを考える必要がありそうです。 特に、フェノタイプ＝１のサンプルが亜分類されて、ある亜集合ではあるアイテムセットによって説明され、別の亜集合では別のアイテムセットで説明され、という複合的な状況を取り出しやすそうなのが、このZDDアイテムセットマイニングっぽいので、そこを確認したいです：　http://www.kamishima.net/archive/freqpat.pdf