【第3回】ゼロから始めるゲノム解析 補足説明と次回予告

f:id:kimoppy126:20201018201441p:plain

先日から隔週で勉強会を開催しています。

bioalgorithms.connpass.com

前回使用した資料は下記で公開しているので是非参考にしてみてください。

第二章のR基礎までは基礎的な内容がほとんどであったため、問題なく理解できた方が多いかと思いますが、
前回の第3回から少し内容が高度(統計の知識が必要)となってきたため、以後の講義はこちらのブログにて補足説明を用意しようと思います。

p.7 Rで正規分布の種々の値

正規分布関連の各関数について、具体的にプロットを示しながら説明します。 Rでは正規分布以外(ポアソン分布、二項分布、負の二項分布、F分布、χ二乗分布、etc...)に関しても同様の関数が使用できます。

dnorm

dnormのdはdensityのd。正規分布の確率密度の値を返します。

# 平均0、標準偏差2の正規分布から、X=-2 における確率密度関数を取得
dnorm(-2, mean=0, sd=2)
## [1] 0.1209854
x <- seq(-10, 10, length=1000)
y <- dnorm(x, mean=0, sd=2)
plot(x,y,type="l")

pnorm

pnormのdはprobabilityのp。確率密度関数の値を受けて、確率(分布関数の値)を返します。

# 平均0、標準偏差2の正規分布から、P(X =< -2)を取得
pnorm(-2, mean=0, sd=2)
## [1] 0.1586553
x <- seq(-10, 10, length=1000)
y <- dnorm(x, mean=0, sd=2)
plot(x,y,type="l")
x <- seq(-10, -2, length=400)
y <- dnorm(x, mean=0, sd=2)
polygon(c(-10, x, -2), c(0, y, 0), col="red")
text(-3, 0.035,"0.1586553")

qnorm

qnormのqはquantileのq。確率(分布関数の値)から正規分布の分位数を返します。

# 平均0、標準偏差2の正規分布からP(X =< t)=0.15を満たすtを取得
qnorm(0.15, mean=0, sd=2)
## [1] -2.072867
x <- seq(-10, 10, length=1000)
y <- dnorm(x, mean=0, sd=2)
plot(x,y,type="l")
x <- seq(-10, -2, length=400)
y <- dnorm(x, mean=0, sd=2)
polygon(c(-10, x, -2), c(0, y, 0), col="red")
text(-3, 0.035,"0.15")
arrows(0, 0.05, -1.9, 0,length=.15)
text(1.6, 0.06, "-2.072867")

rnorm

rnormのrはrandomのr。正規分布の確率密度から取得した乱数を返します。

# 平均0、標準偏差2の正規分布から5つの乱数を取得
rnorm(5, mean=0, sd=2)
## [1] -0.05621737 -1.64373955  1.84492780  2.03735084  0.58431167
x <- seq(-10, 10, length=1000)
y <- dnorm(x, mean=0, sd=2)
plot(x,y,type="l")
arrows(0, 0.1, 3, 0.15, length=.15)
text(4, 0.16, "× 5")

p.11 そもそも信頼区間とは

95%信頼区間、というと、直感的には真の値が95%の確率で含まれる範囲のような印象を受けます。
一方で、実際のところ母数は定hjk数であって確率変数ではないため、そんな範囲は求められません。正確には下記のような定義になります。

x%信頼区間の定義:

母数がx% で含まれる範囲

母集団から標本を取ってきて、その統計量からx%信頼区間を求める、という施行を行った際に、その区間の中に母数が含まれるような確率がx% となるような範囲

p.27 mosaicパッケージについて

帰無分布を求める際に下記のような構文を用いています。
ここでは、ランダムに割り付けたグループ間の差分を算出することを1000回繰り返しています。

# 帰無分布を取得
exp.null <- do(1000) * diff(mosaic::mean(exp ~ shuffle(group), data=gene.df))

mosaic::do

n 回処理を繰り返すことを直感的な構文で書くことができます。

mosaic::mean

mosaicパッケージ利用下で平均を求めるための関数です。
第一引数にformula オブジェクト、dataにデータフレームを指定することで、指定したformula で算出した値をデータフレームとして返してくれます。

※ 同じ関数名でbase::meanがあるため、mosaic::meanとすることで明示的に指定しています。

# データフレームの確認
> head(gene.df)
       exp group
1 2.995615  test
2 4.263062  test
3 3.842166  test
4 5.773570  test
5 4.233943  test
6 4.637260  test

# groupをランダムに割り当てて各グループごとの平均を算出
> mosaic::mean(exp ~ shuffle(group), data=gene.df)
 control     test 
3.451785 2.791092 

p.31, p.32 t分布の名前の由来

Gossetは自分の分布を zと言及していた。1922年から1925年の間にFisher がそれを修正 して、統計的検定の可能性を発展させるのを助けた。Gosset はその修正を t と呼んだ。 t検定 名前の由来

らしいです。

p.31, p.32 t検定における仮定について

サクッとt検定に関してお話ししましたが、正しいt検定を行うには以下の2つの仮定に気を付ける必要があります。

  1. データが正規性を持つ
  2. 母分散が等しい(異なる)

データが正規分布に従うことをどう確かめるか

シャピロ–ウィルク検定を行うことでデータの正規性を検定することができます。

帰無仮説 (H0) は標本分布が正規分布に従うことです。

帰無仮説が棄却された場合、正規分布に従わないとみなされるためt検定を使用することはできません。
R では,デフォルトでインストールされているパッケージの関数 shapiro.test() にて本検定を行うことができます。

> set.seed(100)
> gene1=rnorm(30,mean=4,sd=2)
> gene2=rnorm(30,mean=2,sd=2)

> shapiro.test(x=gene1)

    Shapiro-Wilk normality test

data:  gene1
W = 0.93417, p-value = 0.06341

> shapiro.test(x=gene2)

    Shapiro-Wilk normality test

data:  gene2
W = 0.95753, p-value = 0.2675

p-value $\gt$ 0.05となっているため、各データの正規性は棄却されません(正規分布に従わないとはいえない)。

母分散が違うことをどう確かめるか。

F検定を行うことで2つの分布について、母分散の等分散性を確かめられます。

帰無仮説 (H0) は2群間の母分散は等しいことです。

帰無仮説が棄却された場合、2つの分布の母分散は等しくないとみなされるため、ウェルチのt検定が使用されます。
一方で、帰無仮説が棄却されなかった場合、2つの分布の母分散は等しいとみなされるため、スチューデントのt検定が使用されます。

p.35 t検定の正規性について

z値やt値は平均値を標準化した値なので、結局取り扱っている分布は平均値の分布になります。
中心極限定理から、母集団がどのような分布であっても平均値の分布は正規分布に従うため、サンプルサイズが大きければ多少正規分布から逸脱していてもt検定の妥当性が保証されることになります。

p.38 多重検定を行う分析例

多重検定を行う分析の代表例として、マイクロアレイを用いた発現量変動解析が挙げられます。
異なる2つの条件でDNAマイクロアレイを使ったゲノムワイドな遺伝子発現量の測定を行なって発現量が有意に異なる遺伝子を探し出そうとする実験の場合、DNAマイクロアレイで測定した遺伝子の数だけ統計検定を行なうことになります。

p-value < 0.05 を有意とみなす検定であっても、遺伝子1,000個について発現量の差を比較する場合では1,000回検定することになり、結果平均50回有意性が得られてしまいます。
また、1,000回の検定中、偽陽性が生じる確率は $1-(1-0.05^{10000}) \simeq1.0$となりほぼ100 %の確率で偽陽性が含まれることになります。

p.39 FWER制御法(ボンフェローニ補正、ホルム法など)

FWER制御法ではp-valueに着目した上で、「N回の検定をした際の、1回も偽陽性が生まれない確率」を制御しています。

ボンフェローニ補正

ボンフェローニ補正では「補正後のp-value = 検定数 * 元のp-value」として、p-valueを補正します。
一方で、現実に測定から得られるデータはボンフェローニ補正を行った有意水準を超えるほど極端な値をとりません。
⇒ マイクロアレイなど検定数が多い超多重検定では、有意判定を行える遺伝子がなくなってしまいます。

参考:大規模データの解析における問題点

# p-value 0.001から0.01まで0.001刻みのp-valueのが得られたとする
> p = (1:10) * 0.001
> q.values.bonferroni <- p.adjust(p, "bonferroni")
> q.values.bonferroni
 [1] 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.10
> q.values.bonferroni < 0.05
 [1]  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE

ホルム法

ホルム法では、棄却された仮説は帰無仮説に従わないことを利用し、 p-valueの大きさに従って異なるα水準を用います。

  1. N個の帰無仮説を、p-valueの小さい順に並べます。
  2. 最もp-valueが小さい第1順位の帰無仮説の有意水準をα/Nにします。 「p-value < α/N」であれば、第1順位の帰無仮説を棄却し、対立仮説を採択します。 「p-value > α/N」であれば、第1順位以下のすべての帰無仮説の判定を保留します。
  3. 第1順位の帰無仮説が棄却された場合、第2順位の帰無仮説の有意水準をα/(N-1)にします。 「p-value < α/(N-1)」であれば、第2順位の帰無仮説を棄却し、対立仮説を採択します。 「p-value > α/(N-1)」であれば、第2順位以下のすべての帰無仮説の判定を保留します。
  4. 上記を繰り返します。第(k-1)順位の帰無仮説が棄却されたならば、第k順位の有意水準をα/(N+1-k)にします。p値 < α/(N+1-k)となる最大のkを見つけ出し、第1~k順位の帰無仮説を棄却し、第(k+1)順位以降の帰無仮説の判定を保留します。
# p-value 0.001から0.01まで0.001刻みのp-valueのが得られたとする
> p = (1:10) * 0.001
> q.values.holm <- p.adjust(p, "holm")
> q.values.holm
 [1] 0.010 0.018 0.024 0.028 0.030 0.030 0.030 0.030 0.030 0.030
> q.values.holm < 0.05
 [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

p.40 FDR制御法(Benjamini-Hochberg法、Q-valueなど)

FWER制御法では「1回も偽陽性が生まれない確率」を制御していたのに対し、 FDR制御法では、ある程度の偽陽性を許容するべく、FDR = FP / (FP + TP) に着目します。
これは陽性判定されたもののうちの偽陽性の割合であり、「複数の仮説が棄却されることが想定できる場合」有用な制御法となります。

Benjamini-Hochberg法 (BH法)

Benjamini-Hochberg法 では、仮説のP値は0から1の間で一様に分布しているとした上で、補正済みp-valueを算出します。

# p-value 0.001から0.01まで0.001刻みのp-valueのが得られたとする
> p = (1:10) * 0.001
> q.values.bh <- p.adjust(p, method = "BH")
> q.values.bh
 [1] 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01 0.01
> q.values.bh < 0.05
 [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

Q-value

補正したp-valueを補正前p-valueと区別してq-valueと呼びますが、それとは別にQ-valueという補正方法があります。

先のBH法で補正したq-valueでは、p-valueの分布を一様分布と仮定していました。一方で、実際のp-valueの分布は、一様分布である帰無仮説由来の分布(T0)と0に寄った対立仮説由来の分布(T1)の混合分布だとみなされます。
Q-valueでは、帰無仮説と対立仮説の密度比が $\pi_0:1 - \pi_0$ とみなした上で、BH法のq-valueに帰無仮説が真であるものの割合である$π_0$を掛けた値を用います。

次回予告

次回第4回は、教師なし機械学習による探索的データ分析がテーマとなります。

bioalgorithms.connpass.com

資料の範囲としてはこちらになります。第3回の統計学を踏まえ、いよいよ機械学習を用いた分析に関して扱っていくため、機械学習に馴染みのない方はぜひ復習・予習して臨んで頂けると嬉しいです。

一方で、先日の反省として資料の内容を網羅しようとすると講義自体が分かりにくくなってしまうと感じたため、今後はより基礎的な説明の丁寧さを心がけて説明するように方向修正していきたいと思います。

ここまで参加して頂いた方は引き続き、初めての方はぜひ今回からでもご参加よろしくお願いします。