Go to Qiita Advent Calendar 2024 Top
LoginSignup
3
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 3 years have passed since last update.

@akari5(木村 有歌理)

phyloseqオブジェクトの作成

Last updated at Posted at 2020-06-16

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")
sample04 sample05 sample06 sample10 sample11
06aab2a98e551f6229bdf29ea062c46a 0 0 0 0 0
66722b9fcf4300b9dd1725e8c6a9378e 0 0 0 0 0
7a5de3cf5a9c7c48a8329c204c7ede9b 0 7 0 6 0
708b575f46bb26cc0a46b32e815b21f7 0 0 0 0 0
6e5a27b92e5481892407f2e1b97ace0f 0 0 0 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

3
2
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
3
2

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?