Lock in $30 Savings on PRO—Offer Ends Soon! ⏳

【第3回】ゼロから始めるゲノム解析(R編)

nkimoto
October 09, 2020

 【第3回】ゼロから始めるゲノム解析(R編)

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

nkimoto

October 09, 2020
Tweet

More Decks by nkimoto

Other Decks in Science

Transcript

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

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

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

    var(x) ## [1] 0.2531495 > x=rnorm(20,mean=6,sd=0.7) > sd(x) ## [1] 0.5031397 データを要約する他の指標として、値の変動性を数値化することが考えられる。 「分散」と「標準偏差」はこの指標として用いられる。 rnorm(n) は正規分布からn 個乱数生成する関数 rnorm(n) は正規分布からn 個乱数生成する関数 4
  4. 四分位範囲と分位数 値の大部分がどこに分布しているのかを調べることも役に立つ。 四分位範囲は分散と異なり、外れ値の影響を受けにくいという特徴を持つ。 分位数 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) ## [1] 0.5010954 > boxplot(x,horizontal = T) 四分位範囲 5
  5. Rで正規分布 * norm関数のファミリー(rnorm、dnorm、qnorm、およびpnorm)を用いることで正規分布から得 られる様々な値を取得することができます。 # 平均0、標準偏差2の正規分布から、 X=-2 における確率密度関数を取得 > dnorm(-2,

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

    dbinom(4, n=10, p=0.3) ## [1] 0.2001209 # 期待値λ=1のポアソン分布でを取得 > dpois(1, lambda=1) ## [1] 0.3678794 # 平均0、標準偏差2の負の二項分布から 5つの乱数を取得 > dnbinom(1, size=1, prob=0.1) ## [1] 0.09 # 自由度(10,20)のF分布から、X=1 における確率密度関数の値を取得 > df(1, df1=10, df2=20) ## [1] 0.7143568 # 自由度1のχ二乗分布から、 X=1 における確率密度関数の値を取得 > dchisq(1, df=1) ## [1] 0.2419707 負の二項分布は、シーケンス読み取りカウ ントなどのモデル化に登場! 二項分布は、二値(0, 1)がある確率で生じ る場合に登場!(ex: メチル化) F分布やχ二乗分布は線形モデルや一般 化線形モデルを使用する際に登場! 8
  7. # 必要なパッケージの読み込み 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[1], q[2] ),col="red") text(x=q[1],y=200,round(q[1],3),adj=c(1,0)) text(x=q[2],y=200,round(q[2],3),adj=c(0,0)) 信頼区間の導出(リサンプリング使用) 1000個のブートストラップサンプルを使用した サンプル平均の精度推定 12
  8. 中心極限定理(CLT) 平均 、分散  の母集団 から無作為抽出した 個の標本の平均  は 平均 、分散    の正規分布に近づく。    は標準正規分布に従うため、

    Zスコアが算出可能 Zスコアに変換した値が特定の範囲内に存在する確率か ら、信頼区間を導出することが可能! 標本平均は、母集団分布に関係なく正規分布に従う 平均が0、標準偏差 (SD) が1に なるように変換した得点 13
  9. 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) ## [1] 19.23638 22.00819 ブートストラップで得られる信頼区 間 [19.21, 21.989] 中心極限定理から数理的に 得られる信頼区間 [19.23638, 22.00819] qnorm() はZスコアを算出する関数 15
  10. 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) ## [1] 19.20130 22.04327 qt() はt分布の分位点を算出する関数 t分布から数理的に得られる信頼区 間 [19.20130, 22.04327] 17
  11. 仮説検定の手順 ①帰無仮説と対立仮説を設定する。 「帰無仮説(H 0 )」:検証したい仮説 ⇒今回の場合、帰無仮説は、「サンプル間に差がない」こと。 「対立仮説(H 1 )」:帰無仮説に対する仮説 (本来証明したい仮説)

    ⇒今回の場合、対立仮説は、「サンプル間に差がある」こと。 ②帰無仮説の真偽を検証するための統計量 (今回の場合、差があると言える閾値 )を決定する。 ③統計値を計算し、それを参照値 (今回の場合、平均値の差)と比較して、意味のある差(有意差)がある のかをP値を用いて判断する。 ④結果、帰無仮説H 0 を棄却するか棄却しないかを判断する。 23
  12. # シード値を設定 (毎回同じ結果にするための設定 ) 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[0]," :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
  13. 3.2.2 Using t-test for difference of the means between two

    samples t検定による2つのサンプルセット間の平均値の比較 29
  14. 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
  15. 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
  16. 検定の評価方法 TN/(FP+TN): 特異度 TP/(TP+FN): 感度(検定の検出力) 帰無仮説を 採択した 帰無仮説を 棄却した 帰無仮説が

    真であった 帰無仮説が 偽であった (帰無仮説が真である数 ) (対立仮説が真である数 ) (第2種の過誤) (第1種の過誤) 高くすることが目標 37
  17. 多重検定の修正法② ②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
  18. # ライブラリ読み込み library(qvalue) data(hedenfalk) qvalues <- qvalue(hedenfalk$p)$q bonf.pval=p.adjust(hedenfalk$p,method ="bonferroni") fdr.adj.pval=p.adjust(hedenfalk$p,method

    ="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
  19. 変数間の関係を測定/モデル化したい in Genomics 薬を X 投与! 肝臓である遺伝子 発現量が Y 上昇!

    X 箇所がメチル化! メチル化された遺伝 子発現量が Y 減少! X と Y の関係をモデル化したい! 43
  20. βの値によるコスト関数の大きさ コスト関数/損失関数を用いたアプローチ 「コスト関数」もしくは「損失関数」と呼ばれる関数を最適化する。 一般的にこの関数には、予測値 と観測値 の二乗誤差が使用される。 最適化手順(最急降下法) 1. ランダムな開始点( ,

    )を選択する。 2. コスト関数の偏導関数から、コスト関数をどの方向に最適化するかを決定。 3. コスト関数を最小化する方向に向かって最適化を進める。ステップの大きさはパラ メータとして与える。 4. 収束するまで手順2、3を繰り返す。 48
  21. 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
  22. 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
  23. カテゴリ変数を用いた回帰 # 線形回帰結果を取得 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
  24. 線形回帰における落とし穴 非線形性 真の関係性が線形からかけ離れている場合、予測精度が低下する 場合によって , , 等の変換を検討する必要がある。 説明変数の相関 (多重共線性) 説明変数間の相関が高い場合、係数を適切に推定できない。

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