540

# 【第3回】ゼロから始めるゲノム解析（R編）

2020/10/09 【第3回】ゼロから始めるゲノム解析（R編） 資料
https://bioalgorithms.connpass.com/event/190366/

October 09, 2020

## Transcript

2. ### 統計の背後の考え • 実験の技術的偏り ◦ 技術者の腕 ◦ シーケンサー • 実験環境 ◦

温度 ◦ 湿度 • サンプルの自然変動 ◦ サンプル劣化 FIGURE 3.1: PAX6遺伝子から 20回遺伝子発現量を取得した結果 Why？ 生物学や他の多くの分野では、データは実験によって収集される。 様々な要因により、毎回同じ測定値を取得することはできない。 観測値がばらつく原因 2
3. ### 平均と中央値 平均 median()関数 numeric型ベクトルの中央値を計算 中央値 （n: 奇数） （n: 偶数） mean()関数

numeric型ベクトルの平均を計算 > runif(10) > mean(x) ##  0.3738963 > runif(10) > median(x) ##  0.3277896 生物学や他の多くの分野では、データは実験によって収集される。 中央値は外れ値の影響を受けづらいという特徴がある。 runif(n) は一様分布からn 個乱数生成する関数 runif(n) は一様分布からn 個乱数生成する関数 3
4. ### 分散と標準偏差 分散 sd()関数 numeric型ベクトルの標準偏差を計算 標準偏差 var()関数 numeric型ベクトルの分散を計算 > x=rnorm(20,mean=6,sd=0.7) >

var(x) ##  0.2531495 > x=rnorm(20,mean=6,sd=0.7) > sd(x) ##  0.5031397 データを要約する他の指標として、値の変動性を数値化することが考えられる。 「分散」と「標準偏差」はこの指標として用いられる。 rnorm(n) は正規分布からn 個乱数生成する関数 rnorm(n) は正規分布からn 個乱数生成する関数 4
5. ### 四分位範囲と分位数 値の大部分がどこに分布しているのかを調べることも役に立つ。 四分位範囲は分散と異なり、外れ値の影響を受けにくいという特徴を持つ。 分位数 quantile()関数 numeric型ベクトルの 分位数を計算 四分位範囲 IQR()関数 numeric型ベクトルの

四分位範囲を計算 > x=rnorm(20,mean=6,sd=0.7) > quantile(x) ## 0% 25% ## 5.437119 5.742895 ## 50% 75% 100% ## 5.860302 6.243991 6.558112 > x=rnorm(20,mean=6,sd=0.7) > IQR(x) ##  0.5010954 > boxplot(x,horizontal = T) 四分位範囲 5

7. ### Rで正規分布 * norm関数のファミリー（rnorm、dnorm、qnorm、およびpnorm）を用いることで正規分布から得 られる様々な値を取得することができます。 # 平均0、標準偏差2の正規分布から、 X=-2 における確率密度関数を取得 > dnorm(-2,

mean=0, sd=2) ##  0.1209854 # 平均0、標準偏差2の正規分布から、 P(X =< -2)を取得 > pnorm(-2, mean=0, sd=2) ##  0.1586553 # 平均0、標準偏差2の正規分布から 5つの乱数を取得 > rnorm(5, mean=0, sd=2) ##  -1.8109030 -1.9220710 -0.5146717 0.8216728 -0.7900804 # 平均0、標準偏差2の正規分布から P(X =< t)=0.15 を満たすtを取得 > qnorm(0.15, mean=0, sd=2) ##  -2.072867 7
8. ### その他Rで使える代表的な分布関数 # 成功確率が 0.3 である試行を 10 回行ったときに、 4 回成功する確率を計算。 >

dbinom(4, n=10, p=0.3) ##  0.2001209 # 期待値λ=1のポアソン分布でを取得 > dpois(1, lambda=1) ##  0.3678794 # 平均0、標準偏差2の負の二項分布から 5つの乱数を取得 > dnbinom(1, size=1, prob=0.1) ##  0.09 # 自由度(10,20)のF分布から、X=1 における確率密度関数の値を取得 > df(1, df1=10, df2=20) ##  0.7143568 # 自由度1のχ二乗分布から、 X=1 における確率密度関数の値を取得 > dchisq(1, df=1) ##  0.2419707 負の二項分布は、シーケンス読み取りカウ ントなどのモデル化に登場！ 二項分布は、二値（0, 1）がある確率で生じ る場合に登場！（ex: メチル化） F分布やχ二乗分布は線形モデルや一般 化線形モデルを使用する際に登場！ 8

12. ### # 必要なパッケージの読み込み library(mosaic) # シード値の設定 (毎回同じ結果にするための設定 ) set.seed(21) sample1= rnorm(50,20,5)

# サンプルを取得 # 重複を許容し、ブートストラップサンプリングを実行 boot.means= do(1000) * mean(resample(sample1)) # 上記で得られたサンプルからパーセンタイルを取得 q=quantile(boot.means[,1],p=c(0.025,0.975)) # ヒストグラムをプロット hist(boot.means[,1], col="cornflowerblue", border="white", xlab="sample means") # パーセンタイル値に赤線を引く abline(v=c(q, q ),col="red") text(x=q,y=200,round(q,3),adj=c(1,0)) text(x=q,y=200,round(q,3),adj=c(0,0)) 信頼区間の導出（リサンプリング使用） 1000個のブートストラップサンプルを使用した サンプル平均の精度推定 12
13. ### 中心極限定理（CLT） 平均 、分散　 の母集団 から無作為抽出した 個の標本の平均　 は 平均 、分散　　　　の正規分布に近づく。 　　　は標準正規分布に従うため、

Zスコアが算出可能 Zスコアに変換した値が特定の範囲内に存在する確率か ら、信頼区間を導出することが可能！ 標本平均は、母集団分布に関係なく正規分布に従う 平均が0、標準偏差 (SD) が1に なるように変換した得点 13

15. ### Rで標準正規分布から信頼区間の導出 # 中心極限定理から数理的に 95%信頼区間を推定 sample1= rnorm(50,20,5) # サンプルを取得 alpha=0.05 sd=5

n=50 # 標準正規分布から信頼区間を推定 mean(sample1) + qnorm(c(alpha/2, 1-alpha/2)) * sd/sqrt(n) ##  19.23638 22.00819 ブートストラップで得られる信頼区 間 [19.21, 21.989] 中心極限定理から数理的に 得られる信頼区間 [19.23638, 22.00819] qnorm() はZスコアを算出する関数 15
16. ### 母集団の標準偏差分からない問題 一般に母集団の統計量（母数）は分からない … ↓ でも標本の標準偏差（　）を使用して標準偏差を推定したい！ ↓ 信頼区間の計算に標準正規分布の代わりに t分布（→） を使用すれば、　　　　　　　　から信頼区間が導出可能！ 標本から母集団の信頼区間を求める際は

t分布を使用する　 標本 母集団 推定！ 16
17. ### Rでt分布から信頼区間の導出 ブートストラップで得られる信頼区 間 [19.21, 21.989] # 中心極限定理から数理的に 95%信頼区間を推定 sample1= rnorm(50,20,5)

# サンプルを取得 alpha=0.05 sd=5 n=50 # t分布から信頼区間を推定 mean(sample1) + qt(c(alpha/2, 1-alpha/2), df=n-1) * sd/sqrt(n) ##  19.20130 22.04327 qt() はt分布の分位点を算出する関数 t分布から数理的に得られる信頼区 間 [19.20130, 22.04327] 17

B群の 平均値 20
21. ### どのようにサンプルセット間比較を行うか 最も簡単な方法：2つのサンプルの平均値の差によって判断する。 ⇒ 平均値に差があればサンプルセット同士には差があると判断し、 　 差がなければサンプルセット同士には差がないと判断する。 知っておかなければならない点：測定値にはばらつき(誤差)がある。 (測定器に依存した誤差、測定者に依存した誤差、環境条件に依存した誤差) そのため、サンプル比較する際にはその誤差を考慮する必要がある。 ⇒平均値の差には誤差が含まれることを理解する。

(ある程度の差が生まれてしまうことは仕方がない。) 21

平均値 22
23. ### 仮説検定の手順 ①帰無仮説と対立仮説を設定する。 「帰無仮説(H 0 )」：検証したい仮説 ⇒今回の場合、帰無仮説は、「サンプル間に差がない」こと。 「対立仮説(H 1 )」：帰無仮説に対する仮説 (本来証明したい仮説)

⇒今回の場合、対立仮説は、「サンプル間に差がある」こと。 ②帰無仮説の真偽を検証するための統計量 (今回の場合、差があると言える閾値 )を決定する。 ③統計値を計算し、それを参照値 (今回の場合、平均値の差)と比較して、意味のある差(有意差)がある のかをP値を用いて判断する。 ④結果、帰無仮説H 0 を棄却するか棄却しないかを判断する。 23

25. ### ランダム化を用いた平均値の差の検証方法 概念 サンプルセット間に差がないと思われる場合は、 サンプルセットのラベル(例: “健康” と “病気”)には意味がないはず。 2つのサンプルセットをごちゃ混ぜにして、 再度ランダムに2つのサンプルセットに分けた時の平均値の差と 実際のサンプルセットの平均値の差を比較して、

サンプルセットのラベルには意味がなかったのかを検証する。 ⇒ 実際にサンプルセット間には違いがなかったのか（帰無仮説）を検証する。 25
26. ### ランダム化を用いた平均値の差の検証方法 手順 1) サンプルセットにラベルをランダムに割り当てて平均値 の差を計算することを繰り返す（可能であれば全ての 組み合わせを試す）。 ⇒平均値の差の分布（帰無分布）が計算できる。 2) 帰無仮説が成り立つ場合に、実際の平均値の差の値 を得る確率（頻度）を計算する。

3) 得た確率が設定した有意水準以下であれば、有意な 差があると判断する。（帰無仮説を棄却する。） ランダムセットの 平均値の差の分布(帰無分布) 有意 水準 実際の 平均値の差 26
27. ### # シード値を設定 (毎回同じ結果にするための設定 ) set.seed(100) # 2つの異なる分布から 2つのサンプルセットを取得 gene1=rnorm(30,mean=4,sd=2) gene2=rnorm(30,mean=2,sd=2)

# 平均値の差を計算 org.diff=mean(gene1)- mean(gene2) gene.df=data.frame(exp=c(gene1,gene2), group=c(rep("test",30),rep("control",30))) # 帰無分布を取得 exp.null <- do(1000) * diff(mosaic::mean(exp ~ shuffle(group), data=gene.df)) # ヒストグラムを描画 hist(exp.null[, 1],xlab="null distribution | no difference in samples" , main=expression(paste(H," :no difference in means" )), xlim=c(-2,2),col="cornflowerblue" ,border="white") abline(v=quantile(exp.null[, 1],0.95),col="red" ) abline(v=org.diff,col="blue" ) text(x=quantile(exp.null[, 1],0.95), y=200, "0.05",adj=c(1,0), col="red") text(x=org.diff,y=200,"org. diff." ,adj=c(1,0),col="blue") Rでの実装 27

29. ### 3.2.2 Using t-test for difference of the means between two

samples t検定による2つのサンプルセット間の平均値の比較 29

31. ### スチューデントのt検定 (母分散が等しいと仮定する場合) 母分散が等しいと仮定すると、次の式が成り立つ。 上式は、平均値の差が自由度n 1 +n 2 -2のt分布に従うことを意味する。 スチューデントのt検定はRで実装されており、自由度n 1

+n 2 -2のt分布における、計算し た統計量tの値のP値が取得できる。 , 31

32
33. ### Rで実装 (ウェルチのt検定) # Welch’s t-test stats::t.test(gene1, gene2) # ## ##

Welch Two Sample t-test ## ## data: gene1 and gene2 ## t = 3.7653, df = 47.552, p-value = 0.0004575 ## alternative hypothesis: true difference in means is not equal to 0 ## 95 percent confidence interval: ## 0.872397 2.872761 ## sample estimates: ## mean of x mean of y ## 4.057728 2.185149 33
34. ### Rで実装 (母分散が等しいと仮定したt検定) # t-test with equal variance assumption stats::t.test(gene1, gene2,

var.equal=TRUE ) # 出力 ## ## Two Sample t-test ## ## data: gene1 and gene2 ## t = 3.7653, df = 58, p-value = 0.0003905 ## alternative hypothesis: true difference in means is not equal to 0 ## 95 percent confidence interval: ## 0.8770753 2.8680832 ## sample estimates: ## mean of x mean of y ## 4.057728 2.185149 34

35

37. ### 検定の評価方法 TN/(FP+TN): 特異度 TP/(TP+FN): 感度(検定の検出力) 帰無仮説を 採択した 帰無仮説を 棄却した 帰無仮説が

真であった 帰無仮説が 偽であった (帰無仮説が真である数 ) (対立仮説が真である数 ) (第2種の過誤) (第1種の過誤) 高くすることが目標 37
38. ### 多重検定の修正 検定の数が増えると… 第1種の過誤(False Positives)が増えると予想される。 例えば、5％の有意水準で検定を実行すると、 帰無仮説が真である場合、帰無仮説を誤って棄却する可能性が5%ある。 帰無仮説が全て当てはまる1000回の検定を行うと、 誤った棄却の平均の数は50になる。 また、少なくとも誤った棄却が1回はほぼ100%の確率で発生する。 38

10,000遺伝子についてマイクロ アレイの発現量を比較する場合 など
39. ### 多重検定の修正法① 多重検定における第1種の過誤を防ぐための統計手法は複数ある。 ①ボンフェローニ補正 全検定中でひとつでも偽陽性が得られる確率がある目標値となるように、 　複数の検定から得られたP値をより高い値に補正する。 　具体的には、個々のP値が十分に低い場合に、 　個々のP値(p i )に検定数(m)を掛ける(m×p i

)。 問題点： 陽性判定の条件が厳しくなりすぎるため、帰無仮説を棄却しづらい 39
40. ### 多重検定の修正法② ②Benjamini-Hochberg法 ある程度第1種の過誤を許容して、第2種の過誤を起こす可能性を小さくする方 法の一つとして、False Discovery Rateを調整するという方法が考案された。 P値をランク付けし、検定数(m) × 個々のP値(p i

) をランクiで割る (m×p i / i) 　上記の手順により、「FDR調整済みP値（=q-value）」を取得し、 　特定のテストでの誤検出の割合を見積もることができる。 FDR調整済みP値が0.05であることは、 棄却された帰無仮説の5％が誤りであることを意味する。 問題点： p-valueの分布を一様分布と仮定している。 False Discovery Rate (FDR) : 真の帰無仮説のうち誤って棄却した割合 FP / (TP + FP) 40

="fdr") plot(hedenfalk\$p,qvalues,pch=19,ylim=c(0,1), xlab="raw P-values" ,ylab="adjusted P-values" ) points(hedenfalk\$p,bonf.pval,pch=19,col="red") points(hedenfalk\$p,fdr.adj.pval, pch=19,col="blue") legend("bottomright",legend=c("q-value","FDR (BH)","Bonferroni"), fill=c("black","blue","red")) Rでの実装 出力 41

43. ### 変数間の関係を測定/モデル化したい in Genomics 薬を X 投与！ 肝臓である遺伝子 発現量が Y 上昇！

X 箇所がメチル化！ メチル化された遺伝 子発現量が Y 減少！ X と Y の関係をモデル化したい！ 43

45

48. ### βの値によるコスト関数の大きさ コスト関数/損失関数を用いたアプローチ 「コスト関数」もしくは「損失関数」と呼ばれる関数を最適化する。 一般的にこの関数には、予測値 と観測値 の二乗誤差が使用される。 最適化手順（最急降下法） 1. ランダムな開始点（ ,

）を選択する。 2. コスト関数の偏導関数から、コスト関数をどの方向に最適化するかを決定。 3. コスト関数を最小化する方向に向かって最適化を進める。ステップの大きさはパラ メータとして与える。 4. 収束するまで手順2、3を繰り返す。 48
49. ### 平均 、分散 に従う正規分布を仮定した場合、応答変数を観測する 確率は下記となる。 最尤関数を用いたアプローチ① 最尤推定法では、与えられたデータセット内の全ての応答変数を観測する確率を最大 化するような を求める。 尤度関数 は、全てのデータポイントについて観測される確率の乗算となる。

対数変換 MAX y=ln(x) は連続 49

50
51. ### Rで線形回帰 # 乱数シードを設定 set.seed(32) # 1から100までの50個のX値を取得 x = runif(50,1,100) #

b0とb1、分散（sigma）を設定 b0 = 10 b1 = 2 sigma = 20 # 正規分布から誤差項をシミュレート eps = rnorm(50,0,sigma) ＃線形方程式と誤差項の追加からy値を取得 y = b0 + b1*x+ eps mod1=lm(y~x) # データをプロット plot(x,y,pch=20, ylab="Gene Expression", xlab="Histone modification score") # 回帰直線をプロット abline(mod1,col="blue") 線形回帰によってモデル化された遺伝 子発現とヒストン修飾スコア lmはlinear modelの略 51

52
53. ### 回帰係数の推定標準誤差 標本を使用して係数を推定しているため、 ランダムサンプリングごとに異なる推定値が得られる。 誤差項の分散が一定で平均がゼロであると仮定すると、標本 から母集団の特性を推測するための特性として、 回帰係数の不確実性を下記の通りモデル化することができる （証明は略）。 回帰係数は、ランダムサンプルごとに異なる 回帰係数の変動をヒストグラムで表示 は

RSE、及びXの 分散に依存 53
54. ### Rで線形回帰（結果の確認） # 線形回帰結果を取得 > mod1=lm(y~x) > summary(mod1) ## Call: ##

lm(formula = y ~ x) ## Residuals: ## Min 1Q Median 3Q Max ## -38.50 -12.62 1.66 14.26 37.20 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 0.40951 6.43208 0.064 0.95 ## x 2.08742 0.09775 21.355 <2e-16 *** ## --- ## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 ## Residual standard error: 16.96 on 48 degrees of freedom ## Multiple R-squared: 0.9048, Adjusted R-squared: 0.9028 ## F-statistic: 456 on 1 and 48 DF, p-value: < 2.2e-16 # βの信頼区間を取得 > confint(mod1) ## 2.5 % 97.5 % ## (Intercept) -12.523065 13.342076 ## x 1.890882 2.283952 # βの推定値を抽出 > coef(mod1) ## (Intercept) x ## 0.4095054 2.0874167 推定値、標準誤差、t値、p値 を表示 54

58

60. ### カテゴリ変数を用いた回帰 # 線形回帰結果を取得 set.seed(100) gene1=rnorm(30,mean=4,sd=2) gene2=rnorm(30,mean=2,sd=2) gene.df=data.frame(exp=c(gene1,gene2), group=c(rep(1,30),rep(0,30))) mod2=lm(exp~group,data=gene.df) summary(mod2)

## ## Call: ## lm(formula = exp ~ group, data = gene.df) ## ## Residuals: ## Min 1Q Median 3Q Max ## -4.7290 -1.0664 0.0122 1.3840 4.5629 ## ## Coefficients: ## Estimate Std. Error t value Pr(>|t|) ## (Intercept) 2.1851 0.3517 6.214 6.04e-08 *** ## group 1.8726 0.4973 3.765 0.000391 *** ## --- ## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ## ## Residual standard error: 1.926 on 58 degrees of freedom ## Multiple R-squared: 0.1964, Adjusted R-squared: 0.1826 ## F-statistic: 14.18 on 1 and 58 DF, p-value: 0.0003905 require(mosaic) plotModel(mod2) 60

62. ### 線形回帰における落とし穴 非線形性 真の関係性が線形からかけ離れている場合、予測精度が低下する 場合によって , , 等の変換を検討する必要がある。 説明変数の相関 （多重共線性） 説明変数間の相関が高い場合、係数を適切に推定できない。

係数の推定標準偏差が大きくなる傾向が見られる。 誤差項の相関 誤差項が互いに無相関であることを前提としている。 相関がある場合、係数の信頼区間が狭くなる傾向がある。 誤差項の非一定分散 予測変数の値に関係なく、異なる応答変数の誤差に同じ分散がある。 エラーが一定でない場合（例： X値が増加するとエラーが増加します）、 標準誤差の推定値が信頼できなくなる。 外れ値と高レバレッジ 点 外れ値：Y に関する極値、高レバレッジ点：異常な X値を指す。 これらは、近似直線と標準誤差に影響を与える可能性がある。 場合によってはデータから削除して適合性を高めるべき。 62