初めに
前回の記事でRの開発環境について、いくつか書きました。
wet-to-dry.hatenablog.com
wet-to-dry.hatenablog.com
やっとこさ統計解析の本題に入っていきます。Realtime PCRのデータの統計解析です。
t検定だけならばExcelの関数や手計算で何とかいけていたのですが、ANOVAとTukey-Kramer法を使わなければならなくなったので、Rを使って検定していきます。
Tukey-Kramer法やその他の多重検定についてはこのサイトやこのサイトがわかりやすかったです。
試行錯誤
元データはExcelで持っているので、何とかしてExcelからデータを入力したい。
このサイトを見たところ、
csvで書きだしてから処理
コピペしてclipboardから読み込む
Rパッケージを使う
があるらしいですね。
csvファイルが増えると管理が面倒になるので、できるだけ余分なファイルを作らずにやりたいところ。
一番シンプルでスマートなのはRパッケージかな、と思いましたので、一番有名っぽい"readxl"パッケージをinstallしてみます。
>install.packages("openxl")
WARNING: Rtools is required to build R packages but is not currently installed. Please download and install the appropriate version of Rtools before proceeding:
https://cran.rstudio.com/bin/windows/Rtools/
Warning in install.packages :
'lib = "C:/Program Files/R/R-4.0.2/library"' is not writable
Warning in install.packages :
ディレクトリ 'C:\Users\*******' を作成できません、理由は 'Invalid argument' です
Error in install.packages : unable to create ‘C:/Users/********/R/win-library/4.0’
>
いくつかポップアップダイアログが出た後、このようにエラーで止まってしまいます。
これはこのサイトから、管理者としてRStudioを実行すればいけました。
次に詰まったのが、
> library(readxl)
> xl = readxl::read_excel("C:\Users\excelへのPATH\*****.xlsx", sheet = 1)
エラー: ""C:\U" で始まる文字列の中で 8 進文字なしに '\U' が使われています
>
と、readxlパッケージがPATHを読んでくれません。これはこのサイトから、PATH中のバックスラッシュをスラッシュにしなければならないらしく、
> library(readxl)
> xl = readxl::read_excel("C:/Users/excelへのPATH/*****.xlsx", sheet = 1)
New names:
* `` -> ...4
* AVERAGE -> AVERAGE...8
* SE -> SE...9
* `` -> ...10
* `` -> ...12
* ...
> xl
`Sample Name` `Target Name` `Cт Mean` ΔCT ref RQ AVERAGE...7 SE...8 ...9
<chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <lgl>
1 A-1 Gapdh 29.8 0 0 1 1 0 NA
2 A-1 Aaa 26.6 -3.21 -3.03 1.13 1. 0.111 NA
3 A-1 Bbb 26.8 -2.92 -2.55 1.29 1 0.152 NA
4 A-1 Ccc 27.5 -2.28 -1.92 1.28 1. 0.176 NA
5 A-1 Ddd 27.3 -2.44 -1.92 1.44 1.00 0.220 NA
6 A-1 Eee 26.6 -3.19 -2.79 1.32 1 0.169 NA
7 A-2 Gapdh 29.2 0 NA 1 NA NA NA
8 A-2 Aaa 26.5 -2.67 NA 0.779 NA NA NA
9 A-2 Bbb 27.0 -2.20 NA 0.785 NA NA NA
10 A-2 Ccc 27.8 -1.36 NA 0.679 NA NA NA
>
>
とスラッシュにすればいけました。また、"\"のようなダブルバックスラッシュでもいけました。エスケープ文字だと認識してしまうのかな？
このデータを使って、検定していきます。
ちなみに、データの構造はこの図のようになっています。
一応、僕の使っているデータも置いておくので、練習にでも使ってください。
https://kyotounivcoop-my.sharepoint.com/:x:/g/personal/shinn060425_kyotounivcoop_onmicrosoft_com/EQsqmkfJdVFLqlyYIfj1fs0BHiovt7H6F8Awkcdb58JUzQ?e=10fLNF
僕はrealtime PCRの解析においてΔCT値を正規分布していると考えて検定しますので(参考)、まず"サンプル名"、"遺伝子名"、"ΔCT値"の列を抜き出します。
列の抽出の参考はこちら。
> columnList = c("Sample Name", "Target Name", "ΔCT")
> data = xl[, columnList]
> data[1,1]
`Sample Name`
<chr>
1 A-1
>
そして、次に遺伝子名リストを作成し、それぞれの行を抽出します(データ検索の参考、列抽出の参考)。
> geneList = c("Aaa", "Bbb", "Ccc", "Ddd", "Eee")
> gene1 = data[data$"Target Name"=="Aaa",]
> vx = gene1$ΔCT
> vx
[1] -3.21 -2.67 -3.15 -2.79 -2.50 -2.61 -2.83 -2.53 -3.16
>
このデータを使ってTukey-Kramer法で検定を行います。Rを用いたTukey-Kramer法についてはこちらのサイトを参考にしました。
> fx = factor(rep(c("A","B","C"),c(3,3,3)))
> TukeyHSD(aov(vx~fx))
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = vx ~ fx)
$fx
diff lwr upr p adj
B-A 0.3766667 -0.2835209 1.0368542 0.2633039
C-A 0.1700000 -0.4901875 0.8301875 0.7222893
C-B -0.2066667 -0.8668542 0.4535209 0.6256951
>
これを、それぞれの遺伝子で再帰的に解析します。この際、変数(gene1~gene5)を再帰的に作成する必要があったので、このサイトを参考にしました。
> for (i in 1:length(geneList)) {
+ gene = data[data$"Target Name"==geneList[i],]
+ vx = gene$ΔCT
+ fx = factor(rep(c("A","B","C"),c(3,3,3)))
+ assign(paste("gene", i, sep = ""), TukeyHSD(aov(vx~fx)))
+ }
> gene5
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = vx ~ fx)
$fx
diff lwr upr p adj
B-A 3.323333 1.69604402 4.95062264 0.0018627
C-A 1.666667 0.03937736 3.29395598 0.0456456
C-B -1.656667 -3.28395598 -0.02937736 0.0467114
>
何とかできました。
最後に、データをtxtファイルで出力するコードを書きます。普通に出力してもよかったのですが、Tukey-Kramerの出力がlist形式でわかりづらかったのと、行列の縦がそろったまま見やすく出力したかったので、sink()という関数を使いました。sink()の参考はこちら。
> fileName = "Tukey-Kramer-result.txt"
> path = paste(getwd(), fileName, sep="/")
> sink(path, append = T)
> fileName
> for (i in 1:length(geneList)){
+ print(geneList[i])
+ print(get(paste("gene", i, sep = "")))
+ }
> sink()
>
sink()はコマンドラインに結果が出る代わりに指定したファイルに結果を保存してくれる関数です。
これで結果まで保存できるようになりました。
結果
以上をまとめたものがこちら。
library(readxl)
xl = readxl::read_excel("C:/Users/*****/200922 test3.xlsx", sheet = 1)
columnList = c("Sample Name", "Target Name", "ΔCT")
data = xl[, columnList]
geneList = c("Aaa", "Bbb", "Ccc", "Ddd", "Eee")
for (i in 1:length(geneList)) {
gene = data[data$"Target Name"==geneList[i],]
vx = gene$ΔCT
fx = factor(rep(c("A","B","C"),c(3,3,3)))
assign(paste("gene", i, sep = ""), TukeyHSD(aov(vx~fx)))
}
fileName = "Tukey-Kramer-result.txt"
path = paste(getwd(), fileName, sep="/")
sink(path, append = T)
fileName
for (i in 1:length(geneList)){
print(geneList[i])
print(get(paste("gene", i, sep = "")))
}
sink()
その都度変更しなければならないところは#変更必要と書いてあります。
これの出力はこちら
[1] "Tukey-Kramer-result.txt"
[1] "Aaa"
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = vx ~ fx)
$fx
diff lwr upr p adj
B-A 0.3766667 -0.2835209 1.0368542 0.2633039
C-A 0.1700000 -0.4901875 0.8301875 0.7222893
C-B -0.2066667 -0.8668542 0.4535209 0.6256951
[1] "Bbb"
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = vx ~ fx)
$fx
diff lwr upr p adj
B-A 3.066667 1.6294337 4.5038996 0.0014778
C-A 1.953333 0.5161004 3.3905663 0.0138771
C-B -1.113333 -2.5505663 0.3238996 0.1194498
[1] "Ccc"
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = vx ~ fx)
$fx
diff lwr upr p adj
B-A 2.736667 0.7334353 4.7398980 0.0135548
C-A 1.416667 -0.5865647 3.4198980 0.1555501
C-B -1.320000 -3.3232314 0.6832314 0.1877137
[1] "Ddd"
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = vx ~ fx)
$fx
diff lwr upr p adj
B-A 3.470000 1.6473849 5.2926151 0.0026843
C-A 2.173333 0.3507183 3.9959484 0.0246693
C-B -1.296667 -3.1192817 0.5259484 0.1529928
[1] "Eee"
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = vx ~ fx)
$fx
diff lwr upr p adj
B-A 3.323333 1.69604402 4.95062264 0.0018627
C-A 1.666667 0.03937736 3.29395598 0.0456456
C-B -1.656667 -3.28395598 -0.02937736 0.0467114
BbbのA-B, A-C
CccのA-B
DddのA-B, A-C
EeeのA-B, A-C, B-C
でp<0.05の有意差ありとなっていますね。
まあ、予想通りの結果でした。
まとめ
いい感じでできたので満足です。これを応用すれば他の多重検定も簡単にできると思います。
参考
www.slideshare.net
plaza.umin.ac.jp
keita43a.hatenablog.com
qiita.com
stackoverrun.com
amphipod.hatenablog.com
qiita.com
qiita.com
www.takuro-fujita.com
data-science.gr.jp
futabooo.hatenablog.com
riseki.php.xdomain.jp
