研究室内で、「分類群ごとの積み上げグラフをどうやって作るか」という質問があったので簡単な作り方の例を紹介します。
（こういうやつ↓）
Networks Depicting the Fine-Scale Co-Occurrences of Fungi in Soil Horizons | PLOS ONE
用意するファイルは以下の3つです。
- サンプル×種の表 (行がサンプル、列が種、中身は個体数情報など）
- 分類群の情報 ( 行が分類群の最小単位、種でもいいしNGSデータならOTU・ASVです。列が最小単位の分類群の上位階層）
- あれば環境データ(行がサンプル、列が環境データ）
テストデータです
site=10
sp=5
# サンプル×種の表
mat = matrix(runif(site*sp), nrow=site, ncol=sp,
dimnames=list(paste('S', formatC(1:site,width=2,flag="0"), sep='_'),
paste('Sp', formatC(1:sp,width=2,flag="0"), sep='_')))
# 分類群の情報
taxa = data.frame('Phylum' = formatC(as.integer(runif(sp, 0, 2)),width=2,flag="0"),
'Order' = formatC(as.integer(runif(sp, 0, 4)),width=2,flag="0"),
'Genus' = formatC(as.integer(runif(sp, 0, 8)),width=2,flag="0"),
row.names = colnames(mat))
for(i in 1:ncol(taxa)){ taxa[,i] <- paste(colnames(taxa)[i], taxa[,i], sep='_')}
# 環境データ
env = data.frame('Sample' = rownames(mat),
'Env1' = as.integer(runif(site, 0, 2)),
'Env2' = as.integer(runif(site, 0, 4)),
row.names = rownames(mat))サンプル×種の表の列が、最小分類群単位になっているので、これを見てみたい分類群の階層で、足し合わせます。(例えば、Sp_01 から Sp_05 は、同じ Order_03 という目に入るので、それらを足してあげます）
ただし、これを行うのは、種数が多いと大変です。
なので、関数を用意しました。
# x がサンプル×種の表
# y が分類軍の情報
# taxaLabel がまとめたい分類群の名前
# function は足しわせるのに都合が悪い時用に、用意してあります。
# 注意点は x のcolnames と y のrownamesが対応してある必要があります
Taxa.mat <- function(x, y, taxaLabel, func=function(x){sum(x)}){
colnames(x) <- y[colnames(x), taxaLabel]
summary <- do.call(cbind,
lapply(unique(colnames(x)),
function(a){ num <- which(colnames(x)==a)
apply(as.matrix(x[,num]), 1, func)}) )
colnames(summary) <- unique(colnames(x))
summary
}
あとは、ggplot の geom_bar にかけてあげるだけできます。
library(ggplot2)
library(tidyr)
merge.mat <- Taxa.mat(mat, taxa, 'Genus')
# ロングフォーマットになおす。その際に環境データをくっつける
lf <- gather(cbind(env, merge.mat), key, value, -c(1:ncol(env)))
head(lf)
# bar1 がただの積み上げ。geom_bar の引数のposition = 'fill を指定すると、100 % 積み上げグラフになります
bar1 = ggplot(lf)+
geom_bar(aes(x=Sample, y=value, fill=key), stat='identity')
bar2 = ggplot(lf)+
geom_bar(aes(x=Sample, y=value, fill=key), stat='identity', position='fill')
# ｘ軸を環境データ(Env2 などにしたら、さらに環境でまとめられます）
bar3 = ggplot(lf)+
geom_bar(aes(x=Env2, y=value, fill=key), stat='identity', position='fill')
左がbar1 の結果、右がbar2 の結果x軸ラベルの重なりなどは、以前の記事を参考にしてもらえればと思います。
ggplot2 で、軸ラベルの重なりを防ぐ方法 - 微生物屋のノート
以上が一例になります。
この記事を書きながら、ggplot の段階で色分けする方法もあるなと気づいていしまったので、今回の記事はやりやすい方法を見つける参考にしていただければと思います。
