## 必要なパッケージ
require(tidyverse)
require(plotly)
require(GGally)Rとは
Rは統計解析のためのフリーソフトである。もう少し正確に言うと、Rはコンピュータ言語の名前であり、パソコンにソフトウェアとしてインストールされたRは、R言語を使うための「環境」です。Rは多くの機能を持ち、統計解析からデータの前処理、データの概観、論文用グラフの作成まで幅広く利用することができる。また、「パッケージ」として配布されている拡張プログラムをインストールすることで、様々な解析を簡単に行うことができる。このため、Rを使いこなすスキルは、農学や生命科学の研究者にとって有用なものとなっている。
Rを使った簡単な計算
Rを使った解析では、基本的にコマンドを逐次入力しながら対話的に解析を進めていく（ただし、実際に解析を行う際には、事前に一連のコマンドをRスクリプトとして用意しておくと便利である。Rスクリプトを実行した方が、部分的な修正ができたり、解析の履歴を確認できたりと便利なので覚えておこう)。
まずは、コマンド入力で簡単な計算をしながら、Rに慣れることから始める。
Rの最も簡単な使い方は、簡単な算術式を入力して答えを出すことである。例えば
3 + 5 * 3## [1] 18得られた結果をもとに次の計算を行いたい場合は、次のように何らかの変数に結果を代入する。
x <- 1 + 2
x## [1] 3代入された値は、変数名を通して別の計算に使用することができる。
x + 5 * x## [1] 18関数を用いて様々な計算を行うことができる。
abs(x)## [1] 3sin(x)## [1] 0.14112atan(x)## [1] 1.249046log(x)## [1] 1.098612log10(x)## [1] 0.4771213もう少し複雑な計算をしてみよう。平均が \(mu\) で分散が \(sigma^2\)の正規分布の確率密度関数は、次のようになる。 \[f(x)=\frac{1}{\sqrt{2\pi\sigma^2}}\exp\left(-\frac{(x-\mu)^2}{2\sigma^2}\right)\] これをRを使って計算してみる。
mu <- 3
s2 <- 2
x <- 5
1 / sqrt(2 * pi * s2) * exp(- (x - mu)^2 / (2 * s2)) ## [1] 0.1037769計算結果の確認のために、dnorm関数で正規分布の確率密度を計算しても、同じ結果が得られます。
dnorm(x, mu, sqrt(s2))## [1] 0.1037769
Fig. 1. Normal distribution with mean 3 and variance 2
簡単な統計量の計算
まず、6つの要素からなるベクトルを作ってみよう。これは、6品種/系統の米の粒の長さをmm単位で測定したデータである（データの出典は後述する）。
length <- c(8.1, 7.7, 8.2, 9.7, 7.1, 7.3) # mm scale
length## [1] 8.1 7.7 8.2 9.7 7.1 7.3同一品種の穀物幅を入力し、長さ・幅比を計算する。
width <- c(3.7, 3.0, 2.9, 2.4, 3.3, 2.5)
ratio <- length / width
ratio## [1] 2.189189 2.566667 2.827586 4.041667 2.151515 2.920000まず、サンプル数を数えます。
length(ratio)## [1] 6比率の平均は、mean関数を用いて計算できる。
mean(ratio)## [1] 2.782771平均は、summary関数を用いて、分位数とともに求めることもできる。
summary(ratio)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 2.152 2.284 2.697 2.783 2.897 4.042分散は、var関数を用いて計算できる。
var(ratio)## [1] 0.48063662変数間の共分散は cov関数で計算できる。
cov(length, width)## [1] -0.1773333相関係数は、cor関数を使って計算できる。
cor(length, width)## [1] -0.3901388外部データの取り込みと解析
Rを自分の研究に使う場合、表計算ソフトで整理したデータを読み込んで分析することがほとんどだと思われる。ここでは、他のソフトで作成したデータをRに読み込んで解析する手順について説明する。例として、イネ遺伝資源のゲノムワイド関連研究(Zhao et al. 2011; Nature Communications 2: 467)で使用されたデータを使用する。このデータは、Rice Diversity ( http://www.ricediversity.org/data/ )からダウンロードできる。ここでは、そのデータのサブセットを使用する。
csv 形式で保存されたファイルを読み込むには、read.csv関数を使用する。
rp <- read.csv("rice.csv", row.names = 1, stringsAsFactors = T)取り込んだデータのサイズやデータの一部を確認するには、次のようにする。
dim(rp)## [1] 413 8head(rp)## Accession.Name Sub.population Seed.length Seed.width Plant.height
## 1 Agostano TEJ 8.064117 3.685183 110.9167
## 3 Ai-Chiao-Hong IND 7.705383 2.951458 143.5000
## 4 NSF-TV 4 AUS 8.237575 2.926367 128.0833
## 5 NSF-TV 5 AROMATIC 9.709375 2.382100 153.7500
## 6 ARC 7229 AUS 7.118183 3.281892 148.3333
## 7 Arias TRJ 7.347255 2.528891 119.6000
## Panicle.length Flag.leaf.length Blast.resistance
## 1 20.48182 28.37500 8
## 3 26.83333 39.00833 4
## 4 23.53333 27.68333 3
## 5 28.96364 30.41667 5
## 6 30.91667 36.90833 3
## 7 24.62500 36.99000 2データの分析
多数の変数について、その分布や変数間の関係を調べるために解析を行いたいことがよくある。しかし、計測データには、実験的な理由による欠損項目が含まれていることが多い。ここでは、csvファイルから読み込んだデータを解析してみることとする。
先程と同じ方法で、結晶粒の長さ幅比とその平均値を計算してみる。
ratio <- rp$Seed.length / rp$Seed.width
mean(ratio)## [1] NA平均値を算出することができず、“NA”という値しか得られない。なぜだろうか？
比率に欠損値（RではNAと表現）が含まれているためである。
ratio[1:14]## [1] 2.188254 2.610704 2.814950 4.075973 2.168927 2.905327 3.229055 2.270683
## [9] 2.430681 NA 3.418122 3.061198 4.024255 NAそのような場合は、計算時にオプションとしてna.rm = Tを指定してください。
mean(ratio, na.rm = T)## [1] 2.752084粒の長さと幅の相関係数を計算してみよう。
cor(rp$Seed.length, rp$Seed.width)## [1] NA結果はNAとなる。これは、前回と同様にデータの欠損が原因である。
欠損値を処理するオプションを指定し、再度計算を試みる。相関計数を計算するcorや、共分散を計算するcovなど、2変数以上の関係を表す統計量を計算するときは、欠測値に対する値をna.rm = Tではなくuse =というオプションで指定することに注意する。これは、変数A, B, Cがあった場合、それぞれの変数について異なるサンプルが欠測している場合があるため、欠測の扱いが、欠測を無視する（na.rm = T）という単純な指定だけでは済まなくなるためである。
cor(rp$Seed.length, rp$Seed.width, use = "pair")## [1] -0.2837094今回は、相関関係の算出に成功する。
必要なデータを抽出して計算する
ここでは、tidyverseの一連のパッケージの一つであるdplyrの関数を用いて、データの抽出、選択、新たな変数の作成、並び替えを行う。dplyrは、表型のデータ操作に特化したパッケージである。
まずは、品種名（Accession.Name）、分集団（Sub.population）、種子長（Seed.length）、種子幅（Seed.width）の変数だけを抜き出す。
rp.seed <- rp %>% select(Accession.Name, Sub.population, Seed.length, Seed.width)
head(rp.seed)## Accession.Name Sub.population Seed.length Seed.width
## 1 Agostano TEJ 8.064117 3.685183
## 3 Ai-Chiao-Hong IND 7.705383 2.951458
## 4 NSF-TV 4 AUS 8.237575 2.926367
## 5 NSF-TV 5 AROMATIC 9.709375 2.382100
## 6 ARC 7229 AUS 7.118183 3.281892
## 7 Arias TRJ 7.347255 2.528891つぎに、このうち、分集団が、TRJとTEJのサンプルだけを抽出する。
rp.seed.jp <- rp.seed %>% filter(Sub.population %in% c("TRJ", "TEJ"))
head(rp.seed.jp)## Accession.Name Sub.population Seed.length Seed.width
## 1 Agostano TEJ 8.064117 3.685183
## 2 Arias TRJ 7.347255 2.528891
## 3 Asse Y Pung TRJ 8.572092 2.654675
## 4 Baber TEJ 7.679875 3.382188
## 5 Baghlani Nangarhar TEJ 7.859000 3.233250
## 6 Basmati 217 TRJ 9.301512 2.311362上の2つの操作は、以下のようにすると1度に行うことができる。
rp.seed.jp <- rp %>% select(Accession.Name, Sub.population, Seed.length, Seed.width) %>%
filter(Sub.population %in% c("TRJ", "TEJ"))
head(rp.seed.jp)## Accession.Name Sub.population Seed.length Seed.width
## 1 Agostano TEJ 8.064117 3.685183
## 2 Arias TRJ 7.347255 2.528891
## 3 Asse Y Pung TRJ 8.572092 2.654675
## 4 Baber TEJ 7.679875 3.382188
## 5 Baghlani Nangarhar TEJ 7.859000 3.233250
## 6 Basmati 217 TRJ 9.301512 2.311362種子の縦横比を計算したratioという変数を表型データに加える。
rp.seed.jp.ratio <- rp.seed.jp %>% mutate(ratio = Seed.length / Seed.width)
head(rp.seed.jp.ratio)## Accession.Name Sub.population Seed.length Seed.width ratio
## 1 Agostano TEJ 8.064117 3.685183 2.188254
## 2 Arias TRJ 7.347255 2.528891 2.905327
## 3 Asse Y Pung TRJ 8.572092 2.654675 3.229055
## 4 Baber TEJ 7.679875 3.382188 2.270683
## 5 Baghlani Nangarhar TEJ 7.859000 3.233250 2.430681
## 6 Basmati 217 TRJ 9.301512 2.311362 4.024255新たに加えた変数ratioについて昇順に並び替える。
rp.seed.jp.ratio.sorted <- rp.seed.jp.ratio %>% arrange(ratio)
head(rp.seed.jp.ratio.sorted)## Accession.Name Sub.population Seed.length Seed.width ratio
## 1 Bombon TEJ 7.369758 3.990450 1.846849
## 2 Deokjeokjodo TEJ 6.427725 3.404208 1.888170
## 3 Eh Ia Chiu TEJ 6.708800 3.474967 1.930608
## 4 Edomen Scented TEJ 6.606578 3.403167 1.941303
## 5 F.R. 13A TEJ 7.403492 3.763700 1.967078
## 6 Triomphe Du Maroc TEJ 7.652392 3.863633 1.980621今度は、ratioについて降順に並び替え、最初から10行目までのデータだけを抽出する。つまり、ratioの大きさが上位10個までの品種のデータが抜き出される。
rp.seed.jp.ratio.sorted.top10 <- rp.seed.jp.ratio %>% arrange(desc(ratio)) %>%
slice(1:10)
rp.seed.jp.ratio.sorted.top10## Accession.Name Sub.population Seed.length Seed.width ratio
## 1 Basmati 217 TRJ 9.301512 2.311362 4.024255
## 2 L-202 TRJ 9.970591 2.495973 3.994671
## 3 C57-5043 TRJ 9.690625 2.468067 3.926403
## 4 H256-76-1-1-1 TRJ 8.986817 2.298658 3.909592
## 5 Della TRJ 9.838209 2.522709 3.899859
## 6 S4542A3-49B-2B12 TRJ 9.047510 2.339260 3.867680
## 7 Kaybonnet TRJ 9.068720 2.358720 3.844763
## 8 Katy TRJ 9.044283 2.356800 3.837527
## 9 Tia Bura TRJ 10.014290 2.620810 3.821067
## 10 R 101 TRJ 10.224109 2.697636 3.790025分集団ごとに、種子の縦横比の平均を計算する。
rp.seed.jp.ratio %>% group_by(Sub.population) %>% summarise(mean(ratio, na.rm = T))## # A tibble: 2 × 2
## Sub.population `mean(ratio, na.rm = T)`
## <fct> <dbl>
## 1 TEJ 2.24
## 2 TRJ 3.05データの可視化
実際に統計解析を行う前に、データを様々な角度から見てみることは非常に重要である。例えば、前述した平均や分散といった統計は要約のための統計であり、平均や分散が似ている変数であっても、観測値の分布が大きく異なる場合がある。そのため、データをよく見て、データの特徴を把握することが重要である。また、分析結果を論文として発表するための準備をする際にも、データの可視化が必要となる。ここでは、様々なデータの可視化手法を学ぶ。
まず、可視化手法の説明の前に、データ中の変数を直接呼び出せるようにする。
attach(rp)attachコマンドにより、data$…という指定しなくてもPlant.heightを指定できるようになる。このコマンドを使わない場合は、data$Plant.heightという変数を指定しなければならない。
まずはヒストグラムを描いてみよう。
hist(Plant.height)
次に、箱ひげ図を描いてみる。
boxplot(Plant.height)
次に、いもち病抵抗性スコア（Blast.resistance）のヒストグラムを描く。
hist(Blast.resistance)
分布がうまく可視化されたように見えるが、「落とし穴」がある。
いもち病に対する抵抗性は、9段階（0〜9）のスコアで表現される。そこで、9段階の各レベルにどれだけのアクセッションが含まれているかをまとめてみる。
t <- table(Blast.resistance)
t## Blast.resistance
## 0 1 2 3 4 5 6 7 8 9
## 3 77 23 34 36 24 39 36 52 61描いたヒストグラムは、クラス全体をうまく表現できていないことがわかる。
上記のように「表」機能を使ってまとめたデータから、棒グラフを描くことができる。
plot(t, xlab = "Blast resistance scores", ylab = "Frequency")
棒グラフは「barplot」関数で描くこともできます。ただし、上で描いた棒グラフとは見た目が異なる。
barplot(t)
円グラフを描いてみる。
pie(t)
ここからは、変数間の関係を見ていこう。
plot(Plant.height, Panicle.length)
回帰分析により、データに直線を当てはめる。
plot(Plant.height, Panicle.length)
abline(lm(Panicle.length ~ Plant.height))
Zhao et al. (2011)のファイルから読み込んだデータには、形質だけでなく、遺伝資源の遺伝的背景のデータも含まれている。両データを合わせて可視化し、遺伝的背景と形質の間にどのような関係があるのかを調べてみる。
Sub.populationという変数は、遺伝資源の遺伝的背景の違いを表す。遺伝的背景と草丈・穂丈の間にどのような関係があるのか、可視化して確認してみる。
pop.id <- as.numeric(Sub.population)
plot(Plant.height, Panicle.length, col = pop.id)
levels(Sub.population)## [1] "ADMIX" "AROMATIC" "AUS" "IND" "TEJ" "TRJ"legend("bottomright", levels(Sub.population), col = 1:nlevels(Sub.population), pch = 1)
遺伝的背景の違いによる数値の違いを箱ひげ図で示してみよう。
boxplot(Plant.height ~ Sub.population, border = 1:nlevels(Sub.population))
考えられるすべての組み合わせについて、3つの変数の関係を散布図に描いてみよう。
x <- data.frame(Plant.height, Panicle.length, Flag.leaf.length)
pairs(x, col = pop.id)
グラフをファイルに出力する
論文やプレゼンテーションでグラフを使用するために、グラフをPDFファイルに出力できると便利です。ここでは、その方法について簡単に説明する。
先ほど描いた図形をfig00.pdfというファイルに出力してみる。
pdf("fig00.pdf")
x <- data.frame(Plant.height, Panicle.length, Flag.leaf.length)
pairs(x, col = pop.id)
dev.off()## quartz_off_screen
## 2上記のコマンドを実行すると、map.pdfというファイルがRの作業ディレクトリに出力される。
pdf関数は出力されるグラフの大きさを指定することができる。この図のように大きなサイズが必要な場合は、明示的にグラフの出力サイズを指定した方がよい。
pdf("fig00_large.pdf", width = 20, height = 20)
x <- data.frame(Plant.height, Panicle.length, Flag.leaf.length)
pairs(x, col = pop.id)
dev.off()## quartz_off_screen
## 2また、複数の図を同じpdfファイルに繰り返し出力すると、複数のページを持つpdfファイルに保存される。似たような図や大量の図を出力する場合は、1つのpdfファイルにまとめると便利である。
pdf("fig00_multi.pdf")
x <- data.frame(Plant.height, Panicle.length, Flag.leaf.length)
pairs(x, col = pop.id)
plot(x$Plant.height, x$Panicle.length)
abline(lm(Panicle.length ~ Plant.height, data = x))
dev.off()## quartz_off_screen
## 2インタラクティブな図を描く
plotlyパッケージを使うと、インタラクティブな図を描くことができる。ここでは、先に描いた図をplotlyの関数plot_lyを用いて描いてみよう。
まずは、草丈（Plant.height）について、ヒストグラムを描く。
plot_ly(x = Plant.height, type = "histogram") # plot_ly関数でプロット## Warning: Ignoring 30 observations横向きのヒストグラムも簡単に描ける。
plot_ly(y = Plant.height, type = "histogram")## Warning: Ignoring 30 observations箱ひげ図を描く。
plot_ly(x = "Plant Height", y = Plant.height, type = "box")## Warning: Ignoring 30 observationsSub.populationごとに箱ひげ図を描く。
plot_ly(y = Plant.height, color = Sub.population, type = "box")## Warning: Ignoring 30 observationsいもち病抵抗性（Blast.resistance）のスコアについて、棒グラフを描く。棒グラフを描く前に、先程と同様に、table関数で頻度を集計しておく。
t <- table(Blast.resistance) # まずは集計
plot_ly(x = names(t), y = t, type = "bar")円グラフを描く。
plot_ly(labels = names(t), values = t, type = "pie")次に、2変量（草丈Plant.heightと穂長Panicle.length）の散布図を描く。
plot_ly(x = Plant.height, y = Panicle.length,
type = "scatter", mode = "markers")## Warning: Ignoring 38 observations分集団（Sub.population）で色を付けするとともに、ポップアップで示される文字列（Hoverと呼ばれる）に品種名（Accession.Name）を加える。これにより、例えば、どの品種がどのような値を取っているのかを確認できる。
plot_ly(x = Plant.height, y = Panicle.length,
color = Sub.population, text = Accession.Name,
type = "scatter", mode = "markers")## Warning: Ignoring 38 observations同じ内容をdata.frameを使って描くときには、以下のようにする。
df <- rp %>% select(Accession.Name, Sub.population,
Plant.height, Panicle.length, Flag.leaf.length) %>%
drop_na() # 予め、欠測があるデータを除いておく
plot_ly(data = df, x = ~Plant.height, y = ~Panicle.length,
color = ~Sub.population, text = ~Accession.Name,
type = "scatter", mode = "markers")最後に、3変量間の関係を確認する。
plot_ly(data = df, x = ~Plant.height, y = ~Panicle.length, z = ~Flag.leaf.length,
color = ~Sub.population, text = ~Accession.Name,
type = "scatter3d", mode = "markers")ggplotを使って論文投稿用の図を描く
にtidyverseに含まれるggplot関数を用いてヒストグラムを作成する。ここでは、先に作ったものと同じになるが、改めてグラフ用のdata.frameを作成しておく。
df <- rp %>% select(Accession.Name, Sub.population,
Plant.height, Panicle.length, Flag.leaf.length) %>%
drop_na() # 予め、欠測があるデータを除いておく
ggplot(data = df, mapping = aes(x = Plant.height)) +
geom_histogram(fill = "blue")## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
草丈（Plant.height）について、箱ひげ図を描く。
ggplot(data = df, mapping = aes(y = Plant.height)) +
geom_boxplot(fill = "lightblue")
分集団（Sub.population）ごとに、箱ひげ図を描く。
ggplot(data = df, mapping = aes(x = Sub.population, y = Plant.height, fill = Sub.population)) +
geom_boxplot()
ggplotでは、バイオリンプロットと呼ばれる図を描くこともできる。箱ひげ図に比べ、分布のバターンをより直感的に確認できる。
ggplot(data = df, mapping = aes(x = Sub.population, y = Plant.height, fill = Sub.population)) +
geom_violin()
以下のようにすると、棒グラフを描くことができる。ここでは、改めて、いもち病抵抗性のデータが含まれたdata.frameを準備してから棒グラフを描いている。
df2 <- rp %>% select(Sub.population, Blast.resistance) %>% drop_na()
ggplot(data = df2, mapping = aes(x = Blast.resistance)) +
geom_bar()
2変数（Plant.heightとPanicle.length）の散布図を描いてみる。
ggplot(data = df, mapping = aes(x = Plant.height, y = Panicle.length, color = Sub.population)) +
geom_point()
ggplotでは、geom_smooth関数を付け加えることにより、回帰直線を簡単に加えることができる。
ggplot(data = df, mapping = aes(x = Plant.height, y = Panicle.length)) +
geom_point() +
geom_smooth(method = lm)## `geom_smooth()` using formula 'y ~ x'
最後に、3つの変数（Plant.height, Panicle.length, Flag.leaf.length）の関係に加え、Sub.populationも含めて、4者の関係を一度にグラフにする。
df3 <- df %>% select(-Accession.Name)
ggpairs(data = df3)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
なお、ggplot関数で描いた図は、ggsaveを使ってファイルに出力する。
df3 <- df %>% select(-Accession.Name)
ggpairs(data = df3)## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggsave("fig01.pdf")## Saving 7 x 5 in image
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.