バイオ系wet研究者のRログ (UV–Visスペクトル編)

__🌾本記事は下記の環境に基づいて執筆されました__
- macOS Mojave (ver. 10.14.6)
-- R (ver. 3.6.2)
--- tidyverse (ver. 1.3.0); patchwork (ver. 1.0.0)
 
__🌾本記事の概要__
＃ バイオ系wet研究者が
＃ `R`の`tidyverse`パッケージを利用して
＃ 実際の実験データを解析する

# はじめに
本記事は`R`あるいは`tidyverse`パッケージについて**詳細に解説することを目的としていません**。それらについては，他の多くの記事や書籍がとても丁寧に解説してくれています。

特に自分は[**松村優哉 et al.『RユーザーのためのRStudio [実践] 入門』**](https://gihyo.jp/book/2018/978-4-7741-9853-8)にかなり助けられました。

ただそれらの記事や書籍の多くは，`iris`や`mpg`のような**サンプルデータセット**を利用して解説を行っており，**バイオ系wet研究者 (特に蛋白質関係の分析が専門)** な自分には少しとっつきにくかった (実用までの距離を感じた) ように思います。

そういった視点でネット上の解説記事を眺めてみると，実際の実験データを利用して解説している記事って案外少ないんじゃないかと思います。
やっぱり実験データはおいそれと公開できないからですかね。
(読者層も限定されちゃいますしね 笑)

そこで本記事 (とその続編) では，**私が大学/大学院時代に測定したけどボツった実験データの中でも公開してなんら問題ないデータ**を活用して，研究の中で行っている実作業のログ (に近いもの) を共有できたら良いなと考えています。

「ああ，こういう使い方もできるのね」という感じで，誰かの (特に自分と専門が近い人たちの) モチベーションにつながることがあれば，それはとっても嬉しいなって。

### `R`および`tidyverse`を使い始めたきっかけ
<details><summary>(余談ですので，興味のない方はすっ飛ばしてください)</summary><div>

最も有名なグラフ作成ソフトは何かと訊かれれば，それはもうみんな大好きMS Excel様になるわけですが，正直言って私はMS Excelでグラフを描くのが大の苦手です。特に発表資料に使うような「きれいな」図をMS Excelで作る自信はこれっぽっちもないですし，そんな作業をしたいとはこれっぽっちも思いません。

大学/大学院まではどうしていたかというと，[Igor Pro (WaveMetrics)](https://www.wavemetrics.com/) という有償のソフトウェアを利用していました。実験データの解析 (フィッティングなど) にもこいつは非常に活躍してくれたものです。

が，2019年の3月に無事学位を取得して (隙あらば……) 名残惜しくも大学を去り，企業に勤めることになった今，残念ながら[Igor Pro](https://www.wavemetrics.com/)くんを使うことはできなくなりました。もちろん部署の予算で購入することはできましたが，ライセンスの関係上他の人に解析結果を共有することが難しいなと考え断念しました。

という経緯で，[Igor Pro](https://www.wavemetrics.com/)の代わりになる「きれいな」グラフが描けて，実験の解析もできて，なおかつオープンな環境がないかなと考えていたとき，白羽の矢が立ったのが`R`と`tidyverse`だったわけです。今を時めく`Python`ちゃんにしなかった理由は (自主規制) です。</div>

</details>

### 本記事の構成
本記事は

- 実験の背景説明
- 実験データの説明
- 実験データの解析

で構成されています。『実験の背景説明』は多くの人にとって不要かもしれませんが，ただ書きたかったのです (すみません)。

# 1. 蛋白質のUV–Visスペクトル
### 1.1. 背景を語る
蛋白質関係の分析で最もシンプルな測定といえば，**UV–Visスペクトル測定**ではないでしょうか。
特に**蛋白質濃度を定量する際**には，UV–Visスペクトル測定が広く用いられています。

>フェニル基をもつフェニルアラニン，フェール基をもつチロシンやインドール基をもつトリプトファンも脂肪族側鎖と同様に疎水性で，蛋白質の内部に多く存在する。**蛋白質を定量するときに280 nmの紫外吸収を用いることが多いが，この吸収はほとんどがトリプトファンとチロシンの吸収によるものである。**
>——有坂文雄『バイオサイエンスのための 蛋白質科学入門』より

他にもBCA法なんかで定量したりもできますが，BCA法はどうしても値がブレやすい (実験が下手とも言う) ので，自分はUV–Visスペクトルから濃度定量することが多いです。

UV–Visスペクトルから濃度を求める際には，いわゆる**ランベルト・ベールの法則**を利用します。

$$A_{280} = \epsilon_{280} c l$$

ここで$A_{280}$，$\epsilon_{280}$，$c$，および$l$は，それぞれ波長280 nmにおける吸光度 (Abs.)，波長280 nmにおけるモル吸光係数，モル濃度，および光路長 (通常は1 cm) を表しています。

蛋白質の$\epsilon_{280}$は，蛋白質のアミノ酸配列から容易に推定することが可能です。実験的に求めていたこともありましたが，結局BCA法を経ることになってイマイチなので，ここ数年はもっぱらアミノ酸配列から求めた理論値を使用しています。

ちなみに大学/大学院時代は，$\epsilon_{280}$の推定に[**ProtParam**](https://web.expasy.org/protparam/)をよく利用していました。

ラボによっては蛋白質の濃度定量の際に$A_{280}$だけしか測定しないところもあるようですが，400 nmより長波長側をみると凝集の度合いが確認できたり，大腸菌から持ち込まれた蛋白質以外の低分子 (例えばヘムなど) の吸収が観察されることもあります。なので，自分としては新しく扱う蛋白質の場合は，なるべく広範囲のUV–Visスペクトルを測定するようにしています。

### 1.2. 実験データについて
今回用いる実験データをGoogleスプレッドシートにアップしました。
https://docs.google.com/spreadsheets/d/1cki05r2u6X6dXOn1ziXMqtcKg1o6IwgV9dlD7zqQ0jw/edit?usp=sharing

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/540277/342e3370-36c7-7393-b150-be2957ca5111.png" width=75%>

以降はこのデータが`protein.xlsx`の`Sheet1`に記録されているものとして説明を行います。

`.csv`ファイル等ではなくあえて`.xlsx`ファイルを用いる理由ですが，結局他の実験者からデータを貰う場合はほとんどがMS Excelなんですよね。あと，最低限のデータ加工については，MS Excelで行ってしまう方が早いと考えています。自分が扱う実験データは，基本的にMS Excelで編集できる程度のサイズですからね……。

本データの内容は，**PBS溶液 (pH 7.4) に溶解した濃度未知な某蛋白質のUV–Visスペクトル**です。本蛋白質の$\epsilon_{280}$については，アミノ酸配列から25,900 M<sup>−1</sup> cm<sup>−1</sup>と推定されています。

本UV–Visスペクトルは，UV–Vis分光光度計UV-2600 (Shimadzu) を用いて測定されました。サンプルセルには**光路長1 cm**の石英セルを用いました。

一般に光路長1 cmの分光光度計において信頼できるAbs.の範囲は，おおよそ0.2–2.0だと思われます (もちろん機種によりますが)。今回測定した蛋白質溶液は非常に濃いと予想されたため，**原液を20倍希釈してスペクトルを測定**しました。

### 1.3. 実験データの解析
それでは実験データの解析に移ります。

はじめに，必要なパッケージのインストールと読み込みを行います。

```analysis.R
# パッケージのインストール
install.packages("tidyverse")

# パッケージの読み込み
library(tidyverse)
library(readxl)
```

`tidyverse`には，今後使用する様々なパッケージ (例えば`dplyr`や`ggplot2`など) が含まれています。

`readxl`はExcelファイルを読み込むために使用します。
`readxl`は`tidyverse`に含まれているため個別でのインストールは不要ですが，使用するためには個別に読み込む必要があります。

次に，Excelファイルにまとめた測定データ (`protein.xlsx`) をインポートして`prot_raw`に格納します。

```analysis.R
# 測定データをExcelファイルからインポート
prot_raw <- read_excel("protein.xlsx", sheet = "Sheet1")
```

`prot_raw`の中身はこんな感じです。いい感じに読み込めていますね。

```prot_raw.R
> prot_raw
# A tibble: 521 x 2
   wavelength_nm abs_protein_x20
           <dbl>           <dbl>
 1          240            1.52 
 2          240.           1.39 
 3          241            1.28 
 4          242.           1.17 
 5          242            1.07 
 6          242.           0.979
 7          243            0.898
 8          244.           0.829
 9          244            0.769
10          244.           0.714
# … with 511 more rows
```

ここで，測定の際に蛋白質溶液を20倍希釈したことを考慮し，蛋白質Abs.の実測値である`abs_protein_x20`列を20倍した`abs_protein`列を新たに作成します。

これには`dplyr`パッケージの`mutate()`関数を使用します。

```analysis.R
# 蛋白質のAbs.実測値を20倍した列を追加
prot <- mutate(prot_raw, abs_protein = abs_protein_x20 * 20)
```

生データに対してデータ処理を行ったので，`prot_raw`あらため`prot`というデータフレームに格納しています。

上記の処理は`%>%` (パイプ演算子) を利用すると次のように書けます (慣れましょう)。

```analysis.R
# 蛋白質のAbs.実測値を20倍した列を追加
prot <- prot_raw %>% 
  mutate(abs_protein = abs_protein_x20 * 20)
```

`prot`の中身はこんな感じです。いい感じにデータ編集できました。

```prot.R
> prot
# A tibble: 521 x 3
   wavelength_nm abs_protein_x20 abs_protein
           <dbl>           <dbl>       <dbl>
 1          240            1.52         30.5
 2          240.           1.39         27.8
 3          241            1.28         25.5
 4          242.           1.17         23.3
 5          242            1.07         21.4
 6          242.           0.979        19.6
 7          243            0.898        18.0
 8          244.           0.829        16.6
 9          244            0.769        15.4
10          244.           0.714        14.3
# … with 511 more rows
```

それではいよいよ蛋白質のUV–Visスペクトルをプロットしてみましょう。
グラフの作成には`ggplot2`パッケージを使用します。

シンプルに書くと次のようなコードになります。

```analysis.R
# 蛋白質のUV-Visスペクトルをプロット
prot %>%
  ggplot(mapping = aes(x = wavelength_nm, y = abs_protein)) +
  geom_line()
```

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/540277/71c4716f-a02e-ad4d-063d-4da8cd2a3e3c.png" width=70%>

はい，プロットできました。
ここから少し描画の設定をいじってみます。

`geom_line()`において`colour`および`size`を指定することで，プロットの色および線の太さをそれぞれ変更します。

`scale_x_continuous()`を指定することでx軸の取り方を変更します。
(蛋白質のUV–Visスペクトルでは280 nmの値が重要)

`labs()`で軸名を変更します


```analysis.R
# 蛋白質のUV-Visスペクトルをプロット
prot %>%
  ggplot(mapping = aes(x = wavelength_nm, y = abs_protein)) +
  geom_line(colour = "#00ACC1", size = 1.2) +
  scale_x_continuous(breaks = seq(240, 400, 40), limits = c(240, 400)) +
  labs(x = "Wavelength (nm)", y = "Absorbance")
```

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/540277/c43fdfc0-ec2a-e316-78ac-cf570a1385e4.png" width=70%>

だいぶ良い感じになりました。ただ，`ggplot2`のデフォルトテーマ (グレー背景) だと少しモダン過ぎるように感じるので，`theme_set()`関数を用いて自分なりにテーマを変更します。

```analysis.R
# グラフのテーマ設定 (macOSの場合)
theme_set(
  theme_minimal(base_size = 20, base_family = "Helvetica Neue Bold") +
  theme(
    panel.grid.major = element_line(colour = "#BDBDBD"), 
    panel.grid.minor = element_blank(), 
    axis.title.x = element_text(margin = margin(7, 0, 0, 0)),
    axis.title.y = element_text(margin = margin(0, 8, 0, 0)),
    axis.text = element_text(colour = "#424242"), 
    plot.background = element_rect(fill = "transparent", color = NA),
    legend.title = element_blank(),
    legend.justification = "top"
  )
)

# 蛋白質のUV-Visスペクトルをプロット
prot %>%
  ggplot(mapping = aes(x = wavelength_nm, y = abs_protein)) +
  geom_line(colour = "#00ACC1", size = 1.3) +
  scale_x_continuous(breaks = seq(240, 400, 40), limits = c(240, 400)) +
  labs(x = "Wavelength (nm)", y = "Absorbance")
```

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/540277/6c86be49-a9e8-86b7-2e2a-1bc511a763c5.png" width=70%>

はい，良い感じにグラフが描画できました。もっとクラシカルな表現も可能ですが，自分はいろいろ試行錯誤した結果この設定に落ち着きました (少しモダンな方が社内で目立つという浅はかな考え……)。

themeの設定については細かく説明しません。複雑なことはしていないので，興味のある方は一行ずつ調べてみてください。

ちなみにWindowsの場合はフォント指定の部分のみ少し変わります。

```analysis.R
# グラフのテーマ設定 (Windowsの場合)
windowsFonts(SG = windowsFont("Segoe UI Bold"))        # この行が違う
theme_set(
  theme_minimal(base_size = 20, base_family = "SG") +  # この行も少し違う
  theme(
    panel.grid.major = element_line(colour = "#BDBDBD"), 
    panel.grid.minor = element_blank(), 
    axis.title.x = element_text(margin = margin(7, 0, 0, 0)),
    axis.title.y = element_text(margin = margin(0, 8, 0, 0)),
    axis.text = element_text(colour = "#424242"), 
    plot.background = element_rect(fill = "transparent", color = NA),
    legend.title = element_blank(),
    legend.justification = "top"
  )
)
```

グラフが書けたので`ggsave()`関数を用いて画像として保存しましょう。

```analysis.R
# グラフを保存
ggsave("protein.png",                
  dpi = 300, width = 5.6, height = 4,  
  bg = "transparent")                  
```

ちなみに，本蛋白質の$\epsilon_{280}$は25,900 M<sup>−1</sup> cm<sup>−1</sup>であることを思い出すと，ランベルト・ベールの法則から次のようにして蛋白質濃度を求めることができます。

```analysis.R
# 蛋白質濃度を計算
prot %>% 
  filter(wavelength_nm == 280) %>% 
  mutate(conc_mM = abs_protein / 25900 * 1000)  

# 出力結果
# A tibble: 1 x 4
  wavelength_nm abs_protein_x20 abs_protein conc_mM
          <dbl>           <dbl>       <dbl>   <dbl>
1           280            1.30        26.1    1.01            
```

蛋白質濃度が1.0 mMと求まりましたね。

それでは続いて，低分子化合物のUV–Visスペクトル測定データの解析に移ります。

# 2. 化合物のUV–Visスペクトル (モル吸光係数の決定)
### 2.1. 背景を語る
低分子化合物は蛋白質と異なり，化合物ごとに様々なUV–Visスペクトルを示します。したがって，化合物の濃度定量を分光学的に行うためには，**化合物ごとに固有の$\epsilon$を実験的に決定する必要があります**。

実験的に$\epsilon$を決定するといっても，それほど難しい作業ではありません。

まず，化合物を精密天秤で秤量し，有機溶媒 (基本的にはDMSO) に溶解させます。次に，秤量値と化合物の分子量からこのストック溶液の濃度を決定し，あとは段階希釈したサンプルについてそれぞれUV–Visスペクトルを測定します。得られたUV–Visスペクトルの**ピークトップにおけるAbs.を希釈率から計算した化合物の理論濃度に対してプロットし，線形回帰によって傾きを求めればピークトップの波長における$\epsilon$が求まります**。

出身ラボでは新規の化合物を購入した際は，はじめに100% DMSO中における$\epsilon$を決定する習慣がありました。$\epsilon$を一度決定してしまえば，たとえ吸湿や量的問題から精度に欠ける秤量になったとしても，UV–Visスペクトルを測定すれば濃度定量できるという安心感があるからです。

ただし，溶媒によって$\epsilon$が大きく変わることはざらですし，$\epsilon$決定時の秤量が不正確だと以降の実験がすべてデタラメになってしまいます。あとは購入ロットで少し変わったり……。

### 2.2. 実験データについて
今回用いる実験データをGoogleスプレッドシートにアップしました。
https://docs.google.com/spreadsheets/d/1QDpoPPa0p7IX9-hrqCUXobhNY8QvXVTyfJmyV87FXEo/edit?usp=sharing

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/540277/7f41c1ac-a292-ae35-9367-49b0d340f166.png" width=85%>

以降はこのデータが`drug.xlsx`の`Sheet1`に記録されているものとして説明を行います。

本データの内容は，**DMSOに溶解した濃度既知 (40–140 μM) な某薬剤化合物のUV–Visスペクトル**です。

本UV–Visスペクトルは，UV–Vis分光光度計DU 800 (Beckman Coulter) を用いて測定されました。サンプルセルには**光路長1 cm**の石英セルを用いました。

### 2.3. 実験データの解析
さきほどと同様に，Excelファイルにまとめた測定データ (`drug.xlsx`) をインポートして`prot_raw`に格納します。

```analysis.R
# 実験データをExcelファイルからインポート
drug_raw <- read_excel("drug.xlsx", sheet = "Sheet1")
```

`drug_raw`の中身はこんな感じです。いい感じに読み込めていますね。

```drug_raw.R
> drug_raw
# A tibble: 321 x 7
   wavelength_nm drug_140uM drug_120uM drug_100uM drug_80uM drug_60uM drug_40uM
           <dbl>      <dbl>      <dbl>      <dbl>     <dbl>     <dbl>     <dbl>
 1          240         4.5        4.5        4.5       4.5       4.5       4.5
 2          240.        4.5        4.5        4.5       4.5       4.5       4.5
 3          241         4.5        4.5        4.5       4.5       4.5       4.5
 4          242.        4.5        4.5        4.5       4.5       4.5       4.5
 5          242         4.5        4.5        4.5       4.5       4.5       4.5
 6          242.        4.5        4.5        4.5       4.5       4.5       4.5
 7          243         4.5        4.5        4.5       4.5       4.5       4.5
 8          244.        4.5        4.5        4.5       4.5       4.5       4.5
 9          244         4.5        4.5        4.5       4.5       4.5       4.5
10          244.        4.5        4.5        4.5       4.5       4.5       4.5
# … with 311 more rows
```

次のようにこのまま強引にプロットすることも可能です。

```analysis.R
drug_raw %>% 
  ggplot(mapping = aes(x = wavelength_nm)) +
  geom_line(mapping = aes(y = drug_40uM), colour = "#8D6E63", size = 1.3) +
  geom_line(mapping = aes(y = drug_60uM), colour = "#AB47BC", size = 1.3) +
  geom_line(mapping = aes(y = drug_80uM), colour = "#00ACC1", size = 1.3) +
  geom_line(mapping = aes(y = drug_100uM), colour = "#43A047", size = 1.3) +
  geom_line(mapping = aes(y = drug_120uM), colour = "#FB8C00", size = 1.3) +
  geom_line(mapping = aes(y = drug_140uM), colour = "#EF5350", size = 1.3) +
  labs(x = "Wavelength (nm)", y = "Absorbance")
```

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/540277/4c0361bb-c66e-ba52-bb64-dff2480c334f.png" width=70%>

場合によってはこのように強引にプロットする方が楽かもしれませんが，ちょっとスマートではありませんし，このままだと自動でlabel (凡例) が表示されません。そこで，プロットを行う前に読み込んだデータを**tidy data**に変換する作業を行います。

**tidy data**というと少し難しく感じられますが，今回の場合だと横長データを縦長データに変換するだけです (下図参照)。変換には`dplyr`パッケージの`pivot_longer()`関数を利用します。

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/540277/e9137801-c4ea-7a1f-1a5d-fbf76439da3c.png" width=90%>

```analysis.R
# 実験データをtidy dataに変換
drug_long <- drug_raw %>% 
  pivot_longer(col = -wavelength_nm, names_to = "conc_drug_uM", values_to = "abs_drug")
```

`drug_long`の中身はこんな感じです。うまく縦長データに変換できました。

```drug_long.R
> drug_long
# A tibble: 1,926 x 3
   wavelength_nm conc_drug_uM abs_drug
           <dbl> <chr>           <dbl>
 1          240  drug_140uM        4.5
 2          240  drug_120uM        4.5
 3          240  drug_100uM        4.5
 4          240  drug_80uM         4.5
 5          240  drug_60uM         4.5
 6          240  drug_40uM         4.5
 7          240. drug_140uM        4.5
 8          240. drug_120uM        4.5
 9          240. drug_100uM        4.5
10          240. drug_80uM         4.5
# … with 1,916 more rows
```

`pivot_longer()`については
[『【tidyr】gather?, spread? もう古い。時代はpivot』](https://qiita.com/yanami/items/3775df6c579fd0a2d60c)
を参考にさせていただきました。

次に`dplyr`パッケージの`mutate()`関数を用いて`conc_drug_uM`列を**濃度の値**に変換します。
さらに`arrange()`関数によって`wavelength_nm`および`conc_drug_uM`を昇順でソートします。

```analysis.R
# conc_drug_uM列を対応する濃度の値に変換
drug_long <- drug_long %>%
  mutate(conc_drug_uM = case_when(
    str_detect(conc_drug_uM, fixed("_40")) ~ 40,
    str_detect(conc_drug_uM, fixed("_60")) ~ 60,
    str_detect(conc_drug_uM, fixed("_80")) ~ 80,
    str_detect(conc_drug_uM, fixed("_100")) ~ 100,
    str_detect(conc_drug_uM, fixed("_120")) ~ 120,
    str_detect(conc_drug_uM, fixed("_140")) ~ 140)) %>% 
  arrange(wavelength_nm, conc_drug_uM)
```

これはなんというか，もっとスマートな書き方がありそうですね。
`drug_long`の中身はこのように変換されました。

```drug_long.R
> drug_long
# A tibble: 1,926 x 3
   wavelength_nm conc_drug_uM abs_drug
           <dbl>        <dbl>    <dbl>
 1          240            40      4.5
 2          240            60      4.5
 3          240            80      4.5
 4          240           100      4.5
 5          240           120      4.5
 6          240           140      4.5
 7          240.           40      4.5
 8          240.           60      4.5
 9          240.           80      4.5
10          240.          100      4.5
# … with 1,916 more rows
```

それではtidy dataに変換したデータを用いて，化合物のUV–Visスペクトルをプロットしてみます。

```analysis.R
# 化合物のUV-Visスペクトルをプロット
drug_long %>% 
  ggplot(mapping = aes(x = wavelength_nm, y = abs_drug, colour = factor(conc_drug_uM))) +
  geom_line(size = 1.3) +
  labs(x = "Wavelength (nm)", y = "Absorbance", colour = "Conc. (μM)") +
  theme(legend.title = element_text(size = 16))
```

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/540277/3b12a615-141f-5bfa-a547-bf983696813a.png" width=80%>

labelの表示や配色が自動で行われました。
`ggplot(mapping = aes())`において`colour`に`factor(conc_drug_uM)`を渡していますが，これは`conc_drug_uM`を離散変数と認識させるためです。

グラフからも判るように短波長側は溶媒であるDMSOの影響を強く受けるため，あらためて波長260–380 nmを描画してみましょう。
同時に`scale_colour_manual`を設定し，マニュアルでの配色を行ってみます。

```analysis.R
# 化合物のUV-Visスペクトルをプロット
drug_long %>% 
  ggplot(mapping = aes(x = wavelength_nm, y = abs_drug, colour = factor(conc_drug_uM))) +
  geom_line(size = 1.3) +
  scale_colour_manual(values = c("#8D6E63", "#AB47BC", "#00ACC1",
                                 "#43A047", "#FB8C00", "#EF5350")) +
  scale_x_continuous(breaks = seq(260, 380, 20), limits = c(260, 380)) +
  scale_y_continuous(breaks = seq(0, 1.2, 0.2), limits = c(0, 1.2)) +
  labs(x = "Wavelength (nm)", y = "Absorbance", colour = "Conc. (μM)") +
  theme(legend.title = element_text(size = 16))
```

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/540277/c31b93a6-8cf2-5f9b-837d-9612293b7c74.png" width=80%>

これで満足のいくUV–Visスペクトルが得られました。

それでは続いて，本薬剤化合物の$\epsilon$を求めてみましょう。

まず，各濃度のUV–VisスペクトルにおけるピークトップのAbs.を抽出します。
このとき，DMSOの影響を受ける短波長側 (本データでは240–260 nm) を除いた範囲でピークトップを抽出することに留意してください。

```analysis.R
# 波長260-400 nmにおけるピークトップを抽出
peak_top <- drug_long %>% 
  filter(wavelength_nm >= 260) %>% 
  group_by(conc_drug_uM) %>% 
  slice(which.max(abs_drug))
```

`slice(which.max())`については
[『本日の懺悔：dplyrでsliceをもっと活用しようね』](https://dichika.hateblo.jp/entry/2019/02/04/144333)
を参考にさせていただきました。

`peak_top`の中身を確認すると……

```peak_top.R
> peak_top
# A tibble: 6 x 3
# Groups:   conc_drug_uM [6]
  wavelength_nm conc_drug_uM abs_drug
          <dbl>        <dbl>    <dbl>
1          286.           40    0.32 
2          286.           60    0.470
3          286.           80    0.620
4          286.          100    0.765
5          286.          120    0.932
6          287           140    1.09 
```

各濃度におけるピークトップのAbs.を抽出することができました。
ピークトップは286–287 nm付近のようです。さきほどのグラフを観かえすと，どうやら問題なさそうだなということが判ります。

続いて本データに対して線形回帰を実行し，$\epsilon_{286}$を決定します。
線形回帰には`lm()`関数を用いて次のように行います。

```analysis.R
# 線形回帰によるモル吸光係数 (e286) の決定
fit_e286 <- peak_top %>% 
  lm(abs_drug ~ 0 + conc_drug_uM, data = .)
```
線形回帰の結果を`summary()`関数で確認します。

```fit_e286.R
> summary(fit_e286)

Call:
lm(formula = abs_drug ~ 0 + conc_drug_uM, data = .)

Residuals:
        1         2         3         4         5         6 
 0.009727  0.004990 -0.001047 -0.010383  0.001580  0.001743 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
conc_drug_uM 7.757e-03  2.901e-05   267.4 1.39e-11 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.00684 on 5 degrees of freedom
Multiple R-squared:  0.9999,	Adjusted R-squared:  0.9999 
F-statistic: 7.15e+04 on 1 and 5 DF,  p-value: 1.388e-11
```

以上より，本化合物の$\epsilon_{286}$は7760 M<sup>−1</sup> cm<sup>−1</sup> (R<sup>2</sup> = 0.9999) と求まりました (有効数字3桁)。

続いて線形回帰の結果を可視化してみましょう。
まず`ggplot2`に読み込ませるために，線形回帰の結果を基に関数を定義します。

```analysis.R
# 決定したe286を基に関数を作成
fit_curve_e286 <- function(x) {
  e286 <- coef(fit_e286)[1]
  return(e286 * x)
}
```

上記コードの`fit_e286[1]`は，線形回帰によって決定した$\epsilon_{286}$のことです。
それでは，`ggplot2`で可視化してみましょう。

```analysis.R
# モル吸光係数に関する解析結果をプロット
peak_top %>% 
  ggplot(mapping = aes(x = conc_drug_uM, y = abs_drug)) +
  stat_function(fun = fit_curve_e286, size = 1.3, colour = "black") +
  geom_point(
    shape = 21, color = "black", fill = "#00ACC1",
    size = 5, stroke = 1.5) +
  scale_x_continuous(breaks = seq(0, 180, 40), limits = c(0, 180)) +
  scale_y_continuous(breaks = seq(0, 1.5, 0.4), limits = c(0, 1.5)) +
  labs(x = "Conc. of drug (μM)", y = "Abs. at 286 nm")
```

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/540277/9ccace5e-2ce0-3fc8-d980-c5a1e6235251.png" width=50%>

はい，できました！

グラフの中にテキストで$\epsilon_{286}$の値を追記することも可能ですが，正直`R`で位置を調整するのは骨が折れるので，自分はMS PowerPoint上で済ませることが多いです (邪悪ですみません)。

最後に化合物のUV–Visスペクトルと線形回帰のグラフを並べて表示させてみましょう。
追加のパッケージ`patchwork`をインストール & ロードします。

```analysis.R
# patchworkをインストール & ロード
install.packages("patchwork")
library(patchwork)
```

さきほど作成したUV–Visスペクトルと線形回帰のグラフを，それぞれ`g_drug_spe`および`g_drug_e286`というオブジェクトに格納し，`patchwork`を使って足し合わせます。

```analysis.R
# 化合物のUV-Visスペクトルをプロット
g_drug_spe <- drug_long %>% 
  ggplot(mapping = aes(x = wavelength_nm, y = abs_drug, colour = factor(conc_drug_uM))) +
  geom_line(size = 1.3) +
  geom_vline(xintercept = 286.5, linetype = "dashed", size = 1) +
  scale_colour_manual(values = c("#8D6E63", "#AB47BC", "#00ACC1",
                                 "#43A047", "#FB8C00", "#EF5350")) +
  scale_x_continuous(breaks = seq(260, 380, 20), limits = c(260, 380)) +
  scale_y_continuous(breaks = seq(0, 1.2, 0.2), limits = c(0, 1.2)) +
  labs(x = "Wavelength (nm)", y = "Absorbance", colour = "Conc. (μM)") +
  theme(legend.title = element_text(size = 16))

# モル吸光係数に関する解析結果をプロット
g_drug_e286 <- peak_top %>% 
  ggplot(mapping = aes(x = conc_drug_uM, y = abs_drug)) +
  stat_function(fun = fit_curve_e286, size = 1.3, colour = "black") +
  geom_point(
    shape = 21, color = "black", fill = "#00ACC1",
    size = 5, stroke = 1.5) +
  scale_x_continuous(breaks = seq(0, 180, 40), limits = c(0, 180)) +
  scale_y_continuous(breaks = seq(0, 1.5, 0.4), limits = c(0, 1.5)) +
  labs(x = "Conc. of drug (μM)", y = "Abs. at 286 nm")

# patchworkパッケージを用いてグラフをがっちゃんこ
g_drug_spe + g_drug_e286 +
  plot_annotation(tag_levels = "A")
```

<img src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/540277/1a07da6b-00bd-ebbd-60bc-ee1e0300dc5e.png" width=100%>

おお！ これはなんというかすごい (記事を書くまで`patchwork`を使ったことがなかった人)。
論文で使用するFigureを作成するのに非常に便利なのではないだろうか。積極的に使っていこう……。

これはおまけ程度ですが，`g_drug_spe`については`geom_vline()`を追記し，ピークトップの波長を示す破線を追加しています。

`patchwork`については
[『patchworkを使って複数のggplotを組み合わせる』](https://qiita.com/nozma/items/4512623bea296ccb74ba)
を参考にさせていただきました。

以上で今回はおしまいです。

# 次回予告
次回は**トリプトファン蛍光消光試験による蛋白質–低分子の相互作用解析**を予定しています。