phyloseqオブジェクトの作成

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

``` r
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は“；”で区切られている。

``` r
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が対応するようなテーブルを生成する  

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

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

<table>
<colgroup>
<col style="width: 18%" />
<col style="width: 4%" />
<col style="width: 8%" />
<col style="width: 16%" />
<col style="width: 11%" />
<col style="width: 11%" />
<col style="width: 11%" />
<col style="width: 11%" />
<col style="width: 5%" />
</colgroup>
<thead>
<tr class="header">
<th style="text-align: left;">Feature.ID</th>
<th style="text-align: left;">Kingdom</th>
<th style="text-align: left;">Phylum</th>
<th style="text-align: left;">Class</th>
<th style="text-align: left;">Order</th>
<th style="text-align: left;">Family</th>
<th style="text-align: left;">Genus</th>
<th style="text-align: left;">Species</th>
<th style="text-align: right;">Consensus</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td style="text-align: left;">b5f5167bc021a2b0e18f93d765f303b7</td>
<td style="text-align: left;">Bacteria</td>
<td style="text-align: left;">Proteobacteria</td>
<td style="text-align: left;">Alphaproteobacteria</td>
<td style="text-align: left;">Sphingomonadales</td>
<td style="text-align: left;">Sphingomonadaceae</td>
<td style="text-align: left;">Sphingomonas</td>
<td style="text-align: left;">uncultured bacterium</td>
<td style="text-align: right;">0.8</td>
</tr>
<tr class="even">
<td style="text-align: left;">cee5e095946d1f90373e37e3275355d6</td>
<td style="text-align: left;">Archaea</td>
<td style="text-align: left;">Thaumarchaeota</td>
<td style="text-align: left;">Soil Crenarchaeotic Group(SCG)</td>
<td style="text-align: left;">NA</td>
<td style="text-align: left;">NA</td>
<td style="text-align: left;">NA</td>
<td style="text-align: left;">NA</td>
<td style="text-align: right;">0.9</td>
</tr>
<tr class="odd">
<td style="text-align: left;">a45ebb3a1adc92b2a15285afbefa78bf</td>
<td style="text-align: left;">Bacteria</td>
<td style="text-align: left;">Chloroflexi</td>
<td style="text-align: left;">KD4-96</td>
<td style="text-align: left;">uncultured bacterium</td>
<td style="text-align: left;">NA</td>
<td style="text-align: left;">NA</td>
<td style="text-align: left;">NA</td>
<td style="text-align: right;">0.9</td>
</tr>
<tr class="even">
<td style="text-align: left;">59cfcf40b213772307e70ee970f231f0</td>
<td style="text-align: left;">Bacteria</td>
<td style="text-align: left;">Verrucomicrobia</td>
<td style="text-align: left;">Spartobacteria</td>
<td style="text-align: left;">Chthoniobacterales</td>
<td style="text-align: left;">DA101 soil group</td>
<td style="text-align: left;">uncultured bacterium</td>
<td style="text-align: left;">NA</td>
<td style="text-align: right;">0.8</td>
</tr>
<tr class="odd">
<td style="text-align: left;">af6972e30898eb9a809d47a8e94809dc</td>
<td style="text-align: left;">Bacteria</td>
<td style="text-align: left;">Acidobacteria</td>
<td style="text-align: left;">Acidobacteria</td>
<td style="text-align: left;">Subgroup 6</td>
<td style="text-align: left;">uncultured bacterium</td>
<td style="text-align: left;">NA</td>
<td style="text-align: left;">NA</td>
<td style="text-align: right;">0.6</td>
</tr>
<tr class="even">
<td style="text-align: left;">a742c9230d6676357d2c79ac0a42c785</td>
<td style="text-align: left;">Bacteria</td>
<td style="text-align: left;">Proteobacteria</td>
<td style="text-align: left;">Betaproteobacteria</td>
<td style="text-align: left;">Nitrosomonadales</td>
<td style="text-align: left;">Nitrosomonadaceae</td>
<td style="text-align: left;">uncultured</td>
<td style="text-align: left;">uncultured bacterium</td>
<td style="text-align: right;">0.7</td>
</tr>
</tbody>
</table>

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

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

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

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

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

``` r
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を読み込む

``` r
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を読み込む

``` r
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オブジェクトの中身を確認する

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

<table>
<colgroup>
<col style="width: 25%" />
<col style="width: 6%" />
<col style="width: 10%" />
<col style="width: 9%" />
<col style="width: 11%" />
<col style="width: 12%" />
<col style="width: 8%" />
<col style="width: 15%" />
</colgroup>
<thead>
<tr class="header">
<th style="text-align: left;"></th>
<th style="text-align: left;">Kingdom</th>
<th style="text-align: left;">Phylum</th>
<th style="text-align: left;">Class</th>
<th style="text-align: left;">Order</th>
<th style="text-align: left;">Family</th>
<th style="text-align: left;">Genus</th>
<th style="text-align: left;">Species</th>
</tr>
</thead>
<tbody>
<tr class="odd">
<td style="text-align: left;">06aab2a98e551f6229bdf29ea062c46a</td>
<td style="text-align: left;">Bacteria</td>
<td style="text-align: left;">Acidobacteria</td>
<td style="text-align: left;">Holophagae</td>
<td style="text-align: left;">Subgroup 7</td>
<td style="text-align: left;">NA</td>
<td style="text-align: left;">NA</td>
<td style="text-align: left;">NA</td>
</tr>
<tr class="even">
<td style="text-align: left;">66722b9fcf4300b9dd1725e8c6a9378e</td>
<td style="text-align: left;">Bacteria</td>
<td style="text-align: left;">Chloroflexi</td>
<td style="text-align: left;">Anaerolineae</td>
<td style="text-align: left;">Anaerolineales</td>
<td style="text-align: left;">Anaerolineaceae</td>
<td style="text-align: left;">uncultured</td>
<td style="text-align: left;">uncultured bacterium</td>
</tr>
<tr class="odd">
<td style="text-align: left;">7a5de3cf5a9c7c48a8329c204c7ede9b</td>
<td style="text-align: left;">Bacteria</td>
<td style="text-align: left;">Chloroflexi</td>
<td style="text-align: left;">Anaerolineae</td>
<td style="text-align: left;">Anaerolineales</td>
<td style="text-align: left;">NA</td>
<td style="text-align: left;">NA</td>
<td style="text-align: left;">NA</td>
</tr>
</tbody>
</table>

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

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

``` r
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テーブルを読み込む。

``` r
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オブジェクトを作成した。

``` r
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"))
  )
```

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

``` r
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

最小リード数の確認。

``` r
min(colSums(otu_table(physeq_unrarefied)))
```

    ## [1] 41095

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

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

最小リード数を確認する。
``` r
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する方法もある。

``` r
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)
}
```

実行

``` r
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オブジェクトの作成は以上です。
分かりにくいところなどあればコメント欄に質問してください！
### 参考

<a href="https://rstudio-pubs-static.s3.amazonaws.com/330760_8bba830836324bf6b100d4e76f49e3d2.html" class="uri">https://rstudio-pubs-static.s3.amazonaws.com/330760_8bba830836324bf6b100d4e76f49e3d2.html</a>

<a href="https://rdrr.io/bioc/phyloseq/man/rarefy_even_depth.html" class="uri">https://rdrr.io/bioc/phyloseq/man/rarefy_even_depth.html</a>