Qiime2で生成したデータからデータの可視化に便利なphyloseqオブジェクト作成する。phyloseqパケージを使えばqzaファイルの読み込みが可能である。

library(phyloseq)
library(qiime2R)
library(knitr)
library(tidyverse)

## -- Attaching packages --------------------------------------- tidyverse 1.3.0 --

## √ ggplot2 3.3.0     √ purrr   0.3.4
## √ tibble  3.0.1     √ dplyr   0.8.5
## √ tidyr   1.1.0     √ stringr 1.4.0
## √ readr   1.3.1     √ forcats 0.5.0

## -- Conflicts ------------------------------------------ tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

taxnomyデータを読み込む

今回silvaからのデータベースを用いたため、Kingdom、Phylum、Class、Genus、Speciesは“；”で区切られている。

setwd("C:/Users/akari/Desktop/fastq2")
taxonomy<-read_qza("taxonomy-silva.qza")
#データを確認
head(taxonomy$data)%>%kable()

Feature.ID	Taxon	Consensus
b5f5167bc021a2b0e18f93d765f303b7	Bacteria;Proteobacteria;Alphaproteobacteria;Sphingomonadales;Sphingomonadaceae;Sphingomonas;uncultured bacterium	0.8
cee5e095946d1f90373e37e3275355d6	Archaea;Thaumarchaeota;Soil Crenarchaeotic Group(SCG)	0.9
a45ebb3a1adc92b2a15285afbefa78bf	Bacteria;Chloroflexi;KD4-96;uncultured bacterium	0.9
59cfcf40b213772307e70ee970f231f0	Bacteria;Verrucomicrobia;Spartobacteria;Chthoniobacterales;DA101 soil group;uncultured bacterium	0.8
af6972e30898eb9a809d47a8e94809dc	Bacteria;Acidobacteria;Acidobacteria;Subgroup 6;uncultured bacterium	0.6
a742c9230d6676357d2c79ac0a42c785	Bacteria;Proteobacteria;Betaproteobacteria;Nitrosomonadales;Nitrosomonadaceae;uncultured;uncultured bacterium	0.7

taxonomyデータから、“;”で区切られたバクテリアの名前とfeature.IDが対応するようなテーブルを生成する

taxtable<-taxonomy$data %>% as.tibble() %>% separate(Taxon, sep=";", c("Kingdom","Phylum","Class","Order","Family","Genus","Species"))

head(taxtable)%>%kable(format ="markdown")

Feature.ID	Kingdom	Phylum	Class	Order	Family	Genus	Species	Consensus
b5f5167bc021a2b0e18f93d765f303b7	Bacteria	Proteobacteria	Alphaproteobacteria	Sphingomonadales	Sphingomonadaceae	Sphingomonas	uncultured bacterium	0.8
cee5e095946d1f90373e37e3275355d6	Archaea	Thaumarchaeota	Soil Crenarchaeotic Group(SCG)	NA	NA	NA	NA	0.9
a45ebb3a1adc92b2a15285afbefa78bf	Bacteria	Chloroflexi	KD4-96	uncultured bacterium	NA	NA	NA	0.9
59cfcf40b213772307e70ee970f231f0	Bacteria	Verrucomicrobia	Spartobacteria	Chthoniobacterales	DA101 soil group	uncultured bacterium	NA	0.8
af6972e30898eb9a809d47a8e94809dc	Bacteria	Acidobacteria	Acidobacteria	Subgroup 6	uncultured bacterium	NA	NA	0.6
a742c9230d6676357d2c79ac0a42c785	Bacteria	Proteobacteria	Betaproteobacteria	Nitrosomonadales	Nitrosomonadaceae	uncultured	uncultured bacterium	0.7

rooted-tree（系統樹）データを読み込む

Qiime2で生成したrooted-tree.qzaも読み込む。（rooted-treeが無くても、phyloseqオブジェクトの作成は可能である。relative abundance等の一部のデータは可視化できるが、系統樹は作成できない。）

setwd("C:/Users/akari/Desktop/fastq2")
tree<-read_qza("rooted-tree.qza")

サンプルOTUデータを読み込む(rarefy済みデータ）

既にrarefy済みのOTUテーブルを使う。Rを使ったrarefyの方法については後述する。

setwd("C:/Users/akari/Desktop/fastq2/core-metrics-results")
otu <- read_qza("rarefied_table.qza")
head(otu$data[1:5,1:5])

##                                  sample04 sample05 sample06 sample10 sample11
## b5f5167bc021a2b0e18f93d765f303b7      459      657      856      826      978
## cee5e095946d1f90373e37e3275355d6      481      295      347      433      346
## a45ebb3a1adc92b2a15285afbefa78bf      381      362      396      366      320
## 59cfcf40b213772307e70ee970f231f0      500      287      362      324      336
## af6972e30898eb9a809d47a8e94809dc      311      482      308      331      304

metadataを読み込む

Qiime2で使用したmetadataを読み込む

setwd("C:/Users/akari/Desktop/fastq2")
metadata <-read.table("sample_metadata2.tsv",  sep='\t', header=T,  comment="")
head(metadata)

##   X.SampleID  number       plant  fertilizer      system
## 1  #q2:types numeric categorical categorical categorical
## 2   sample04       4  VignaMaize         Ctr         Mix
## 3   sample05       5  VignaMaize         Ctr         Mix
## 4   sample06       6  VignaMaize         Ctr         Mix
## 5   sample10      10 MucunaMaize         Ctr         Mix
## 6   sample11      11 MucunaMaize         Ctr         Mix

phyloseqオブジェクトを作成する

作成したtaxtable、otu、metadataを読み込む

physeq<-phyloseq(
  otu_table(otu$data, taxa_are_rows = T),
  phy_tree(tree$data),
  tax_table(as.data.frame(taxtable) %>% select(-Consensus) %>% column_to_rownames("Feature.ID") %>% as.matrix()), 
  sample_data(metadata %>% as.data.frame() %>% column_to_rownames("X.SampleID"))
  )

phyloseqオブジェクトの中身を確認する

physeq@tax_table[1:3,1:7] %>%kable(format ="markdown")

	Kingdom	Phylum	Class	Order	Family	Genus	Species
06aab2a98e551f6229bdf29ea062c46a	Bacteria	Acidobacteria	Holophagae	Subgroup 7	NA	NA	NA
66722b9fcf4300b9dd1725e8c6a9378e	Bacteria	Chloroflexi	Anaerolineae	Anaerolineales	Anaerolineaceae	uncultured	uncultured bacterium
7a5de3cf5a9c7c48a8329c204c7ede9b	Bacteria	Chloroflexi	Anaerolineae	Anaerolineales	NA	NA	NA

physeq@sam_data[1:3,1:3]

##          number      plant fertilizer
## sample04      4 VignaMaize        Ctr
## sample05      5 VignaMaize        Ctr
## sample06      6 VignaMaize        Ctr

physeq@otu_table[1:5,1:5] %>%kable(format ="markdown")

	sample05	sample10
06aab2a98e551f6229bdf29ea062c46a	0	0
66722b9fcf4300b9dd1725e8c6a9378e	0	0
7a5de3cf5a9c7c48a8329c204c7ede9b	7	6
708b575f46bb26cc0a46b32e815b21f7	0	0
6e5a27b92e5481892407f2e1b97ace0f	0	0

リード数を最小リードにそろえる（Rarefy）

リード数が最小値に調整されていないOTUテーブルからphyloseqオブジェクトを作成した場合、後からrarefy_even_depth関数を使ってrarefyできる。
rarefyされていないOTUテーブルを読み込む。

setwd("C:/Users/akari/Desktop/fastq2")
otu_unrarefied <- read_qza("table.qza")
head(otu$data[1:5,1:5])

##                                  sample04 sample05 sample06 sample10 sample11
## b5f5167bc021a2b0e18f93d765f303b7      459      657      856      826      978
## cee5e095946d1f90373e37e3275355d6      481      295      347      433      346
## a45ebb3a1adc92b2a15285afbefa78bf      381      362      396      366      320
## 59cfcf40b213772307e70ee970f231f0      500      287      362      324      336
## af6972e30898eb9a809d47a8e94809dc      311      482      308      331      304

同様にphyloseqオブジェクトを作成した。

physeq_unrarefied<-phyloseq(
  otu_table(otu_unrarefied$data, taxa_are_rows = T),
  phy_tree(tree$data),
  tax_table(as.data.frame(taxtable) %>% select(-Consensus) %>% column_to_rownames("Feature.ID") %>% as.matrix()), 
  sample_data(metadata %>% as.data.frame() %>% column_to_rownames("X.SampleID"))
  )

各サンプルのリード数を確認する。

colSums(otu_table(physeq_unrarefied))

## sample04 sample05 sample06 sample10 sample11 sample12 sample16 sample17 sample18 sample19 sample20 
##    54541    69549    57381    62342    60117    57082    56844    59873    48362    48639    50972 
## sample21 sample25 sample26 sample27 sample31 sample32 sample33 sample37 sample38 sample39 sample40 
##    51361    49143    54346    53944    62389    59152    53505    54991    62433    59392    52108 
## sample41 sample42 sample46 sample47 sample48 sample52 sample53 sample54 sample58 sample59 sample60 
##    62283    61249    51743    75021    53190    76523    51661    41095    54015    72819    66860 
## sample61 sample62 sample63 
##    58771    51131    41266

最小リード数の確認。

min(colSums(otu_table(physeq_unrarefied)))

## [1] 41095

以下のコマンドを使ってrarefyする。

physeq_rarefy <- rarefy_even_depth(physeq_unrarefied, sample.size = 41095)

最小リード数を確認する。

colSums(otu_table(physeq_rarefy))

## sample04 sample05 sample06 sample10 sample11 sample12 sample16 sample17 sample18 sample19 sample20 
##    41095    41095    41095    41095    41095    41095    41095    41095    41095    41095    41095 
## sample21 sample25 sample26 sample27 sample31 sample32 sample33 sample37 sample38 sample39 sample40 
##    41095    41095    41095    41095    41095    41095    41095    41095    41095    41095    41095 
## sample41 sample42 sample46 sample47 sample48 sample52 sample53 sample54 sample58 sample59 sample60 
##    41095    41095    41095    41095    41095    41095    41095    41095    41095    41095    41095 
## sample61 sample62 sample63 
##    41095    41095    41095

補足

rarefyの他にもscale_reads関数を使ってサンプルをnormalizationする方法もある。

scale_reads <- function(physeq,n){
  physeq.scale <- transform_sample_counts(physeq, function(x) {round((n*x/sum(x)))})
  #otu_table(physeq.scale) = round(otu_table(physeq.scale))
  physeq.scale = prune_taxa(taxa_sums(physeq.scale) > 0, physeq.scale)
  return(physeq.scale)
}

実行

physeq_scaled <- scale_reads(physeq_unrarefied, 41095)
colSums(otu_table(physeq_scaled))

## sample04 sample05 sample06 sample10 sample11 sample12 sample16 sample17 sample18 sample19 sample20 
##    41118    41057    41087    41115    41082    41090    41114    41107    41095    41101    41097 
## sample21 sample25 sample26 sample27 sample31 sample32 sample33 sample37 sample38 sample39 sample40 
##    41092    41115    41095    41100    41092    41100    41096    41073    41109    41116    41105 
## sample41 sample42 sample46 sample47 sample48 sample52 sample53 sample54 sample58 sample59 sample60 
##    41101    41053    41089    41123    41083    41094    41108    41095    41087    41091    41120 
## sample61 sample62 sample63 
##    41064    41098    41181

Phyloseqオブジェクトの作成は以上です。
分かりにくいところなどあればコメント欄に質問してください！

参考

https://rstudio-pubs-static.s3.amazonaws.com/330760_8bba830836324bf6b100d4e76f49e3d2.html

https://rdrr.io/bioc/phyloseq/man/rarefy_even_depth.html