p.adjust()関数を使う

• 昨日の続き
• マイクロアレイデータをクラスタリングした
• フェノタイプと検定して、FDR補正してみる

# サンプル数
Ns <- 500
# マーカー数
Nm <- 1000
# サンプルのパターン数(群数)
Ns.pt <- 10
# マーカーのパターン数(群数)
Nm.pt <- 10
# サンプル・マーカーの存在位置を多次元空間酔歩の道として作る
trail <- matrix(rnorm(Ns.pt*Nm.pt),Ns.pt,Nm.pt)
trail <- apply(trail,2,cumsum)
# ３次元分だけ見てみよう
library(rgl)
plot3d(trail[,1:3])
matplot(trail,type="l")
# パターン数(群数)ごとにいくつのサンプル、いくつのマーカーを帰属させるかをランダムに決める
library(MCMCpack)
ps <- rdirichlet(1,rep(1,Ns.pt))
pm <- rdirichlet(1,rep(1,Nm.pt))
ss <- sample(1:Ns.pt,Ns,replace=TRUE,prob=ps)
sm <- sample(1:Nm.pt,Nm,replace=TRUE,prob=pm)
# 行数=サンプル数、列数=マーカー数の行列
M <- trail[ss,sm]
# 少し乱す
M <- jitter(M,1000)

• フェノタイプをMに依存して作る

# Mに依存したphenotypeを作る
# 標本ごとに全マーカーの値の和を取る
phenotype <- apply(M,1,sum)
# その値と平均との大小関係で0か1かを割り振る
phenotype <- (sign(phenotype - mean(phenotype)) + 1)/2
# 一部のサンプル(4/5)のphenotypeをランダムに振りなおす
phenotype[1:(4*Ns/5)] <- phenotype[sample(1:(4*Ns/5))]
#phenotype <- sample(1:2,Ns,replace=TRUE)
# ２群に属するサンプルIDをとりだす
group1 <- which(phenotype==0)
group2 <- which(phenotype==1)

• 各マーカー対フェノタイプでt検定する

# t.testする
p <- apply(M,2,function(x){t.test(x[group1],x[group2])[[3]]})
#p2 <- rep(0,length(M[,1]))
#for(i in 1:length(M[1,])){
#	p2[i] <- t.test(M[group1,i],M[group2,i])[[3]]
#}

• 出てきたp値の分布をチェックする
  • 小さ目のp値が多いことがわかる

# p値をソートしてプロットすると、p値の偏り具合がわかる
plot(sort(p))
#plot(p, p2)

• fdr補正する

# fdr法で補正する
fdr.p <- p.adjust(p,"fdr")

• 補正前後のp値を比較する
  • 補正されてp値が大きくなっている

# 補正前後のp値を比較する
plot(p,fdr.p)

• 補正の具合を図を使って調べる(こちらから)
  • 図、右端にごくわずか、補正後も有意な検定があることがわかる

# 図を使ってfdr補正の意味合いを見る
# 補正前のp値が低い方が補正後のp値も低い
# カットオフp値を、補正後p値が下回る比率が制御される

peven<-ppoints(length(p),a=0)

a=0.05
plot(peven,sort(peven,decreasing=TRUE),xlim=c(0,1),ylim=c(0,1),type="l")
par(new=TRUE)
#この線がFDR(BH法)の基準線。これより下なら採択。
plot(peven,a*((length(p)-(1:length(p))+1)/length(p)),xlim=c(0,1),ylim=c(0,1),type="l",col="red")
par(new=TRUE)
plot(peven,sort(fdr.p,decreasing=TRUE),xlim=c(0,1),ylim=c(0,1),type="l",col="blue")
abline(h=a)
#p.adjust()関数の返り値は基準直線としてaをいくつに設定したら交わるか、の値を返す。
#水平直線が0.05の高さになっていることがわかる
abline(v=peven[length(which(fdr.p>a))])