以前のブログで DADA2 と Claident を用いたアンプリコンシーケンスのデータ解析について触れました。今回は DADA2 と Claident を用いて「分類群 × サンプルのテーブル」ができた後の解析について解説します。DADA2 は Amplicon Sequence Variant (ASV) を生成するので、ここではこの「分類群 × サンプルのテーブル」を ASV テーブルと呼ぶことにします。分類群 = Operational Taxonomic Unit (OTU) の場合は OTU テーブルと呼びます。
DADA2 後の統計解析
DADA2 で ASV テーブルを作成した後はいくつかお決まりの解析が待っています。ここでは以下の 5 点について、R と phyloseq パッケージを利用した実装とともに解説します。
- 群集組成バターンを棒グラフで可視化
- 種多様性パターンの可視化
- Rarefaction カーブの作成
- 次元圧縮によるパターンの可視化
- ASV の再クラスタリング
phyloseq は McMurdie & Holmes (2013) PLoS ONE で発表されたアンプリコンシーケンスのデータの統計解析を行うための R のパッケージです。2021年4月時点で5000回以上引用されているようです。公式チュートリアルはこちら。
デモデータ
今回も Ushio (2019) Methods in Ecology and Evolution のデータを使います。全配列データは DRA (こことここ) に、全解析コードやメタデータは Zenodo (ここ) にアーカイブしてありますが、今回のデモでは一部のデータだけ使用します。
研究は DNA 抽出法の検討をするために行いました。複数のサイトから環境水をサンプリングし、16S メタバーコーディングを行っています。DNA 抽出法は 3 種類試しました。デモ用のデータは、２サイト (Site = “sea”, “pond”)、３抽出法 (Method = “w_beads”, “wo_beads”, “destruction”)、５反復 + ネガコン のデータを含みます。
DADA2 や Claident の部分は飛ばして、ASV table からスタートします。Github にデモデータやスクリプトを置いたので、自分で再現したい方はダウンロードしてください (https://github.com/ong8181/random-scripts/tree/master/05_DADA2Phyloseq)。
準備
まずは以下のコードで必要なパッケージをロードします。
# Load library and functions
library(tidyverse); packageVersion("tidyverse") # 1.3.0
library(phyloseq); packageVersion("phyloseq") # 1.34.0
library(cowplot); packageVersion("cowplot") # 1.1.1
library(ggsci); packageVersion("ggsci") # 2.9
theme_set(theme_cowplot())
tidyverse はデータ整形を効率的に行うため、cowplot と ggsci は可視化のための ggplot2 パッケージの補助パッケージです。ggplot2 は tidyverse をロードするときに同時にロードされています。theme_set(theme_cowplot()) で図のテーマを指定しています (これはただの好みです)。
パッケージのロードができたら以下のコードでサンプルメタデータ、ASV テーブル、 分類群情報を読み込みます。分類群情報は Claident で推定したものです。
sample_sheet <- read.csv("data_subset/sample_sheet.csv", row.names = 1)
asv_sheet <- read.csv("data_subset/asv_sheet.csv", row.names = 1)
tax_sheet <- read.csv("data_subset/tax_sheet.csv", row.names = 1)
phyloseq へのデータの読み込み
デモデータの読み込みができたら、まずはデータがきちんと phyloseq オブジェクトに変換できる形になっているかどうかを確認します。
dim(sample_sheet); dim(asv_sheet); dim(tax_sheet)
all(rownames(asv_sheet) == rownames(sample_sheet))
all(colnames(asv_sheet) == rownames(tax_sheet))
以下の条件を満たしている必要があります:
asv_sheet の行数 = サンプル数 = sample_sheet の行数asv_sheet の行名 = サンプル ID = sample_sheet の行名asv_sheet の列数 = 分類群数 = tax_sheet の行数asv_sheet の列名 = 分類群 ID = tax_sheet の行名
今回は ASV テーブルの行 = サンプル、列 = 分類群、となっていますが、行と列を逆にしたテーブルも読み込み可能です。データが phyloseq オブジェクトに変換できる形になっていることを確認できたら、以下のコードで phyloseq オブジェクトを生成します。
ps_all <- phyloseq(otu_table(asv_sheet, taxa_are_rows = FALSE),
sample_data(sample_sheet),
tax_table(as.matrix(tax_sheet)))
ASV テーブルの行 = 分類群、列 = サンプル、となっている場合は taxa_are_rows = TRUE としてください。tax_sheet は as.matrix() をかまさないと読み込めないのですが、理由は分かりません。また、系統樹情報も同時に読み込むことができますが、今回は利用しません。
phyloseq にはいくつものデータハンドリングのための関数が用意されています。
## Select "sea" samples
subset_samples(ps_all, Site == "sea")
## Select samples with total reads > 10^4
prune_samples(sample_sums(ps_all) > 10^4, ps_all)
## Select "Proteobacteria" phylum
subset_taxa(ps_all, phylum == "Proteobacteria")
## Select top 50 taxa
top_taxa <- names(sort(taxa_sums(ps_all), decreasing = TRUE)[1:50])
prune_taxa(top_taxa, ps_all)
subset_xxx() と prune_xxx() が似ていますが、選びたいサンプル名や分類群名のリストがある場合は prune_xxx() が便利です。例えば、select_samples <- c("R001", "R002", "R003") とある場合は prune_samples(select_samples, ps_all) として選べます。ps_all を入れるべき位置が subset_xxx() と prune_xxx() で違う理由は分かりません。
## Transform sample counts
transform_sample_counts(ps_all, function(x) 100* x / sum(x))
transform_sample_counts(ps_all, log)
上では配列を相対値に変換、log を取る、という操作をしています。
全体のパターンを棒グラフで可視化
群集組成の可視化
plot_bar() で全体のパターンを確認できます。また、 + を続けることで ggplot2 の機能が利用できます。scale_fill_igv() は ggsci パッケージ内のカラーパレットです。
# plot_bar
plot_bar(ps_all, x = "Sample", fill = "phylum") + scale_fill_igv()
facet_wrap() を使って、サイトと手法ごとに分けて結果を表示。
plot_bar(ps_all, x = "replicate", fill = "phylum") +
facet_wrap(.~ Site + Method) + scale_fill_igv()
種多様性パターンの可視化
plot_richness() で多様性パターンを可視化できます。以下ではシンプルに検出 ASV の数をプロットしています。
# plot_richness
plot_richness(ps_all, measures = "Observed")
measures の引数を変更するとシンプソンやシャノンの多様度指数も図示することができます。
psmelt() と ggplot() による可視化
psmelt() を使うと元データをバラすことができ、より自由度の高い可視化が可能になります。例えば、以下は最も優占していた分類群の配列数を箱ひげ図 + ジッタープロットで表示したものです。
# Visualize patterns (scatterplot)
ps_m1 <- psmelt(ps_all) %>%
select(OTU, Sample, Abundance, Site, Method) %>%
filter(OTU == top_taxa[1], Site == "pond")
ggplot(ps_m1, aes(y = Abundance, x = Method)) +
geom_boxplot(outlier.shape = NA, width = 0.2) +
geom_jitter(height = 0, width = 0.2)
Rarefaction カーブ
｢自分の配列データが多様性をどの程度を捉えているのか｣ を配列数が増える毎にどのくらい検出 ASV の数が増加していくのかで把握します。この「検出 ASV 数 v.s. 配列数」の曲線は rarefaction カーブと呼ばれ、微生物群集関連の研究では頻繁に登場します。以下では、ここの記事を参考に rarefaction カーブを ggplot で書いています。
まず、群集生態学関連のパッケージ vegan でデータを rarefy しています。ちなみにこの時点でも図が出力されます。これ以後のコードは単に ggplot で描くためだけのものです。
## From
# Rarefy samples
ps_rr <- vegan::rarecurve(as.matrix(otu_table(ps_all)), step = 50, label = TRUE)
rarecurve で得た出力 (ps_rr) を ggplot 向けに整形します。
rare <- lapply(ps_rr, function(x){
b <- as.data.frame(x)
b <- data.frame(ASV = b[,1], raw_read = rownames(b))
b$raw_read <- as.numeric(gsub("N", "", b$raw_read))
return(b)
})
# Label list
names(rare) <- rownames(sample_data(ps_all))
# Convert to the data frame
rare_df <- map_dfr(rare, function(x) return(data.frame(x)), .id = "sample")
最後に ggplot で図示します。概ね 5,000 – 10,000 配列以上取得しても出現する ASV の数に変化はないようです。ただし、DADA2 による解析結果は原理上シングルトンが出現しないので解釈には注意が必要です。
# Visualize with ggplot
ggplot(rare_df, aes(x = raw_read, y = ASV, color = sample)) +
geom_line() + scale_color_igv() +
xlab("Reads") + ylab("The number of ASV")
次元圧縮によるパターンの可視化
次元圧縮 (次元削減) による高次元データの可視化は全体のパターンを把握する上で非常に有効です。以下では、NMDS、t-SNE、UMAP を適用して微生物群集データを 2 次元で可視化しています (各手法の原理の詳細はここでは解説しません)。
ここではデータのうち海から取得した配列データ (Site == "sea" のサンプル) に対して次元圧縮を適用してみます。また、ネガコンを除いたサンプルのみ (sample_nc == "sample") を解析します。
ps_sea <- subset_samples(ps_all, Site == "sea" & sample_nc == "sample")
NMDS
Nonmetric Multidimensional Scaling (NMDS; 非計量多次元尺度法) は次元圧縮の古典的な手法です。phyloseq では ordinate() と plot_ordination() で解析を実行できます。ランダムな過程を含むので、再現性確保のために乱数のシードを設定しています。
# Nonmetric multidimensional scaling
set.seed(1234)
ps_bray <- ordinate(ps_sea, "NMDS", "bray")
plot_ordination(ps_sea, ps_bray, color = "Method", shape = "Method") +
stat_ellipse(geom = "polygon", alpha = 0.1, aes(fill=Method)) +
geom_point(size = 2) + scale_color_startrek() + scale_fill_startrek() +
xlab("Axis 1") + ylab("Axis 2") + ggtitle("NMDS")
上のコードでは stat_ellipse() を使って各グループの信頼区間 (範囲?) も表示しています。scale_xxxx_startrek() は ggsci パッケージの関数で、カラーパレットを指定しています。主成分分析 (PCA) なども ordinate() 関数を用いた同じコードで実行できます。
t-SNE
t-distributed Stochastic Neighbor Embedding (t-SNE) は確率分布に基づいて次元圧縮する手法で 2008年に提案されています (Maaten & Hinton 2008)。アンプリコンシーケンスの微生物群集データよりは RNA-seq データの解析などでよく見る印象です。phyloseq そのものには t-SNE は実装されていませんが、tsnemicroboita パッケージを使うことで phyloseq オブジェクトのまま t-SNE を実行できます (単に t-SNE を実行するだけのパッケージは他にも色々あります)。
# t-SNE
#library(tsnemicrobiota)
tsne_res <- tsnemicrobiota::tsne_phyloseq(ps_sea, distance='bray', perplexity = 5, rng_seed = 1234)
tsnemicrobiota::plot_tsne_phyloseq(ps_sea, tsne_res, color = "Method", shape = "Method") +
stat_ellipse(geom = "polygon", alpha = 0.1, aes(fill=Method)) +
geom_point(size = 2) + scale_color_startrek() + scale_fill_startrek() +
xlab("Axis 1") + ylab("Axis 2") + ggtitle("t-SNE")
UMAP
Uniform Manifold Approximation (UMAP) は 2018 年に提案された手法です (Mcinnes et al. 2018 arXiv)。t-SNE よりもデータをきれいに分けることができ、かつ高速とされています。ただし、今回のデモデータでは分類群の数もサンプルの数も非常に少ないので、その効果は実感できません。アンプリコンシーケンスの微生物群集データに適用されているのはほとんど見たことがありません (自分が知らないだけかもですが)。phyloseq オブジェクトをそのまま扱って UMAP を実行するパッケージは見つけられませんでした。
# UMAP by uwot package
#library(uwot)
set.seed(1234)
umap_res <- uwot::tumap(otu_table(ps_sea)@.Data, n_neighbors = 5, n_components = 2)
umap_df <- cbind(data.frame(sample_data(ps_sea)), umap_res)
colnames(umap_df)[(ncol(umap_df)-1):ncol(umap_df)] <- c("UMAP1", "UMAP2")
ggplot(umap_df, aes(x = UMAP1, y = UMAP2, color = Method, shape = Method)) +
stat_ellipse(geom = "polygon", alpha = 0.1, aes(fill=Method)) +
geom_point(size = 2) + scale_color_startrek() + scale_fill_startrek() +
xlab("Axis 1") + ylab("Axis 2") + ggtitle("UMAP")
ざっと並べてみると、サンプル数や ASV 数がそれほど多くないのであれば NMDS で十分な印象を受けます。データサイズが大きくなれば (RNA-seq のデータとか)、t-SNE や UMAP が良いかもしれません。いずれにしてもどちらも試してみることをおすすめします。
ASV の再クラスタリング
ASV はアンプリコン同士の違いを 1 塩基レベルで識別する非常に強力な方法ですが、種内多型のような変異も拾うため、多様性 (= ASV の数) が増大し、ときに見たいパターンが見えなくなるときがあります。そのような場合、DADA2 でエラー除去した ASV を OTU として再クラスタリングすることが有効な場合があります。
以下では、R を用いて DADA2 の出力を OTU に再クラスタリングする方法を紹介します。以下のコードは DADA2 の issues のページで議論されていたものを参考にしています。Biostrings と DECIPHER (http://www2.decipher.codes/) というパッケージを使用しています。また、R 内で解析を完結されることにこだわらないのであれば vsearch などのツールを使っても同様のことができるはずです。
クラスタリングプロセスは非常にシンプルです。以下では 97% 一致で OTU を作成しています。
# Preparation
n_cores <- parallel::detectCores() # set to number of cpus/processors to use for the clustering
dna <- Biostrings::DNAStringSet(tax_sheet$seq)
# Find clusters of ASVs to form the new OTUs
aln <- DECIPHER::AlignSeqs(dna, processors = n_cores)
aln_dist <- DECIPHER::DistanceMatrix(aln, processors = n_cores)
# use `cutoff = 0.03` for a 97% OTU
clusters <- DECIPHER::IdClusters(aln_dist, method = "complete",
processors = n_cores, type = "clusters",
cutoff = 0.03, showPlot = F)
クラスタリングができたらクラスター情報を利用して ASV テーブルを変換します。
colnames(clusters) <- "cluster_id"
clusters$cluster_id <- factor(clusters$cluster_id, levels = unique(clusters$cluster_id))
length(unique(clusters$cluster_id))
# Use dplyr to merge the columns of the asv_sheet matrix for ASVs in the same OTU
merged_asv_sheet <- asv_sheet %>% t %>%
rowsum(clusters$cluster_id) %>% t
また、phyloseq オブジェクトにクラスタリング結果を適用することもできます。speedyseq パッケージ内の merge_taxa_vec() という関数を使用しています。
# Make OTU-based phyloseq objects
# Import DADA2 ASV output to phyloseq
asv_sheet2 <- asv_sheet; colnames(asv_sheet2) <- tax_sheet$seq
ps_tmp <- phyloseq(otu_table(asv_sheet2, taxa_are_rows = FALSE))
# Merge taxa in "ps" using cluster information
ps_otu <- speedyseq::merge_taxa_vec(ps_tmp, group = clusters$cluster_id)
以下は解析そのものとは関係ないですが、ASV と OTU の対応を確認するためのサマリーを作成しています。
# Get OTU sequences and OTU table
otu_seqs <- colnames(otu_table(ps_otu))
otu_sheet <- otu_table(ps_otu)@.Data
# Sanity check (make clustering summary data frame)
otu_only <- data.frame(taxa_id = sprintf("OTU%05d", 1:length(otu_seqs)),
seq = otu_seqs)
clusters$seq_sum <- colSums(asv_sheet)
clusters$asv_seq <- tax_sheet$seq
asv_1st <- clusters %>% group_by(cluster_id) %>% summarize(otu_seqs = asv_seq[[1]])
clusters$otu_seq <- unlist(asv_1st[match(clusters$cluster_id, asv_1st$cluster_id), "otu_seqs"])
お役立ちパッケージ
その他、phyloseq にはいろいろな拡張パッケージがあります。一部を紹介します。
#library(speedyseq)
system.time(psmelt(ps_all)) # 0.103 sec elapsed
system.time(speedyseq::psmelt(ps_all)) # 0.025 sec elapsed
#library(fantaxtic)
ps_top20 <- fantaxtic::get_top_taxa(ps_all, n = 20, relative = TRUE,
discard_other = FALSE, other_label = "Other")
fantaxtic::fantaxtic_bar(ps_top20, color_by = "phylum",
label_by = "phylum", other_label = "Other", order_alg = "alph")
まとめ
微生物群集のアンプリコンシーケンス解析では、サンプルメタデータ･配列データ･ASV (or OTU) テーブル･分類群データ･系統樹など、複数のデータを組み合わせて解析することが求められます。phyloseq はそのようなデータを整然と扱えるようにデザインされていますので、うまく使って再現性のある解析を進めていただければと思います。
References
Written with StackEdit.
