ELISAの測定結果をRで計算する 2. 4パラメーターロジスティック回帰で吸光度から濃度を求める
0. はじめに
それでは、4パラメーターロジスティック曲線を描いて、回帰式からELISAの値を算出していきます。前回記事の続きです。
pediatricsurgery.hatenadiary.jp
この記事のコンテンツは
となります。尚、本記事は前ブログの記事を大幅加筆修正したものです。前ブログでの誤りなども記載しておりますので、ご参考いただければと思います。
1. 必要なパッケージのインストール
Rで4パラメーターロジスティック曲線を描く方法はいくつかあるようですが、{nplr}と{drc}を使用します
まずdrcをダウンロードしてインストールしておきます。ついでに{ggplot2}と{dplyr}や{tidyr}を使うため、これらが1つの包括された{tidyverse}も立ち上げておきましょう。{ggplot2}と{dplyr}とは何ぞ？という人はこちらもどうぞ。データ解析には必須のパッケージです。
http://stat.biopapyrus.net/ggplot/
# install.packages("drc", dependencies = TRUE) library(drc) # install.packages("nplr") library(nplr) # install.packages("tidyverse") library(tidyverse)
2. ELISAデータをGitHubから入手する
ここでELISAのデータを使用したいのですが、ちょうどいいデータがなかったので、コチラのデータを使用して下さい。練習で使った96wellのデータを少しいじったものをGitHubにアップロードしました。GitHubからRに直接ダウンロードするには一手間必要です。ちょっと面倒くさいけど、下記の様な手順を踏む必要があるようです。
# install.packages("RCurl", dependencies = TRUE) library("RCurl") getURL("https://raw.githubusercontent.com/Razumall/NOPS/master/Data/ELISA.csv") %>% read.csv(text = ., header = TRUE) %>% tbl_df() -> df
と入力すると
| X1 | X2 | X3 | X4 | X5 | X6 | X7 | X8 | X9 | X10 | X11 | X12 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1.1197 | 1.0665 | 0.0920 | 0.0932 | 0.1470 | 0.1299 | 0.1355 | 0.1485 | 0.2746 | 0.2875 | 0.0572 | 0.0576 |
| 0.6696 | 0.6609 | 0.0910 | 0.0891 | 0.1700 | 0.1606 | 0.1160 | 0.1205 | 0.3440 | 0.3598 | 0.0579 | 0.0580 |
| 0.3759 | 0.3399 | 0.0880 | 0.0852 | 0.1272 | 0.1170 | 0.0887 | 0.0892 | 0.3433 | 0.3631 | 0.0575 | 0.0576 |
| 0.1898 | 0.1869 | 0.0696 | 0.0678 | 0.1175 | 0.1165 | 0.0667 | 0.0662 | 0.4164 | 0.4589 | 0.0691 | 0.0565 |
| 0.1158 | 0.1110 | 0.0554 | 0.0545 | 0.2451 | 0.2428 | 0.1200 | 0.1226 | 0.2679 | 0.2908 | 0.0567 | 0.0626 |
| 0.0812 | 0.0817 | 0.1431 | 0.1382 | 0.0903 | 0.0879 | 0.1168 | 0.1276 | 0.3529 | 0.3694 | 0.0525 | 0.0561 |
| 0.0682 | 0.0665 | 0.1015 | 0.1008 | 0.1631 | 0.1700 | 0.2103 | 0.2110 | 0.0507 | 0.0504 | 0.0565 | 0.0568 |
| 0.0529 | 0.0529 | 0.0836 | 0.0780 | 0.0693 | 0.0723 | 0.1936 | 0.1885 | 0.0504 | 0.0511 | 0.0572 | 0.0574 |
こんな感じのデータが出力されると思います。8 * 12 = 96 wellになっていて、2列ずつのduplicateになっており、左の2列がstandardになっており、1000, 500, 250, 125, 62.5, 31.25, 15.625, 0 ( blank ) となっています。残りの6列は測定したいサンプルの吸光度を示します。
3. バックグラウンドの値を差し引く
空のwellでも吸光度は0にならないので、何も入れなかったwellの吸光度をバックグラウンドとして全体の値から差し引く必要があります。1列目と2列目の8行目がバックグラウンドにしてあります。
df %>% mutate_all(funs(. - as.numeric((df[8,1]+df[8,2])/2 ))) -> df2 # バックグラウンドの平均をとって全体から差し引く
| X1 | X2 | X3 | X4 | X5 | X6 | X7 | X8 | X9 | X10 | X11 | X12 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 1.0668 | 1.0136 | 0.0391 | 0.0403 | 0.0941 | 0.0770 | 0.0826 | 0.0956 | 0.2217 | 0.2346 | 0.0043 | 0.0047 |
| 0.6167 | 0.6080 | 0.0381 | 0.0362 | 0.1171 | 0.1077 | 0.0631 | 0.0676 | 0.2911 | 0.3069 | 0.0050 | 0.0051 |
| 0.3230 | 0.2870 | 0.0351 | 0.0323 | 0.0743 | 0.0641 | 0.0358 | 0.0363 | 0.2904 | 0.3102 | 0.0046 | 0.0047 |
| 0.1369 | 0.1340 | 0.0167 | 0.0149 | 0.0646 | 0.0636 | 0.0138 | 0.0133 | 0.3635 | 0.4060 | 0.0162 | 0.0036 |
| 0.0629 | 0.0581 | 0.0025 | 0.0016 | 0.1922 | 0.1899 | 0.0671 | 0.0697 | 0.2150 | 0.2379 | 0.0038 | 0.0097 |
| 0.0283 | 0.0288 | 0.0902 | 0.0853 | 0.0374 | 0.0350 | 0.0639 | 0.0747 | 0.3000 | 0.3165 | -0.0004 | 0.0032 |
| 0.0153 | 0.0136 | 0.0486 | 0.0479 | 0.1102 | 0.1171 | 0.1574 | 0.1581 | -0.0022 | -0.0025 | 0.0036 | 0.0039 |
| 0.0000 | 0.0000 | 0.0307 | 0.0251 | 0.0164 | 0.0194 | 0.1407 | 0.1356 | -0.0025 | -0.0018 | 0.0043 | 0.0045 |
4. スタンダードの散布図を描く
僕は最初の2列 ( X1とX2 ) にduplicateでスタンダードを置くようにしています。4パラメーターロジスティック曲線をまず書くのに必要なのはスタンダードの値なので、それを抜き出してデータフレームにしておきます。
df2 %>% dplyr::select(X1, X2) %>% mutate(Concentration = c(1000*0.5^(0:6), 0)) %>% gather(Column, Absorbance, -Concentration) %>% filter(Concentration!=0) %>% # バックグラウンドは除去する dplyr::select(-Column) -> Standard # kable(Standard, align="c")
| Concentration | Absorbance |
|---|---|
| 1000.000 | 1.0668 |
| 500.000 | 0.6167 |
| 250.000 | 0.3230 |
| 125.000 | 0.1369 |
| 62.500 | 0.0629 |
| 31.250 | 0.0283 |
| 15.625 | 0.0153 |
| 1000.000 | 1.0136 |
| 500.000 | 0.6080 |
| 250.000 | 0.2870 |
| 125.000 | 0.1340 |
| 62.500 | 0.0581 |
| 31.250 | 0.0288 |
| 15.625 | 0.0136 |
上記のコードを実行するとStandardには上記の様なデータが格納されていると思います。これを{ggplot2}で散布図にしてみましょう。
ggplot(Standard, aes(x=Concentration, y=Absorbance))+ geom_point()+ scale_x_log10() # X軸をlog scaleにする
するとこんな感じになると思います。確かにシグモイド曲線になっているのがわかると思います。
5. 4パラメーターロジスティック回帰の回帰式を求めて予測値を求める
5-1. {nplr}: おすすめ
{nplr}のvignetteを参考にしてください。nplr()という関数を使用します。
モデルを作成する
# install.packages("nplr") library(nplr) nplr <- nplr( x = Standard$Concentration, y = Standard$Absorbance )
> nplr Instance of class nplr Call: nplr(x = Standard$Concentration, y = Standard$Absorbance) 4-P logistic model Bottom asymptote: 0.006489387 Top asymptote: 1.936195 Inflexion point at (x, y): 2.953461 0.9713423 Goodness of fit: 0.9987979 Weighted Goodness of fit: 0.9999795 Standard error: 0.01331459 0.006026113
実行するとという結果が返って来ます。どうやら自動でパラメーターをいくつに設定するかを計算してくれるようですね。きちんと、ELISA kitのマニュアルに書かれている4-P logistic modelを選択してくれています。しかもGoodness of fit: 0.9987979であり、モデルの適合度も優れていそうです。
回帰曲線をプロットする
plot( nplr, xlab="Concenttration", ylab="Absorbance")
こんな感じで回帰曲線がプロットされます。
測定値を得る
df2 %>% dplyr::select(-c(X1, X2)) %>% # Standardの列を除外する gather(Column, Absorbance) %>% mutate(Sample = c(rep(1:8, 2), rep(9:16, 2), rep(17:24, 2), rep(25:32, 2), rep(33:40, 2)) ) %>% # 各wellに入れたサンプルの番号をつける mutate(Sample = as.factor(Sample)) %>% group_by(Sample) %>% summarise(Mean = mean(Absorbance), SD = sd(Absorbance)) -> Measured # kable(Measured, align="c")
| Sample | Mean | SD |
|---|---|---|
| 1 | 0.03970 | 0.0008485 |
| 2 | 0.03715 | 0.0013435 |
| 3 | 0.03370 | 0.0019799 |
| 4 | 0.01580 | 0.0012728 |
| 5 | 0.00205 | 0.0006364 |
| 6 | 0.08775 | 0.0034648 |
| 7 | 0.04825 | 0.0004950 |
| 8 | 0.02790 | 0.0039598 |
| 9 | 0.08555 | 0.0120915 |
| 10 | 0.11240 | 0.0066468 |
| 11 | 0.06920 | 0.0072125 |
| 12 | 0.06410 | 0.0007071 |
| 13 | 0.19105 | 0.0016263 |
| 14 | 0.03620 | 0.0016971 |
| 15 | 0.11365 | 0.0048790 |
| 16 | 0.01790 | 0.0021213 |
| 17 | 0.08910 | 0.0091924 |
| 18 | 0.06535 | 0.0031820 |
| 19 | 0.03605 | 0.0003536 |
| 20 | 0.01355 | 0.0003536 |
| 21 | 0.06840 | 0.0018385 |
| 22 | 0.06930 | 0.0076368 |
| 23 | 0.15775 | 0.0004950 |
| 24 | 0.13815 | 0.0036062 |
| 25 | 0.22815 | 0.0091217 |
| 26 | 0.29900 | 0.0111723 |
| 27 | 0.30030 | 0.0140007 |
| 28 | 0.38475 | 0.0300520 |
| 29 | 0.22645 | 0.0161927 |
| 30 | 0.30825 | 0.0116673 |
| 31 | -0.00235 | 0.0002121 |
| 32 | -0.00215 | 0.0004950 |
| 33 | 0.00450 | 0.0002828 |
| 34 | 0.00505 | 0.0000707 |
| 35 | 0.00465 | 0.0000707 |
| 36 | 0.00990 | 0.0089095 |
| 37 | 0.00675 | 0.0041719 |
| 38 | 0.00140 | 0.0025456 |
| 39 | 0.00375 | 0.0002121 |
| 40 | 0.00440 | 0.0001414 |
ここで大事なのはサンプルの入れ間違いなどなどによって、duplicateの値が大きく異なっている場合、もしくはサンプルの入れ忘れが発生していることがあります。Meanがマイナスになっているものは少なくともサンプルの入れ忘れでしょうし、SDが大きいものはちょっと再測定を検討した方がいいかもしれません。ここではマイナスになっているものだけ除去しておきます。
Measured %>% filter(Mean > 0) -> Measured2
Concentrationを予測する
{nplr}のgetEstimates()を使用します。
getEstimates(nplr, Measured2$Mean) %>% tbl_df() %>% dplyr::rename(Absorbance=y, Concentration=x) -> Estimates_nplr # kable(Estimates_nplr %>% head(), align="c")
| Absorbance | x.025 | Concentration | x.975 |
|---|---|---|---|
| 0.0397000 | 20.134198 | 43.21722 | 63.56545 |
| 0.0371500 | 17.233111 | 40.66222 | 62.53931 |
| 0.0337000 | 14.231966 | 37.12865 | 58.59670 |
| 0.0158000 | 2.736750 | 16.49146 | 42.45722 |
| 0.0151566 | 2.513367 | 15.62500 | 40.78934 |
| 0.0877500 | 67.788432 | 86.20568 | 103.56121 |
とすると95%信頼区間を含めた吸光度からの濃度の予測値を立ててくれます。これで欲しい値を手に入れることができました。めでたしめでたし。
予測値のプロットを描いて希釈濃度が妥当か検討する
ggplot(Estimates_nplr, aes(x=Concentration, y=Absorbance))+ geom_point()+ scale_x_log10()+ geom_hline(aes(yintercept=0.1), color="gray", linetype="dashed")
とするとこんな感じの曲線が得られます。実際に吸光度から濃度を計算する際に重要なのは、信頼性の高い予測値はシグモイド曲線が直線になっている範囲でしか得られないとされているようですので水平線を引いた、Absorbance < 0.1のサンプルは希釈濃度を調整して再測定した方がよいのかもしれません。
5-2. {drc}
前ブログで記載していた{drc}のdrm()に関しても記載しておきます。Rで用量反応曲線を参考にさせていただきました。drm()で使用できるモデルの一覧が上記にまとまっていましたので引用させていただきます。
drc パッケージ内の関数とモデルの対応を一部下記に示します。
| 関数 | モデル |
|---|---|
| G.4() | 4 パラメータゴンペルツ |
| L.4() | 4 パラメータロジスティック |
| LL.4() | 4 パラメータ log ロジスティック |
| W1.4() | 4 パラメータワイブル 1 |
| W2.4() | 4 パラメータワイブル 2 |
| BC.4() | 4 パラメータ Brain-Cousens |
モデルを作成する
drc <- drm( Concentration ~ Absorbance, data = Standard, fct = L.4() ) # L.4: 4パラメータロジスティック
これを実行すると各パラメーターが算出されます。
> drc
A 'drc' model.
Call:
drm(formula = Concentration ~ Absorbance, data = Standard, fct = L.4())
Coefficients:
b:(Intercept) c:(Intercept) d:(Intercept) e:(Intercept)
-1.257 -558.225 3577.330 1.445
回帰曲線をプロットする
plot(drc, type="all", col="gray", pch=19)
{nplr}の結果と比べてX軸とY軸が逆になっていることに注意してください。
Concentrationを予測する
一般的な回帰式で使うpredict()を使用します。
Measured2 %>% mutate(Concentration = predict(drc, newdata = .[,2] )) %>% # predict()でConcentrationを予測 dplyr::rename(Absorbance = Mean) -> Estimates_drc # kable(Estimates_drc %>% head(), align="c")
| Sample | Absorbance | SD | Concentration |
|---|---|---|---|
| 1 | 0.03970 | 0.0008485 | 45.23250 |
| 2 | 0.03715 | 0.0013435 | 43.58224 |
| 3 | 0.03370 | 0.0019799 | 41.35549 |
| 4 | 0.01580 | 0.0012728 | 29.91191 |
| 5 | 0.00205 | 0.0006364 | 21.24566 |
| 6 | 0.08775 | 0.0034648 | 77.03440 |
とすると、吸光度から濃度が予測されます。
予測値のプロットを描いて希釈濃度が妥当か検討する
ggplot(Estimates_drc, aes(x=Concentration, y=Absorbance))+ geom_point()+ scale_x_log10()+ geom_hline(aes(yintercept=0.1), color="gray", linetype="dashed")
同じ計算をしているはずなのに、{nplr}の結果と微妙に違うのは…ちょっとわかりません。
番外編: 前ブログ記事の訂正
尚、前ブログの記事で
余談ですが、実際はfctの設定は " L.4 " が4パラメーターロジスティックのoptionなのですが、ここでは " LL.4 " の4パラメータlogロジスティックを使っています。" log " がつくのとつかないのでどういう違いがあるのかわかりませんでしたが、
- 参考にしたサイトでもLL.4を使用している
- Softmax Proで自動で計算してくれる、係数に近かったのはlogありの方
- 実際に吸光度から濃度計算した値もlogありの方が信用できる
というところから fct=LL.4のoptionを選択しました。あとは散布図を描けばOKです。
と記載しており、これは誤りだと訂正させてください。 ( 当たり前ですよね、よく調べず、しかも長期間放置してすみませんでした ) 何人かの方からコメントをいただいていましたが、ほとんど前ブログを放置していたので、気づきませんでした。お詫び申し上げます。
