- 2015年度の修士向け講義は、「基礎の基礎」に戻ります(こちら)
- 日本人類遺伝学会の教育講演も、それに沿った話にします
- そのためのメモ
- こちらに図などがうまくいっていないかもしれないepub(文書の散逸防止のためにkindleにも置いています。以下のRmdファイルからhtml化・epub化は可能です→こちらを参考に)
この文章は、日本人類遺伝学会第59回大会、日本遺伝子診療学会第21回大会での教育講演『遺伝子多型のためのp値の話』をまとめたものです。
こんな2x2分割表があったとする。
```{r}
n <- 100
g <- sample(1:2,n,replace=TRUE)
g[1:(n/2)] <- sort(g[1:(n/2)],decreasing=TRUE)
p <- sample(0:1,n,replace=TRUE)
tab <- table(data.frame(p,g))
tab
```
カイ二乗検定をして、縦軸と横軸が無関係かどうかの検定をする。
```{r}
chisq.test(tab)
```
p値が帰ってきて、棄却するべきかどうかの情報が得られる。
p値が小さいほど、無関係という仮説（帰無仮説）を棄却することになる。
「基準を定めて」p < 0.05 なら棄却する、ということにすれば、棄却するかどうかの決断ができる。
「その基準を有意水準として、帰無仮説を棄却する」と説明し、「関連があった」と主張する。
その逆なら「関連がなかった」と主張する。
その心は、「帰無仮説が正しいと仮定したときに、このテーブルと同じかそれよりもっと偏ったテーブルを観察する確率はpくらい小さいから（pくらい大きいから）」。
オッズ比を計算する。
オッズ比は、相対危険度の推定値として使えるから。
オッズ比〜相対危険度が１であるときが、「無関係」に相当する。
推定値なので、「この値らしい」という点推定値があり、また、推定には幅があるから、だいたいこのくらいの範囲という意味で信頼区間がある（細かいことはここでは省略する）。
「だいたいこのくらいの範囲」の中に、「無関係の相対危険度〜１」が入っているかどうかで、およその帰無仮説を棄却するかどうかのめども立つ。
```{r}
library(vcd)
or <- oddsratio(tab,log=FALSE)
or
confint(or)
```
このようにデータが得られたら、たいていの場合、「検定」と「推定」をする。
「検定」するのは、「帰無仮説」を捨てるため。
「推定」するのは、相対危険度など、「因子の強さ」の程度に見当をつけるため。
たいていの場合は、「検定」をして、「関連あり」を選んだ後、「その関連の強さを数値として理解する」というプロセスをとる。
たくさんの2x2表を適当に作り、カイ二乗検定p値とオッズ比の95%信頼区間を出して、その関係をプロットしてみる。
横軸にp値、縦軸に95%信頼区間をとると、オッズ比の推定値とその信頼区間は１より大きい側になることもあれば、小さい側になることもあり、また、p値の大小の順番がオッズ比の信頼区間の上下限の値の大小の順番と完璧に一致するわけではないので、多少のでこぼこはあるが、総じて、p値の大小とオッズ比の信頼区間の大小には滑らかな関係があり、p値が0.05のあたりで、オッズ比の95%信頼区間に1が入るかどうかが分離できることが見て取れる。
```{r,warning=FALSE}
library(MCMCpack)
n <- 100
n.iter <- 1000
fracs <- rdirichlet(n.iter,rep(50,4))
ps <- rep(0,n.iter)
ors <- matrix(0,n.iter,2)
for(i in 1:n.iter){
tmp.tab <- matrix(rmultinom(1,n,fracs[i,]),2,2)
ps[i] <- chisq.test(tmp.tab)$p.value
ors[i,] <- confint(oddsratio(tmp.tab,log=FALSE))
}
plot(ps,ors[,1],pch=20,cex=0.5,xlim=c(0,0.2),xlab="p-value",ylab="95%CI",main="p値とオッズ比の95%信頼区間の関係")
points(ps,ors[,2],pch=20,cex=0.5)
for(i in 1:n.iter){
segments(ps[i],ors[i,1],ps[i],ors[i,2])
}
abline(h=1,col=3)
abline(v=0.05,col=3)
```
起きやすさを0から1までの数値で表したもの。
色々な検定方法などでは、異なる名前の統計量が計算されるけれど、最終的には、p値に変換することで、解釈がしやすくなる。
p値が小さければ小さいほど、観測データと同じかそれよりも起こりにくいデータは生じにくいことを意味する。
したがって、p値が小さいとき、その観測データを引き起こす前提として設定した仮説（帰無仮説）は信用しにくいことになる。
p値が小さければ小さいほど、「観測データと同じかそれよりも起こりにくいデータが得られる確率の総和」は小さくなる。
しかしながら、「観測データそのものが起きる確率」自体が小さいわけではない。
実際、p値は0から1のどの値も同程度に起きやすく起きにくい。
それを『帰無仮説が成り立つとき、p値は一様分布に従う』と言う。
p値が仮説を信じるかどうかにおいて、わかりやすいのは、この『一様分布に従う』という性質のお蔭。
一様分布とは、次のヒストグラムに示すように、等確率の分布のこと。
ヒストグラムを描いてみる。
```{r}
n <- 10000
ps <- runif(n)
hist(ps,xlab="p value")
```
分布が一様であることを示すには、小さい順に並べてプロットすることでも示せる。
```{r}
plot(sort(ps))
```
このように、「観察されたp値は一様分布である」ということを、「グラフにすると対角線上に並ぶ」ということを利用したのがQQプロット。
p値・一様分布の場合には、ただ、小さい順に並べてプロットすればよいけれど、それ以外(たとえばカイ二乗値)の場合には、小さい順に並べてプロットすると直線ではなくてカーブを描く。
「たくさんの観測されたカイ二乗値が、カイ二乗分布になるはず」と思っているときは、横軸に「カイ二乗分布からの値が同じ数だけランダムに得られたとして、それを小さい順に並べたときの理論的な値」を、縦軸に「観察された値」を取ってプロットしたものが、QQプロット。
対角線上に並んでいれば、「理論通り」と言える。
QQプロットは、「帰無仮説が成り立つときに、たくさんのp値が得られたとしたときに、得られるべきp値の理論値を横軸に、得られたp値の観測値を縦軸にしたプロット」。
観測が理論通りなら、QQプロットは対角線上、に一直線になる。
下の図はその様子。
```{r,warning=FALSE}
qqplot(qunif(ppoints(n)),ps,xlab="Theoretical",ylab="Observed",main="一様分布に従うp値のQQ plot")
```
SNVのほとんどは、2アレル(M,m)。
常染色体上のSNVは3ジェノタイプ(MM,Mm,mm)。
病気 vs. 非病気のとき、フェノタイプは0/1の2値。
```{r}
table2columns <- function(tab){
m <- matrix(0,0,2)
for(i in 1:length(tab[,1])){
for(j in 1:length(tab[1,])){
m <- rbind(m,cbind(rep(i,tab[i,j]),rep(j,tab[i,j])))
}
}
m
}
tab <- matrix(c(10,10,10,20,10,40),2,3)
tab
gp <- table2columns(tab)
p <- gp[,1] - 1
g <- gp[,2]
#gp
#n <- 1000
#a <- 0.3
#g.f <- c(a^2,2*a*(1-a),(1-a)^2)
#g <- sample(0:2, n, replace = TRUE,prob=g.f)
#p <- rep(0,n)
#pr <- c(0.35,0.4,0.45)
#for(i in 1:3){
# p[which(g==i-1)] <- sample(0:1,length(which(g==i-1)),replace=TRUE,prob=c(pr[i],1-pr[i]))
#}
#g[1:(n/8)] <- sort(g[1:(n/8)], decreasing = TRUE)
#p[1:(n/8)] <- sort(p[1:(n/8)], decreasing = TRUE)
#tab <- table(data.frame(p,g))
#tab
```
3つのジェノタイプはあるけれど、アレルがリスクを生み出しており、アレル本数が0,1,2本と増えるにしたがってリスクが増えると考えるのは、よい考え方だ。
その考え方に沿った検定がトレント検定とロジスティック回帰検定です。
まず、テーブルを作ってみる。
3ジェのタイプのコントロールの数は等しく、ケースの数が増えている。
```{r}
tab <- matrix(c(10,5,10,20,10,30),2,3)
#tab <- matrix(c(40,80,50,50,50,50),2,3)*3
tab
```
トレンド検定して、オッズ比も計算する。
```{r}
gp <- table2columns(tab)
p <- gp[,1] - 1
g <- gp[,2]
prop.trend.test(tab[1,],apply(tab,2,sum))$p.value
# オッズ比
or2vs1 <- oddsratio(tab[,1:2],log=FALSE)
or2vs1
#confint(or2vs1)
or3vs1 <- oddsratio(tab[,c(1,3)],log=FALSE)
or3vs1
```
オッズ比は、1.5、2と計算できましたが、『トレンド検定をするという気持ち』に照らすと、このオッズ比はちょっと違います。
トレンド検定をするというのは、線形回帰をしているのと同じです。
線形回帰では、ジェノタイプごとのケースの割合が直線に乗るように推定します。
```{r,warning=FALSE}
# 線形回帰と同じこと
lm.out <- lm(p~g)
lm.out$coefficients
a <- lm.out$coefficients
plot(matrix(c(1,0,2,0,3,0,1,1,2,1,3,1),byrow=TRUE,ncol=2),pch=20,cex=t(tab)*0.5,xlim = c(-2,6),ylim=c(-1,2),xlab="genotype",ylab="phenotype",main="トレンド検定と線形回帰")
abline(a[1],a[2],col=2)
p.prop <- a[1] + (1:3)*a[2]
points(cbind(1:3,p.prop),col=2,pch=20,cex=3)
abline(h=c(0,1))
abline(v=c(1,2,3))
```
このように、ジェノタイプごとのフェノタイプ率が推定できるので、基準とするジェノタイプを決めれば、それに対するオッズ比も計算できる。
トレンド検定・線形回帰では、このようにして計算したオッズ比を想定していたことになる。
トレンド検定のp値を用いて「関連あり」と信じるのであれば、トレンド検定が支持しているオッズ比は、1.5,2ではなく、線形回帰から計算されるオッズ比になる。
```{r}
# オッズとオッズ比
odds <- p.prop/(1-p.prop)
(or2vs1.2 <- odds[2]/odds[1])
(or3vs1.2 <- odds[3]/odds[1])
```
ロジスティック回帰もトレンド検定と同様によく用いられる。
やってみる。
p値と
```{r}
glm.out <- glm(p~g,family=binomial)
summary(glm.out)$coefficients[2,4]
```
2ジェノタイプのオッズ比が得られる。
```{r,warning=FALSE}
exp(summary(glm(p~g,family=binomial))$coefficients[2,1])^(1:2)
#summary(glm(p~g,family=binomial))$coefficients
```
ロジスティック回帰によって推定されるオッズ比の特徴は、ホモ接合体のそれがヘテロ接合体のそれの二乗になっていることである。
トレンド検定では直線を引いたが、ロジスティック回帰では、シグモイドカーブを引く。
```{r,warning=FALSE}
plot(matrix(c(1,0,2,0,3,0,1,1,2,1,3,1),byrow=TRUE,ncol=2),pch=20,cex=t(tab)*0.5,xlim = c(-2,6),ylim=c(-1,2),xlab="genotype",ylab="phenotype",main="ロジスティック回帰検定")
x <- seq(from=-2,to=6,length=100)
tmp <- exp(glm.out$coefficients[1]+glm.out$coefficients[2]*x)
y <- tmp/(1+tmp)
points(x,y,col=3,type="l")
abline(h=c(0,1))
abline(v=c(1,2,3))
```
トレンド検定とロジスティック回帰とを併せて描いてみる。
直線と曲線はそれなりに違うけれど、3ジェノタイプに対応するx軸上の位置はほとんど同じであることが分かる。
```{r,warning=FALSE}
plot(matrix(c(1,0,2,0,3,0,1,1,2,1,3,1),byrow=TRUE,ncol=2),pch=20,cex=t(tab)*0.5,xlim = c(-2,6),ylim=c(-1,2),xlab="genotype",ylab="phenotype",main="トレンド検定とロジスティック回帰検定")
abline(lm.out$coefficients[1],lm.out$coefficients[2],col=2)
points(x,y,col=3,type="l")
abline(h=c(0,1))
abline(v=c(1,2,3))
```
式に表して整理してみる。
トレンド検定では、
$$
\text{フェノタイプ1の割合} = a \times \text{ジェノタイプ} + b
$$
という式、
ロジスティック回帰では
$$
\frac{\text{フェノタイプ1の割合}}{\text{フェノタイプ0の割合}} = b\times e^{a \times \text{ジェタノタイプ}}
$$
となるように線を引いている。
```{r,warning=FALSE}
n.iter <- 100
p.prop <- p.log <- rep(0,n.iter)
for(i in 1:n.iter){
a <- runif(1,0.2,0.8)
g.f <- c(a^2,2*a*(1-a),(1-a)^2)
g <- sample(0:2, n, replace = TRUE,prob=g.f)
p <- sample(0:1, n, replace = TRUE)
tab <- table(data.frame(p,g))
p.prop[i] <- prop.trend.test(tab[1,],apply(tab,2,sum))$p.value
p.log[i] <- summary(glm(p~g,family=binomial))$coefficients[2,4]
}
plot(log(p.prop),log(p.log),xlab="トレンド",ylab="ロジスティック",main="トレンドとロジスティックのlog(p)")
abline(0,1,col=2)
```
トレンド検定は
> 相加モデルに相当する
> 計算が簡単
> 正確検定もできる(低アレル頻度・少サンプル数の場合の対処法がある)
> 個人別のデータが不要
> 「線形回帰」に基づくオッズ比が結果に付いていないことが多い
> 共変量を組み込めない(線形回帰にすればよいが、それならロジスティック回帰をすればよい)
ロジスティック回帰検定は
> 相乗モデルに相当する
> 計算が面倒ではある(けれど計算機がやってくれるので問題はないが、検定個数が大量になるとそれなりに影響してくる)
> オッズ比が結果についてくることが多い(ただし、係数は対数で返ってくることが通例)
> 年齢・性別などの共変量を組み込みやすい
トレンド検定では、ジェノタイプ別の発病率の上昇分がリスクアレルの本数に比例していると仮定していたので相加モデルと呼ばれ、ロジスティック回帰では、
ジェノタイプ別のオッズ比がリスクアレルの本数のべき乗になっていると仮定したので相乗モデルと呼ばれます。
これらはいずれも、リスクアレル本数が１のとき(ヘテロ接合体)のリスクを
二つのホモ接合体のリスクの中間的な値とするモデルです。
これ以外に、単一遺伝子病などに認められてきた遺伝形式に優性形式と劣性形式があります。
優性形式、劣性形式に合わせた検定というものがあります。
優性形式では、リスクアレルを0本持つか、1本以上持つかの2つに分けて検定しますし、劣性形式では、その逆で、2本持つか、そうでないかに分けて検定します。
実際にはリスクアレルがSNVのどちらのアレルなのかの予想がついていないことがほとんどなので、0 vs.(1,2) と (0,1) vs. 2の両方の検定を実施して、結果として、優性形式と劣性形式とについて検定することになります。
まず初めに考えることは、「優性・劣性形式の検定をするのかしないのか」です。
その次に考えることは、「優性・劣性形式の検定をしたとして、結果をどう解釈するか」「優性・劣性形式の検定をしなかったとして、結果をどう解釈するか」です。
することのメリットとデメリット
* 優性・劣性形式に照らした結果がわかる
* 1つの2x3表に複数の検定をすると、複数のp値が得られる。複数のp値が出たら、そのp値はには「補正」をしないといけない
しないことのメリットとデメリット
* 1つのp値しか出ていなければ、複数のp値の補正について悩む必要はない
* 優性・劣性形式に照らしての判断になっていない
相加・相乗モデルでの検定から得られる1個のp値で、「関連ありなし」を判断することになるので、
本当は「優性、または、劣性」な影響があるとすると、偽陰性が増える。
ジェノタイプ別のオッズ比を計算して、「優性形式っぽい」ときには、「優性なのかも」と思い、「劣性形式っぽい」ときには、「劣性なのかも」と思う、という考え方もある。
ただし、トレンド検定・ロジスティック回帰をしたときには、「そのモデルに基づいたジェノタイプ別のオッズ比は『すでに出ている』」わけで、ジェノタイプ別にオッズ比を出しなおして、遺伝形式を問題にしないほうがよい、という考え方もある。
相加・相乗モデルでの検定で、ひとまず「関連があるかないか」の選別をして、そのうえで、「相加・相乗モデル」と「優性（劣性）モデル」のどちらが、よりもっともらしいのか、を数字に変えることは可能だけれど、通常、そのようなことはやらない。
では、少し戻って、「優性形式が真」のときに、「相加モデルだけ」で検定すると、どれくらい偽陰性になるのかを見てみることにする。
たとえば、頻度が0.3のアレルがリスクアレルであって、優性形式でフェノタイプに影響しているとする。
今、リスクアレルを1本または2本持つと、リスクが1.2倍になるとし、
ケース・コントロール、ともに1000人ずつで、関連を調べるとする。
相加モデルで検定するか、優性モデルで検定するかで、p値はどのようになるか、その違いが問題となる。
優性モデルで検定すると有意になるのは30%程度なのに対して、相加モデルで検定すると、23%程度となる。
このパワーの差が、真のモデルに沿って検定するかどうかの違いである。
どちらの検定でもp値が0.01未満になる場合、どちらの検定でもならない場合、どちらか片方の検定でのみp値が0.01未満になる場合の4通りのうち、優性モデルで検定すれば、有意になるのに、相加モデルで検定したがために有意にならなかった場合が問題である。
図ではそれを赤点で表している。
```{r,warning=FALSE}
a <- 0.3
r <- c(1.2,1.2,1)
n.case <- n.cont <- 1000
g.cont <- c(a^2,2*a*(1-a),(1-a)^2)
g.case <- g.cont * r
g.case <- g.case/sum(g.case)
n.iter <- 1000
v.case <- rmultinom(n.iter,n.case,g.case)
v.cont <- rmultinom(n.iter,n.cont,g.cont)
p.dom <- p.add <- rep(0,n.iter)
for(i in 1:n.iter){
p.dom[i] <- prop.trend.test(v.case[,i],v.case[,i]+v.cont[,i],score=c(2,2,1))$p.value
p.add[i] <- prop.trend.test(v.case[,i],v.case[,i]+v.cont[,i])$p.value
}
col <- rep(1,n.iter)
col[which(p.dom < 0.01 & p.add >= 0.01)] <- 2
plot(log10(p.add),log10(p.dom),pch=20,cex=0.1,col=col)
abline(0,1,col=3)
abline(h=-2,col=3)
abline(v=-2,col=3)
length(which(p.add<0.01))/n.iter
length(which(p.dom<0.01))/n.iter
```
『相加モデルと優性モデルの両方で検定してやればよいのでは…!?』
先ほどの例で、相加モデルと優性モデルとの両方で検定をして、どちらでもよいから、有意なp値が得られたら、「関連あり」とみなすことにすれば、パワーが上がってよいだろう、ということになる。
図で示せば、色のついたところのすべてを関連ありとする、ということである。
```{r}
col <- rep(1,n.iter)
col[which(p.dom < 0.01 | p.add < 0.01)] <- 4
col[which(p.dom < 0.01 & p.add >= 0.01)] <- 2
plot(log10(p.add),log10(p.dom),pch=20,cex=0.1,col=col)
abline(0,1,col=3)
abline(h=-2,col=3)
abline(v=-2,col=3)
```
しかし、残念なことに、今度は、パワーではなくて、偽陽性の方に注意しなくてはいけなくなる。
リスクがないSNVで同じことをやってみる。
2つの検定のどちらかでp < 0.01 になったら、「関連あり」とみなす、という作戦なので、1.4% で「関連あり」となってしまうことがわかる。
1% でだけ「関連あり」となることを許容しているのに、1.4%も「関連あり」としてしまうのは、偽陽性が1.4倍に増えるから問題である。
```{r,warning=FALSE}
a <- 0.3
r <- c(1,1,1)
n.case <- n.cont <- 1000
g.cont <- c(a^2,2*a*(1-a),(1-a)^2)
g.case <- g.cont * r
g.case <- g.case/sum(g.case)
n.iter <- 1000
v.case <- rmultinom(n.iter,n.case,g.case)
v.cont <- rmultinom(n.iter,n.cont,g.cont)
p.dom.2 <- p.add.2 <- rep(0,n.iter)
for(i in 1:n.iter){
p.dom.2[i] <- prop.trend.test(v.case[,i],v.case[,i]+v.cont[,i],score=c(2,2,1))$p.value
p.add.2[i] <- prop.trend.test(v.case[,i],v.case[,i]+v.cont[,i])$p.value
}
col.2 <- rep(1,n.iter)
col.2[which(p.dom.2 < 0.01 | p.add.2 < 0.01)] <- 4
plot(log10(p.add.2),log10(p.dom.2),pch=20,cex=0.1,col=col.2)
abline(0,1,col=3)
abline(h=-2,col=3)
abline(v=-2,col=3)
length(which(p.add.2<0.01 | p.dom.2 < 0.01))/n.iter
```
優性モデルにもパワーを出しつつ、偽陽性を増やさないようにしたい。
そのためには、優性モデル検定と相加モデル検定の両方を実施して、そのどちらかで小さいp値が出たら「関連あり」とすることにしつつ、その「関連あり」とみなすための判定基準のp値を厳しめにする、という作戦をとるとよい。
図で示せば、基準となる線を引きなおして、偽陽性が1%になるようなところを見つければよい。
```{r}
par(mfcol=c(1,2))
plot(log10(p.add),log10(p.dom),pch=20,cex=0.1,col=col)
abline(0,1,col=3)
abline(h=-2,col=3)
abline(v=-2,col=3)
col.3 <- col.2
col.3[which(p.add.2 < -2.2 | p.dom.2 < -2.2)] <- 5
plot(log10(p.add.2),log10(p.dom.2),pch=20,cex=0.1,col=col.3)
abline(0,1,col=3)
abline(h=-2,col=3)
abline(v=-2,col=3)
abline(h=-2.2,col=4,lw=2)
abline(v=-2.2,col=4,lw=2)
par(mfcol=c(1,1))
```
力技で見つけてみる。
```{r,warning=FALSE}
p.check <- seq(from=0,to=1,length=1000)
f <- rep(0,length(p.check))
for(i in 1:length(f)){
f[i] <- length(which(p.add.2 < p.check[i] | p.dom.2 < p.check[i]))/length(p.add.2)
}
plot(log10(p.check),log10(f),xlab="基準の値(対数)",ylab="有意の割合(対数)")
abline(h=-2,col=2)
clst <- which(abs(log10(f)+2) == min(abs(log10(f)+2)))[1]
abline(v=log10(p.check[clst]),col=2)
print(p.check[clst])
```
このように、「帰無仮説」のときに、複数のp値がどういう組み合わせで出るのかがわかっていれば、「関連あり」とみなすp値の基準を決めて、偽陽性が0.01になるようにできることがわかります。
問題は、「複数のp値がどういう組み合わせで出るのか」がわかっていないことです。
2つの検定のp値を縦軸と横軸に取ります。
2つの検定のp値がお互いに関係がないものとします。
こんな図になります。
優性モデルと相加モデルのときとは違います。
p=0.01はわかりにくいのでp=0.1で考えてみます。
```{r,warning=FALSE}
par(mfcol=c(1,2))
n <- 1000
plot(p.add.2[1:n],p.dom.2[1:n],pch=20,main="優性・相加")
p1 <- runif(n)
p2 <- runif(n)
plot(p1,p2,pch=20,main="2つの独立な検定")
par(mfcol=c(1,1))
```
```{r,warning=FALSE}
par(mfcol=c(1,2))
n <- 1000
col <- rep(1,n)
col[which(p.add.2[1:n] < 0.1)] <- 2
col[which(p.dom.2[1:n] < 0.1)] <- 3
col[which(p.add.2[1:n] < 0.1 & p.dom.2[1:n] < 0.1)] <- 4
plot(p.add.2[1:n],p.dom.2[1:n],pch=20,col=col,main="優性・相加")
p1 <- runif(n)
p2 <- runif(n)
abline(v=0.1,col=2)
abline(h=0.1,col=2)
col2 <- rep(1,n)
col2[which(p1 < 0.1)] <- 2
col2[which(p2 < 0.1)] <- 3
col2[which(p1 < 0.1 & p2 < 0.1)] <- 4
plot(p1,p2,pch=20,col=col2,main="2つの独立な検定")
abline(v=0.1,col=2)
abline(h=0.1,col=2)
#plot(log10(p1),log10(p2),pch=20)
par(mfcol=c(1,1))
```
優性・相加の方が、独立な2つの場合に比べて
> 黒は多い
> 青も多い
> 赤と緑は少ない
```{r,warning=FALSE}
print(c("黒","赤","緑","青"))
table(col)
table(col2)
```
2つの検定を実施したら、p値を2倍して解釈する。n個の検定を実施したら、p値をn倍して解釈する。
黒が90%くらいになるとよいのですが、90%より多いです。
「関連あり」は少なすぎるということです。
ボンフェロニ法は「関連を少なくし過ぎる」ので「保守的すぎる」と言います。
「2つの独立な検定」の場合で言えば、黒い点がある部分の面積が、$$0.9^2=0.81$$と、0.9より小さすぎたので、$$0.95^2=0.9025$$にしてみたら、今度は、0.9より大きくなってしまった、ということです。
$$\sqrt(0.9)^2 = 0.9 = (0.9487)^2$$にすると、ちょうどよいはずなので、そうすることもできます。
Sidakの補正方法です。
こうすれば、2つの検定が独立であれば、この方法でうまくいきます。
どちらか片方でp<0.1であるものを「関連あり」とすると、10%だとよいところが15%を「関連あり」とみなしてしまって、多すぎます。
ボンフェロニ法にすると、7.7%まで減ってしまって、保守的にすぎます。
2つの検定が独立だったときは、Sidakはよかったのですが、優性と相加のように相互に近い関係だろ、Sidak法でもボンフェロニ法とほとんど変わらず(7.7%)、まだまだ保守的です。
```{r,warning=FALSE}
par(mfcol=c(1,2))
n <- 1000
col <- rep(1,n)
col[which(p.add.2[1:n] < 0.05)] <- 2
col[which(p.dom.2[1:n] < 0.05)] <- 3
col[which(p.add.2[1:n] < 0.05 & p.dom.2[1:n] < 0.05)] <- 4
plot(p.add.2[1:n],p.dom.2[1:n],pch=20,col=col,main="優性・相加")
p1 <- runif(n)
p2 <- runif(n)
abline(v=0.1,col=2)
abline(h=0.1,col=2)
abline(v=0.05,col=3)
abline(h=0.05,col=3)
col2 <- rep(1,n)
col2[which(p1 < 0.05)] <- 2
col2[which(p2 < 0.05)] <- 3
col2[which(p1 < 0.05 & p2 < 0.05)] <- 4
plot(p1,p2,pch=20,col=col2,main="2つの独立な検定")
abline(v=0.1,col=2)
abline(h=0.1,col=2)
abline(v=0.05,col=3)
abline(h=0.05,col=3)
#plot(log10(p1),log10(p2),pch=20)
par(mfcol=c(1,1))
table(col)
table(col2)
```
```{r}
col.nominal <- rep(1,n)
col.nominal[which(p.add.2 < 0.1)] <- 2
col.nominal[which(p.dom.2 < 0.1)] <- 3
col.nominal[which(p.add.2 < 0.1 & p.dom.2 < 0.1)] <- 4
col.b <- rep(1,n)
col.b[which(p.add.2 < 0.05)] <- 2
col.b[which(p.dom.2 < 0.05)] <- 3
col.b[which(p.add.2 < 0.05 & p.dom.2 < 0.05)] <- 4
col.s <- rep(1,n)
col.s[which(p.add.2 < 1-sqrt(0.9))] <- 2
col.s[which(p.dom.2 < 1-sqrt(0.9))] <- 3
col.s[which(p.add.2 < 1-sqrt(0.9) & p.dom.2 < 1-sqrt(0.9))] <- 4
table(col.nominal)
table(col.b)
table(col.s)
```
分布がわかれば、補正の仕方はわかるということはわかりました。
実際には分布がわからないから困っているわけで、分布がわからないときにどうするかが問題です。
方法は2つあります。
> 分布がわからないままに、補正する
> わからない分布を調べてから、それに基づいて補正する
の2つです。
補正をしないままに「関連あり」と主張すると、「それは言い過ぎ」と反論されます。
逆に、ボンフェロニ法・Sidak法で補正したp値を使って「関連あり」と主張するのは、2つの検定が独立でも「言い過ぎ」ではないし、まして、2つの検定が独立でないときには、「もっと安心」なので、「関連あり」と言ってもよい、と言うことになります。
逆に「関連なし〜関連ないと信じても悪くない」と主張する場合には、ボンフェロニ法・Sidak法で補正したp値を使った場合には、「それは言い過ぎ」ということになりますし、補正しないままのp値で判断するなら、「関連なし」と主張する場合には、「安心」ということになります。
たいていの場合は「関連あり」と主張するので、「ボンフェロニ法・Sidak法で補正しているから安心」ということになります。
このように「言い過ぎ」はダメで「安心であればＯＫ」ということを「保守的な判断」と言いますが、
『こうだったら、論文がアクセプトになって好都合なのに』と思ったら、『自分に不都合に』判断することです。『ストイックに・禁欲的に』ですね。
わからないなら、調べればよいです。
たとえば、こんな2x3表が得られたとします。
もちろん、この表に対して優性検定・相加検定の実施は可能です。
今、知りたいのは、「関連がない」との仮定の下で、ケースとコントロールを、それぞれ55人、30人のサンプリングして、『このSNV』のタイピングをした結果の表をたくさん作ること、そして、その表に優性・相加の検定を実施してp値を求めることです。
そうは言っても、『このSNV』の集団でのジェノタイプ頻度の正確なところはわかっていませんから、仮想の表をたくさん作るのは難しいです。
そういうときは、分割表の周辺度数を固定して、周辺度数が同じ表をたくさん作るという方法もあります。
いわゆる、分割表の正確確率検定で、『周辺度数を等しくする表のすべての生起確率を計算』して、検定p値を決めるのと同じ方法です。
可能なすべての表を全部列挙してその生起確率を計算することは、できなくはないですが、結構、大変です。
その大変さは、表の自由度に応じてどんどん大きくなります。
2x3表の自由度は2なので、それほど大変ではないので、正確に確率を計算することはできますが、一般的なマルチプルテスティング対策としては、現実的ではありません。
そこで乱数を発生させて、正確確率法と同じような分布を得る方法がランダマイゼーション法・(乱数を使った)パーミュテーション法です。
```{r}
tab <- matrix(c(10,5,10,20,10,30),2,3)
tab
```
ランダマイゼーション法・パーミュテーション法は簡単です。
今、2x3表では、サンプル全員のフェノタイプとジェノタイプが観察されていますが、ランダマイゼーション法では、誰にどのフェノタイプが観測されたかを、『なかったこと』にして、ランダムに割り付けなおします。
そうすると、フェノタイプが観察と変わりますから、ジェノタイプと突き合わせてできる2x3表が変わります。
変わりますが、ケース・コントロールの人数、ジェノタイプ別の人数は変わりません。
周辺度数は変わっていないということです。
この割り付けを$$n!$$通り、すべて実施すれば、正確確率法と同じになりますし、$$n!$$は多すぎるから、適当な値でやめにすれば、適当な数の表がランダムに作れたことになります。
作った表で優性・相加の検定をしてp値を出してやれば、目的は達成されます。
分割表の周辺度数が等しいという条件の下、ある分割表を観察する確率を正確に計算したり、それに準ずる方法としてランダマイゼーション・パーミュテーション法を用いることで、互いに独立でない複数の検定のp値のマルチプルテスティング補正ができることがわかりました。
ここで、正確確率検定に話を戻します。
分割表があったときに、トレンド検定を実施することもできますが、正確確率検定を実施することもできるからです。
2x3表だと話が少し複雑になるので、2x2表に戻って考えることにします。
>「分割表のセルの期待値が小さめのときには、表にカイ二乗検定をするのではなく、フィッシャーの正確確率検定をする方がよい」
と言われています。
どうしてか？
> セルの期待値が小さめのときには、カイ二乗検定のp値は『不正確』だから
> 正確確率検定の方が『保守的』だから
例を見ておきます。
```{r}
tab <- matrix(c(18,2,4,10),2,2)
tab
chisq.test(tab,correct=FALSE)$p.value
fisher.test(tab)$p.value
```
これらは正しいので、表を観察したら、この原則に従うのが良いのです。
良いのですが、GWASなどでは、山ほど検定を行います。そうすると、この原則に従うとどうなるかというと…
帰無仮説の棄却水準を0.05として、割合が3%の因子について100人 vs. 100人でたくさんの検定をするとすると、
実際にp<0.05となる割合は、普通のカイ二乗検定ならは0.05くらいですが、フィッシャーの正確確率検定をすると0.012くらいになります。
なるべく偽陽性を出さないためには、よいことですが、パワーのことを考えると、随分、かい離した値と言えます。これは、表の期待値が小さめのときだけに正確確率検定を適用した場合でも、ほとんど変わりません。
また、p値の分布も、『一様』からは程遠いです。
このような条件でのp値はp<0.05の割合が少ないだけでなく、
p値の平均も随分と違います。
カイ二乗検定ならばp値の平均は0.5程度となり、一様分布らしさがありますが、正確検定を用いると、平均が0.7くらいになります。
たくさんの検定の結果からp値の分布と取り出し、それが一様分布になっているとよい、という考えの下、たくさんのp値の処理をするようなときに、この正確検定の『保守性』が邪魔になることもあります。
正確検定ならいつもよい、というわけではない、という話としてとらえておくことにします。
```{r,warning=FALSE}
n.iter <- 100
n.ca <- 100
n.co <- 100
a.ca <- 0.03
a.co <- 0.03
ca <- rmultinom(n.iter,n.ca,prob=c(a.ca,1-a.ca))
co <- rmultinom(n.iter,n.co,prob=c(a.co,1-a.co))
p.chi <- p.fi <- p.select <- min.exp <- rep(0,n.iter)
for(i in 1:n.iter){
tab <- matrix(c(ca[,i],co[,i]),2,2)
min.exp[i] <- min(matrix(apply(tab,1,sum),ncol=1) %*% matrix(apply(tab,2,sum),nrow=1)/sum(tab))
p.chi[i] <- chisq.test(tab,correct=FALSE)$p.value
p.fi[i] <- fisher.test(tab)$p.value
p.select[i] <- p.chi[i]
if(min.exp[i] < 5){
p.select[i] <- p.fi[i]
}
}
length(which(p.chi < 0.05))/n.iter
length(which(p.fi < 0.05))/n.iter
length(which(p.select < 0.05))/n.iter
```
```{r,warning=FALSE}
par(mfcol=c(1,2))
hist(p.chi,main="カイ二乗分布のp値")
hist(p.select,main="期待値が小さい時だけ正確検定したときのp値")
par(mfcol=c(1,1))
mean(p.chi[which(!is.na(p.chi))])
mean(p.fi)
mean(p.select)
```
一つのSNVに１つの検定をするか複数の検定をするか、という話をしてきました。
一つの遺伝子に複数のSNVがあり、その遺伝子とフェノタイプとに関連があるのかを調べることと似ています。
1つのSNV ------------1つの遺伝子
複数の検定-----------複数のSNVに対してい実施する複数の検定
SNVを2個観察すると、組み合わせジェノタイプは9種類になるので、
結局2x9表を観察することになります。
今、2個のSNVについて、それぞれ検定をしてフェノタイプとの関係を調べるとします。
よくやることです。
2x3表の3ジェノタイプに対して、相加モデルであったなら、1,2,3という重みをつけ、優性モデルであったなら、0,1,1という重みを付けたわけですが、
9ジェノタイプがあって、片方のSNVについて相加モデルで検定するというのは
```{r}
(m33 <- matrix(c(1,2,3,1,2,3,1,2,3),3,3))
```
という重みを付けることだし、もう一つのSNVについて相加モデルで検定するというのは
```{r}
t(m33)
```
という重みを付けることですから、
2x9表に2検定のマルチプルテスティングをしていることになって、先ほどの議論と同じです。
もし、2x9表の2つのSNVのそれぞれに相加モデルと優性モデルとをそれぞれやるのであれば、次のような4つの重みづけでマルチプルテスティングをしていることになります。
```{r}
m33
t(m33)
(m33.2 <- matrix(c(0,1,1,0,1,1,0,1,1),3,3))
t(m33.2)
```
検定間が独立であれば、Sidak法が有用で、そうでなければ、補正をしないなら、『保守的に』、補正をするなら『正確に確率を計算』するか『ランダマイゼーション法・パーミュテーション法』が有用なのでした。
かなりの強さのLDのある領域にn.snv=20個のSNVがあるとして、100人 vs. 100人のケースコントロール関連検定を実施するとします。
個々のSNVについて、相加モデル検定をすることにすれば、n.snv=20個の検定となります。
今、有意水準0.01での棄却をするには、20個のp値の中の最小のp値の分布を調べて、ある確率0.01未満でしか観測されないくらい小さいp値とはどれほどかを調べることになります。
LDのプロットではこれくらい強い領域であるとすると、それは約0.001と10分の1くらいに小さくする必要があることがランダマイゼーション法・パーミュテーション法から示されます。
```{r}
n.snv <- 20 # SNVの数
n.h <- 10 # そこにあるハプロタイプ種類数を限定的にしてLDを生じさせる
# 適当にハプロタイプ頻度を決めて、HWE仮定の下、ジェノタイプデータを作る
H <- matrix(sample(0:1,n.snv*n.h,replace=TRUE),ncol=n.snv)
library(MCMCpack)
library(genetics)
h.freq <- rdirichlet(1,rep(1,n.h))
# ケース・コントロール人数
n.case <- 100
n.cont <- 100
h.case <- matrix(sample(1:n.h,2*n.case,replace=TRUE,prob=h.freq),ncol=2)
h.cont <- matrix(sample(1:n.h,2*n.cont,replace=TRUE,prob=h.freq),ncol=2)
H.case.1 <- H[h.case[,1],]
H.case.2 <- H[h.case[,2],]
H.cont.1 <- H[h.cont[,1],]
H.cont.2 <- H[h.cont[,2],]
G.case <- H.case.1+H.case.2
G.cont <- H.cont.1+H.cont.2
G.case.cont <- rbind(G.case,G.cont)
P <- c(rep(1,n.case),rep(2,n.cont))
```
```{r}
# geneticsライブラリを使って、LDプロットを作る(そのためにデータ変換をする)
alleles <- c("A","T")
g.case <- g.cont <- g.case.cont <- list()
for(i in 1:n.snv){
g.case[[i]] <- as.genotype.allele.count(G.case[,i],alleles)
g.cont[[i]] <- as.genotype.allele.count(G.cont[,i],alleles)
g.case.cont[[i]] <- as.genotype.allele.count(c(G.case[,i],G.cont[,i]),alleles)
}
data <- data.frame(g.case.cont)
data <- makeGenotypes(data)
ldt <- LD(data)
LDtable(ldt)
```
フェノタイプレコードをランダムにシャッフルして、p値分布を作る
```{r,warning=FALSE}
n.iter <- 500
ps <- matrix(0,n.iter,n.snv)
for(i in 1:n.iter){
P.sh <- sample(P)
for(j in 1:n.snv){
tab <- table(data.frame(P.sh,G.case.cont[,j]))
ps[i,j] <- chisq.test(tab)$p.value
}
}
# 最小p値を取り出し、それをalpha水準に照らして、閾値を決める
hist(apply(ps,1,min))
alpha <- seq(from=0,to=1,by=0.01)
alpha <- alpha[-1]
q <- quantile(apply(ps,1,min),alpha)
correct.rate <- alpha/q
# alpha=0.01に相当するp値
q[1]
# 補正率は、p値が大きいほど小さくなる
plot(alpha,correct.rate)
# 独立検定何個分に相当するか。これを実際に実施した検定数と比較してみる
print(correct.rate[1])
```
ここでは、ある特定のLD領域にある複数のSNVのことを扱いましたが、これをゲノムワイドに広げても考え方は同じです。
50万から100万のSNVについて1SNVあたり1検定をするとして、有意水準を0.01にするとなると
ボンフェロニ補正するなら
$$
0.01 \times \frac{1}{10^6} = \frac{1}{10^8}
$$
くらい小さいp値が出て初めて、「帰無仮説が正しいとしたと仮定したら、珍しい」と言える、ことになりますが、
LDもあるし、少し割り引いて、$$10^{-7} - 10^{-8}$$ くらいでどう？
というのが、ゲノムワイドな補正p値の目安、と言うことです。
SNVの数がを十倍の一千万個に増やした場合には、
> 百万個の場合と同じ領域をカバーしながら、個数だけ増やしたのなら、マルチプルテスティングで目指すべきp値はそれほど下げなくてもよい
> 百万個の場合にカバーしなかった場所をカバーすることにしたのなら、マルチプルテスティングで目指すべきp値は増やした分だけ下げる(十分の１にする)ことになり
ます。
連鎖不平衡にある複数のSNVを用いた解析におけるマルチプルテスティングについて述べてきましたが、連鎖不平衡を用いた関連解析について少し確認しておくことにします。
連鎖不平衡解析では、「真の遺伝因子」があって、その近傍に「真の遺伝因子と連鎖不平衡関係にある代用マーカー」があるので、代用マーカーを調べることで、「真の遺伝子の位置」の見当をつけます。
連鎖不平衡とは、A/aというバリアントとB/bというバリアントがあったときに、A-B,A-b,a-B,a-bの4つのペアの割合について、2x2分割表を作り、それについての独立性検定をするようなものです。サンプル総数が1の2x2表のカイ二乗値の最大値は1で、最小値は0となり、それが連鎖不平衡係数$$r^$$に相当します。
念のためやっておきます。
```{r,warning=FALSE}
# 2SNVのアレル頻度と、SNV間のr(r^2の平方根、ただし正負あり)を与えて
# 4ハプロタイプ頻度を2x2行列で返す
my.2SNV.hap <- function(pA,pB,r){
pH.e <- outer(c(pA,1-pA),c(pB,1-pB),"*")
r.squared <- r^2
pH.d <- pH.e + sqrt(prod(c(pA,1-pA,pB,1-pB)))*c(1,-1,-1,1)*r
return(pH.d)
}
pA <- 0.8
pB <- 0.8
r <- 0.9
r.squared <- r^2
pH.d <- my.2SNV.hap(pA,pB,r)
pH.d
# カイ二乗テストをすると、そのカイ二乗値はr^2
r^2
chisq.test(pH.d,correct=FALSE)$statistic
```
では、2SNVがあって、その片方が真のリスクを有しているとして、真のリスクSNVとLD関係にあるSNVとで相加モデル検定をしたときに、p値がどんな風になるかを見てみます。
低頻度のアレルをリスクアレルとし、リスクアレルの保有コピー数が0,1,2のときに有病率を6%,8%,10%として見ます。
このようなときに、ケース・コントロール1000人数ずつで解析をすると、
2SNVのp値はいずれも低めに出て、両者には相関があります。
真のリスクSNVのp値の方が少し低めです。
```{r,warning=FALSE}
my.LD.sim <- function(h.freq,r.gen.1,n.case=1000,n.cont=1000,n.pop=10^6,n.iter=1000){
n.snv <- 2 # SNVの数
n.h <- 4
H <- matrix(c(0,0,0,1,1,0,1,1),byrow=TRUE,ncol=n.snv)
h.pop <- matrix(sample(1:n.h,2*n.pop,replace=TRUE,prob=h.freq),ncol=2)
H.pop.1 <- H[h.pop[,1],]
H.pop.2 <- H[h.pop[,2],]
G.pop <- H.pop.1+H.pop.2
R.pop <- G.pop[,1]
P.pop <- rep(0,n.pop)
for(i in 0:2){
P.pop[which(R.pop==i)] <- sample(0:1,replace=TRUE,length(which(R.pop==i)),prob=c(1-r.gen.1[i+1],r.gen.1[i+1]))
}
p1 <- p2 <- rep(0,n.iter)
for(i in 1:n.iter){
sample.case <- sample(which(P.pop==1),n.case)
sample.cont <- sample(which(P.pop==0),n.cont)
tmp.df1 <- data.frame(g=G.pop[c(sample.case,sample.cont),1],p=P.pop[c(sample.case,sample.cont)])
tmp.df2 <- data.frame(g=G.pop[c(sample.case,sample.cont),2],p=P.pop[c(sample.case,sample.cont)])
tab1 <- table(tmp.df1)
tab2 <- table(tmp.df2)
p1[i] <- prop.trend.test(tab1[,1],apply(tab1,1,sum))$p.value
p2[i] <- prop.trend.test(tab2[,1],apply(tab2,1,sum))$p.value
}
return(list(p1=p1,p2=p2))
}
pA <- 0.8
pB <- 0.8
r <- 0.9
r.squared <- r^2
pH.d <- my.2SNV.hap(pA,pB,r)
r.a <- -0.02
r.b <- 0.1
r.gen.1 <- (0:2)*r.a+r.b
r.gen.1
p.out <- my.LD.sim(h.freq=c(pH.d),r.gen.1 =r.gen.1)
par(mfcol=c(2,2))
plot(p.out$p1,p.out$p2)
plot(log10(p.out$p1),log10(p.out$p2))
h <- hist(log10(c(p.out$p1,p.out$p2)),col=gray(0.99),border=FALSE)
hist(log10(p.out$p1),density=17,add=TRUE,breaks=h$breaks)
hist(log10(p.out$p2),col=2,density=20,add=TRUE,breaks=h$breaks)
par(mfcol=c(1,1))
cor(p.out$p1,p.out$p2)
print(r.squared)
mean(p.out$p1)
mean(p.out$p2)
```
では同じ条件で、$$r^2$$だけを変化させるとどうなるかを見てみます。
2SNVのp値の相関が$$r^2$$に応じて変化します。相関係数と$$r^2$$とが同じ値です。
また、真のリスクSNVのp値がLD関係にあるSNVのp値以下になる頻度も$$r^2$$に応じて高くなることがわかります。$$r^2=0.9$$くらいでは、半分くらいがそうなっています。
このように、$$r^2$$で計られたLDの強さは、関連検定のp値の相関係数と関連があり、
また、強いLDにあれば、真のSNVよりp値が小さくなる確率も上がることがわかります。
```{r,warning=FALSE}
rs <- sqrt(0.1 * (1:10))
cr <- rep(0,length(rs))
n.iter <- 500
p1s <- p2s <- matrix(0,length(rs),n.iter)
for(i in 1:length(rs)){
r <- rs[i]
pH.d <- my.2SNV.hap(pA,pB,r)
p.out <- my.LD.sim(h.freq=c(pH.d),r.gen.1 =r.gen.1,n.iter=n.iter)
cr[i] <- cor(p.out$p1,p.out$p2)
p1s[i,] <- p.out$p1
p2s[i,] <- p.out$p2
}
plot(rs^2,cr,xlab="r.squared",ylab="correlation of p values")
abline(0,1,col=2)
par(mfrow=c(3,4))
frac.p1.less.than.p2 <- rep(0,length(rs))
for(i in 1:length(rs)){
ttl <- paste("r.squared=",rs[i]^2)
plot(log10(p1s[i,]),log10(p2s[i,]),main=ttl)
abline(0,1,col=2)
frac.p1.less.than.p2[i] <- length(which(p1s[i,]-p2s[i,] >= 0))/n.iter
}
par(mfcol=c(1,1))
plot(rs^2,frac.p1.less.than.p2,xlab="r.squared",ylab="fraction of p1 >= p2")
boxplot(t(log10(p2s)))
```
あるアレル・あるハプロタイプに着目すると、個人は、そのアレルを0本持つか、1本持つか、2本持つかの3通りあります。
1個のSNVについてこれを調べるには、相加モデル検定(かロジスティック回帰検定)をすればよいです。
SNVを組み合わせてハプロタイプに分解するには、SNPごとのジェノタイプを実験で調べるだけでは無理なので、ハプロタイプ推定をすることもあります。
この推定の大規模な場合がHLA領域のハプロタイプ推定になります。
シークエンシングを用いて、うまくハプロタイプごと検定することもできるかもしれません。その場合は「推定」は不要です。
アレル本数で検定する、というのは
```{r}
m <- matrix(c(20,29,24,18,6,3),byrow=FALSE,nrow=2)
m
```
というような2x3表を
```{r}
a <- matrix(c(m[1,1]*2+m[1,2],m[1,2]+2*m[1,3],m[2,1]*2+m[2,2],m[2,2]+2*m[2,3]),byrow=TRUE,nrow=2)
a
```
というような2x2表にして検定する、ということです。
カイ二乗検定してみます。
```{r}
prop.trend.test(m[1,],apply(m,2,sum))
chisq.test(a,correct=FALSE)
```
同じです。
じゃあ、トレンド検定なんていう面倒くさいことしなければよいということになりそうです。
少し数字を変えてみます。
```{r}
m <- matrix(c(23,26,24,18,6,3),byrow=FALSE,nrow=2)
m
a <- matrix(c(m[1,1]*2+m[1,2],m[1,2]+2*m[1,3],m[2,1]*2+m[2,2],m[2,2]+2*m[2,3]),byrow=TRUE,nrow=2)
a
prop.trend.test(m[1,],apply(m,2,sum))
chisq.test(a,correct=FALSE)
```
また、同じです。
```{r}
m <- matrix(c(25,26,22,18,6,3),byrow=FALSE,nrow=2)
m
a <- matrix(c(m[1,1]*2+m[1,2],m[1,2]+2*m[1,3],m[2,1]*2+m[2,2],m[2,2]+2*m[2,3]),byrow=TRUE,nrow=2)
a
prop.trend.test(m[1,],apply(m,2,sum))
chisq.test(a,correct=FALSE)
```
今度は違います。
やはり、検定手法の名前が違うということは、トレンド検定と2x2表化したカイ二乗検定は『同じなわけではない』ということですね。
では、同じになった場合と同じにならなかった場合との違いはなんなのでしょう。
同じになった場合の周辺度数を調べると、ケースとコントロールとを合わせた3ジェノタイプの人数が、
```{r}
m <- matrix(c(23,26,24,18,6,3),byrow=FALSE,nrow=2)
apply(m,2,sum)
```
となっています。
これは、
$$0.7^2 = 0.49$$
$$2\times 0.7 \times 0.3 = 0.42$$
$$0.3^2=0.09$$
という比率になっています。
これを遺伝学分野ではハーディ・ワインバーグ平衡(HWE)にある、と言いますが、HWEにあるかそうでないかが影響します。
ですから、ケース・コントロールを合わせたときのジェノタイプカウントがHWEに近ければどちらの方法でやっても結果に大差はないが、HWEからの逸脱が大きいときには、2つの手法の結果は違う、違うからにはどちらかを採用するのがよいことになる、それはトレンド検定である、ということになり、いつもトレンド検定をできるのならばそうしておくのが良いでしょう。
ハプロタイプの場合は、アレルの数がハプロタイプの種類数となりますが、考え方は同じです。
ハプロタイプの本数をケース・コントロール別に数えて検定することは、ハプロタイプ全体でHWEになっていそうかどうかの影響を受けるということになります。
ただし、実際のところ、ハプロタイプの本数推定しかできないこともあり、その場合には、ハプロタイプのペアがどうなっているのかまでわからないこともありますから、ハプロタイプ本数での検定をせざるを得ないのが現実ですし、それが許容されてもいます。
ハプロタイプ本数をケース・コントロール間で比較するときも同様です。
また、ハプロタイプ本数での検定をするときには、「ハプロタイプの種類数だけ検定を繰り返し」ますから、その分のマルチプルテスティングのことを考慮する必要があることは、もうわかると思います。
```{r}
n.ca <- 1000
n.co <- 1000
a1 <- 0.4
a2 <- 0.3
g1 <- c(a1^2,2*a1*(1-a1),(1-a1)^2)
g2 <- c(a2^2,2*a2*(1-a2),(1-a2)^2)
g12 <- matrix(g1,ncol=1) %*% matrix(g2,nrow=1)
g.ca <- rmultinom(1,n.ca,prob=c(g12))
g.co <- rmultinom(1,n.co,prob=c(g12))
matrix(c(g.ca,g.co),byrow=TRUE,nrow=2)
```
では、そのHWEについて考えてみます。
HWEとは、SNVの3ジェノタイプの割合の偏りがあるかどうかの検定です。
偏りがないというのは、2アレルの組合わせのうちホモ接合体(MM,mm)が妙に多すぎたり、少なすぎたりしないということです。
どういうときに偏るかと言えば、Mアレルが多い集団とmアレルが多い集団があって、2つの集団が別々である場合です。極端な例では、Mアレルしかない地域とmアレルしかない地域では、それぞれ、全員がMM,mmですから、この2地域からのサンプルを合わせると、ヘテロ接合体Mmが一人もいないことになります。
実際にはこれほど極端なことはないですが、地域差があって、サンプルがその混合であるときにはホモ接合体の割合が多くなります。
偏りがないときというのは、地域差がなく、よく混ざり合っているときのことで、ランダム・メイティングしているとき、と言い換えることもできます。
アレル頻度が M=0.7,m=0.3の場合には
$$
MM = 0.7 \times 0.7 = 0.49\\
Mm = 2 \times 0.7 \times 0.3 = 0.42\\
mm = 0.3 \times 0.3 = 0.09
$$
となっているのが、偏りがない状態(平衡状態)です。
HWE検定では、確かに、集団の偏り具合(集団の構造化)の有無がわかるのですが、たくさんのSNVを使ったGWASのようなスタディでは、集団の構造化については、個々のHWE検定では判断しません。あるSNVのデータを見ると、構造化していて、別のSNVのデータを見ると構造化していない、ということになったら、個人の集まりとしての構造化はどうなっているのかに対して、答えが山ほど出てしまって判断できなくなるからです。
さて、HWE検定では、もう一つ別のことを調べます。
3ジェノタイプの観察人数をHWE検定したところ、HWEから外れていると判断できたとします。
たとえば、次の例v.hwD.1 のように、ホモ接合体が多い方に外れていたり、v.hwD.2のようにホモ接合体が少ない方に外れていたりした場合です。
```{r}
v.hwE <- c(49,42,9)
v.hwD.1 <- c(55,30,15)
v.hwD.2 <- c(44,52,4)
```
検定をすると、確かにp値が小さめです。
```{r}
(p1 <- pchisq(sum((v.hwD.1-v.hwE)^2/v.hwE),1,lower.tail=FALSE))
(p2 <- pchisq(sum((v.hwD.2-v.hwE)^2/v.hwE),1,lower.tail=FALSE))
```
集団構造化によるHWEからの逸脱であれば、ホモ接合体が多めになるはずですが、ホモ接合体が少な目になることの説明はつきません。
この場合は、「本当のジェノタイプはHWEなのだけれど、実験したら、何かの理由でホモ接合体のサンプルがヘテロ接合体とコールされてしまった」からかもしれません。
実験エラーによるHWEからのずれは、ホモ接合体のコール数が多い方への偏りも起こし得るので、ホモ接合体が多く観察された場合には、集団構造化と実験エラーの両方を疑う必要があります。
このようにHWE検定では集団構造化ではなく、実験エラーを見つけるために使うこともできます。
2x3表でケース・コントロールの関連を検定したときには、p値のほかに、検定の強さを表すオッズ比のことを気にしました。
HWE検定でもHWEからの偏りの強さのことを気にします。
HWEに合致しているときが0、ヘテロ接合体が0人のときが1になるような指標があるとわかりやすいです。また、ホモ接合体が多い方に偏っているのか、少ない方に偏っているかもわかると便利です。
$$
1- \frac{\text{観測されたヘテロの割合}}{\text{アレル頻度から期待されるヘテロの割合}}
$$
のような式が好都合です。Fixation指標です。
```{r}
v.hwE <- c(49,42,9)
v.hwD.2 <- c(51,38,11)
v.hwE.3 <- v.hwE * 10
(v.hwD.3 <- v.hwD.2*10)
(f <- 1-(v.hwD.2[2]/100)/(v.hwE[2]/100))
```
HWE検定では
『サンプルは集団構造化がない母集団を代表しており、実験もうまく行っている』
というのが帰無仮説ですから、
もしHWE検定のp値が小さかったときには、
*『サンプルは集団構造化がある母集団を代表している』か
*『実験がうまく行っていない』か
のいずれかの可能性を否定しにくいことがわかります。
後述しますが、ある程度の集団構造化があってもそれに対する方策はありますので、集団構造化だけを理由に特定のSNVを排除しなくてもよいです。
HWE検定を使って排除したいのは、『実験がうまくいっていない」場合です。
また、HWE検定もマルチプルテスティングの影響を受けますから、100万マーカーを検定すれば、HWE検定のp値が100万分の1であっても、それだけで問題にはなりません。
実験の失敗による異常は、HWE検定p値が、そのほかの大多数の「成功していると思われるマーカーのp値の分布」に照らして、「外れている」ことから見当をつけるのがよいでしょう。
p値は、一様分布に従っているから、その値を0.01と聞けば、「あー、0.01的に珍しいことなんだ」とわかるわけですから、p値の本領は一様分布になっていることです。
しかしながら、実際にGWASを実施して、数十万個のp値を算出して、その分布を見てやると、一様分布になっていない。
このときに大きく分けて2つのアプローチがあります。
(1) ほとんどすべてのp値は、帰無仮説が正しいはずだから、「本当は一様分布」なはず。それなのに、一様になっていないのだったら、「一様分布」になるように修正してしまおう、という作戦。
(2) もう一つは、たくさんのp値があるけれど、結構の数の検定が、実は対立仮説が真なのだから、全部のp値が作る分布が一様分布になるわけがない。真に関連がある結果はなるべく拾いつつ、本当は帰無仮説が正しいものが偽陽性にならないように工夫する。
ケース・コントロール関連をたくさんのマーカーについて検定を繰り返すとき(GWASのように)、ほぼすべてのマーカーが疾患と無関係なので、理想的にはp値は一様分布になります。それはカイ二乗検定をしている場合で言えば、たくさんの検定統計量がカイ二乗分布に従っている、ということです。
ケース群とコントロール群とで「違いがない」のに、カイ二乗値のすべてが0にならず、大きめの値が出たりするか、と言えば、サンプリングのばらつきがあって、「たまたま、ケースとコントロールとで違いが出る」からです。
この、「たまたま出る違い」は、「ケースとコントロールをサンプリングする母集団」が「きれいに均一」な場合には、カイ二乗分布そのものになるのですが、「母集団が、似ている人、似ていない人が混ざっている」と、たとえ、「母集団からケースとコントロールとをランダムにサンプリングし」ても、「たまたま」起きる「ケースサンプルとコントロールサンプル」の違いは少し大きめになります。
こんな風に。
```{r}
n <- 100000
w <- 1.5
df <- 1
k <- rchisq(n,df)*w
k. <- rchisq(n,df)
h <- hist(c(k,k.),col=gray(0.99),border=FALSE)
hist(k,density=17,add=TRUE,breaks=h$breaks)
hist(k.,col=2,density=20,add=TRUE,breaks=h$breaks)
```
p値で見ると、小さいp値が多めになって、このようになります。
```{r}
ps <- pchisq(k,df,lower.tail=FALSE)
ps. <- pchisq(k.,df,lower.tail=FALSE)
h <- hist(c(ps,ps.),col=gray(0.99),border=FALSE)
hist(ps,density=17,add=TRUE,breaks=h$breaks)
hist(ps.,col=2,density=20,add=TRUE,breaks=h$breaks)
```
一様分布ではないp値はそのまま信じるわけにいかないので、補正が必要です。
今の状況は、「ケースとコントロールとは、もともと違いがな」く、「サンプリングも母集団からうまくされている」けれど、「母集団が均質ではなくてばらつきがある」という単純な状態なので、補正も単純にできることが知られています。
カイ二乗分布が大きめにシフトしているので、すべてのカイ二乗値をある値で割ってやって、うまく、p値が一様になるようにできることが知られています。
このうまく補正する値が、ジェノミックコントロールのλ値です。1は補正の必要がないことを意味し、1より大きければ、p値が小さめに出ているので、大きめに補正してやる必要があることを意味します。
補正前と補正後のp値をQQプロットした図を提示すれば次のようになります。補正前の赤いラインはカーブしていますが、補正後の黒い線は対角線になっていて、うまく補正で来ていることになります。
```{r}
median.k <- quantile(k,0.5)
# カイ二乗値の中央値
median.chisq <- qchisq(0.5,df)
lambda <- median.k/median.chisq
lambda
k.corr <- k/lambda
ps <- pchisq(k,df,lower.tail=FALSE)
ps.corr <- pchisq(k.corr,df,lower.tail=FALSE)
hist(ps)
hist(ps.corr)
qqplot(qunif(ppoints(n)),ps.corr,xlab="Theoretical",ylab="Observed",type="l")
points(qunif(ppoints(n)),sort(ps),col=2,type="l")
```
ジェノミック・コントロールのようなタイプの補正では、小さいp値ほど補正によるp値の増大が大きくなります。
ですから、
元のp値(p.nominal),補正後p値(p.corr),両値の比(fold)としてそれを見てみます。
```{r,warning=FALSE}
pp <- 10^(seq(from=0,to=-15,by=-0.1))
k <- qchisq(pp,1,lower.tail=FALSE)
k2 <- k * 1.1
pp2 <- pchisq(k2,1,lower.tail=FALSE)
plot(log10(pp),(pp/pp2),type="l",xlab="p.nominal",ylab="補正後値/補正前値")
print(data.frame(p.nominal = pp2,p.corr=pp,fold=pp/pp2),digits=2)
```
たしかに、GWASでゲノム全体のどこに関連のある遺伝子があるかわからないけれど、あったとしても、ごく一部の遺伝子しか、フェノタイプと関係ないだろう、と考えるときには、得られたp値が一様分布をとると仮定するのが妥当です。
そのうえで、ほとんどすべてのマーカーを関連なしとみなしたいので、「厳しく」棄却するのがよいです。
その目的にかなっているのが、ボンフェロニ法やSidak法でした。
しかしながら、２万数千個の遺伝子の発現量を組織Xと組織Yで比べたり、癌部と非癌部とで比べたりする場合には、そもそも、同じわけがない、と考えていますし、実際に、調べると、かなりの数の遺伝子の発現量が明らかに違うことがわかります。
こういうときに、厳しく切り捨てると、「かなりの『本物の関連』」も切り捨てられてしまうので、目的に適わないという問題が生じます。
False Discovery Rateと呼ばれる方法がそれです。
ボンフェロニ法はSidak法では、「棄却基準」をすべてのp値に等しく適用しました。
せっかく、たくさんの検定をしたのだから、それ全体を使って判断しよう、という考え方です。
次のQQプロットは、帰無仮説が真である検定200個と対立仮説が真である検定40個が混ざっている場合で、黒が帰無仮説を真とする結果、赤が対立仮説を真とする結果を表しています。
```{r}
k.null <- rchisq(200,1)
k.alt <- rchisq(40,1,6)
k <- c(k.null,k.alt)
p <- pchisq(k,1,lower.tail=FALSE)
col <- c(rep(1,length(k.null)),rep(2,length(k.alt)))
ord <- order(k,decreasing=TRUE)
qqplot(qunif(ppoints(length(k))),p,col=col[ord],pch=20,xlab="Theoretical",ylab="Observed")
```
赤のp値は小さい値が多いので、左下に多いです。
黒のp値は一様に分布していますが、赤よりも相対的に大きいので、右に偏って分布しています。
このようなたくさんのp値が得られたときに、赤・黒の情報はないわけですが、赤は、「拾っ」て、黒は「捨て」たいので、そのような方法を作ろう、という考え方です。
ボンフェロニ法やSidak法では、以下のように水平線を引いて、その線より下に来た点を「拾う」方法です。
```{r}
qqplot(qunif(ppoints(length(k))),p,col=col[ord],pch=20,xlab="Theoretical",ylab="Observed")
abline(h=0.01,col=4)
```
FDRは、水平線ではない基準を入れる方法で、たとえば、原点を通る次のような直線を引きます。
```{r}
qqplot(qunif(ppoints(length(k))),p,col=col[ord],pch=20,xlab="Theoretical",ylab="Observed")
abline(0,0.1,col=4)
```
p値が小さい部分を拡大してみます。
```{r}
qqplot(qunif(ppoints(length(k))),p,col=col[ord],pch=20,xlab="Theoretical",ylab="Observed",xlim=c(0,0.2))
abline(0,0.1,col=4)
```
なぜ、水平線にしないかを考えてみます。
たくさんのp値があって、それが一様分布に従っているときに、一番小さいp値は、小さい値を中心に少しばらつきますし、一番大きいp値は、大きい値を中心にばらつきます。
ですから、たくさん得られたp値を小さい順に並べたとき、一番小さいp値の場合にそれが、「帰無仮説を信じるには小さすぎる」という基準の値は、二番目に小さいp値のための基準値よりも小さくしてもよさそうです。
もちろん、三番目に小さいp値の場合はさらに緩めにすることができるでしょう。
このように何番目に小さいかによって基準値を変えてやります。
そうは言っても2番目に小さいp値が2番目用の基準値をクリアしなかったけれど、3番目に小さいp値が3番目用の基準値をクリアしたからと言って、2番目に小さいp値の結果は、関連なしとし、3番目は関連ありとするのは、いかにもおかしいです。
したがって、何番目に小さいかによって基準を定め、そういう基準をクリアしたp値を見つけたら、それ以下のp値はすべて「関連あり」側にする、という約束も付け加えます。
この直線の引き方としては、傾きが0.05(たとえばの水準)を引くことです。
こんな風に考えます。
たくさんの検定のうち、最も大きいp値を出したものについて、それが帰無仮説を棄却する場合を考えるとします。
一番大きいp値を棄却することとしたら、それより小さいすべてのp値の検定の帰無仮説を棄却することとなります。
すると、今、「帰無仮説が真」かもしれない仮説は、問題にしている検定の対象だけです。たった１つです。
したがって、その検定を棄却するかどうかは、単独の検定について考えればよいから、通常の0.05や0.01を使えばよいです。
それ以外のp値(二番目、三番目…に大きいp値)については、別途基準を決める必要がありますが、思い切って、原点からの直線を引いてしまえ、というやり方をここでは示しました。
他にも複数のやり方が提唱されていますが、大まかな発想を理解するのには、この方法で十分でしょう。
