初めに

前回の記事でRの開発環境について、いくつか書きました。

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
# A tibble: 54 x 21
   `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   
# ... with 44 more rows, and 12 more variables: コピー <dbl>, ...11 <dbl>,
#   ...12 <dbl>, t検定 <dbl>, ...14 <chr>, AVERAGE...15 <dbl>, ...16 <dbl>,
#   ...17 <dbl>, ...18 <lgl>, SE...19 <dbl>, ...20 <dbl>, ...21 <dbl>
>
>

とスラッシュにすればいけました。また、"\"のようなダブルバックスラッシュでもいけました。エスケープ文字だと認識してしまうのかな？

このデータを使って、検定していきます。

ちなみに、データの構造はこの図のようになっています。

f:id:wet-to-dry:20200922163041p:plain

一応、僕の使っているデータも置いておくので、練習にでも使ってください。

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") #make columnList
> data = xl[, columnList] #access data
> data[1,1]
# A tibble: 1 x 1
  `Sample Name`
  <chr>        
1 A-1          
>

そして、次に遺伝子名リストを作成し、それぞれの行を抽出します(データ検索の参考、列抽出の参考)。

> geneList = c("Aaa", "Bbb", "Ccc", "Ddd", "Eee")
> gene1 = data[data$"Target Name"=="Aaa",]
> vx = gene1$ΔCT
> vx #遺伝子名"Aaa"の行が抽出されてベクトルで表示されるはず
[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()はコマンドラインに結果が出る代わりに指定したファイルに結果を保存してくれる関数です。

これで結果まで保存できるようになりました。

結果

以上をまとめたものがこちら。

####200922 test3.R#### #変更必要

#install.packages("readxl")

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