必要なパッケージ

require(DECIPHER) # bioconductorからインストールする必要がある
require(seqinr)
require(ape)

はじめに

　ここでは、微生物叢のメタゲノムデータの解析（https://www.biorxiv.org/content/10.1101/2021.10.04.462986v2.full.pdf ）に用いられた、ダイズの根圏微生物の16S rRNA（リボゾーマルRNA）のデータの一部を用いて解析を行う。なお、rRNAは、ウィルスを除く全生物に存在して、進化速度が比較的遅いことから、種分類などによく用いられている。メタゲノム解析では、微生物叢の多様性と系統間の関係を把握するために16S rRNAを用いられることが多い（https://bifidus-fund.jp/keyword/kw009.shtml ）。

配列データの読み込み

最初に、FASTA形式で保存された16S rRNAの配列データを読み込む。

seq <- readDNAStringSet("soy_rmbiom_d.fas", format = "fasta")
seq

## DNAStringSet object of length 100:
##       width seq                                             names               
##   [1]   253 TACGAAGGGGGCTAGCGTTGCT...GAAAGCGTGGGGAGCAAACAGG s001
##   [2]   253 TACGAAGGGGGCTAGCGTTGCT...GAAAGCGTGGGGAGCAAACAGG s002
##   [3]   253 TACGTAGGGTGCAAGCGTTGTC...GAAAGCGTGGGGAGCGAACAGG s003
##   [4]   253 TACGTAGGTGGCAAGCGTTATC...GAAAGCGTGGGGAGCAAACAGG s004
##   [5]   253 TACGTAGGGCGCAAGCGTTGTC...GAAAGCGTGGGGAGCGAACAGG s005
##   ...   ... ...
##  [96]   253 TACGTAGGGTGCGAGCGTTAAT...GAAAGCGTGGGGAGCAAACAGG s096
##  [97]   253 TACGAAGGGGGCTAGCGTTGCT...GAAAGCGTGGGGAGCAAACAGG s097
##  [98]   253 TACGTAGGGCGCAAGCGTTGTC...GAAAGCGTGGGGAGCGAACAGG s098
##  [99]   253 TACGAAGGGGGCTAGCGTTGTT...GAAAGCGTGGGGAGCAAACAGG s099
## [100]   253 TACGTAGGGTGCGAGCGTTGTC...GAAAGCGTGGGGAGCGAACAGG s100

全部で100の16S rRNAの配列がある。

DECIPHERパッケージのAlignSeqs関数を用いて100の配列間でマルチプルアライメントを行う。

seq.aligned <- AlignSeqs(seq)

## Determining distance matrix based on shared 8-mers:
## ================================================================================
## 
## Time difference of 0.07 secs
## 
## Clustering into groups by similarity:
## ================================================================================
## 
## Time difference of 0.02 secs
## 
## Aligning Sequences:
## ================================================================================
## 
## Time difference of 0.28 secs
## 
## Iteration 1 of 2:
## 
## Determining distance matrix based on alignment:
## ================================================================================
## 
## Time difference of 0 secs
## 
## Reclustering into groups by similarity:
## ================================================================================
## 
## Time difference of 0.02 secs
## 
## Realigning Sequences:
## ================================================================================
## 
## Time difference of 0.18 secs
## 
## Iteration 2 of 2:
## 
## Determining distance matrix based on alignment:
## ================================================================================
## 
## Time difference of 0 secs
## 
## Reclustering into groups by similarity:
## ================================================================================
## 
## Time difference of 0.02 secs
## 
## Realigning Sequences:
## ================================================================================
## 
## Time difference of 0.03 secs

アライメントした結果を確認する。

seq.aligned

## DNAStringSet object of length 100:
##       width seq                                             names               
##   [1]   257 TACGAAGGGGGCTAGCGTTGCT...GAAAGCGTGGGGAGCAAACAGG s001
##   [2]   257 TACGAAGGGGGCTAGCGTTGCT...GAAAGCGTGGGGAGCAAACAGG s002
##   [3]   257 TACGTAGGGTGCAAGCGTTGTC...GAAAGCGTGGGGAGCGAACAGG s003
##   [4]   257 TACGTAGGTGGCAAGCGTTATC...GAAAGCGTGGGGAGCAAACAGG s004
##   [5]   257 TACGTAGGGCGCAAGCGTTGTC...GAAAGCGTGGGGAGCGAACAGG s005
##   ...   ... ...
##  [96]   257 TACGTAGGGTGCGAGCGTTAAT...GAAAGCGTGGGGAGCAAACAGG s096
##  [97]   257 TACGAAGGGGGCTAGCGTTGCT...GAAAGCGTGGGGAGCAAACAGG s097
##  [98]   257 TACGTAGGGCGCAAGCGTTGTC...GAAAGCGTGGGGAGCGAACAGG s098
##  [99]   257 TACGAAGGGGGCTAGCGTTGTT...GAAAGCGTGGGGAGCAAACAGG s099
## [100]   257 TACGTAGGGTGCGAGCGTTGTC...GAAAGCGTGGGGAGCGAACAGG s100

上と変わっていないようにも見えるが、配列の長さが253から257になっている。長さの増加は、アライメントの際にギャップが入れられたことによる。

より完全に確認するには、DECIPHERパッケージのBrowseSeqs関数を用いる。

BrowseSeqs(seq.aligned, highlight=1)

アライメントした結果をFASTA形式のファイルに出力する。

seqs <- as.list(as.character(seq.aligned))
names <- names(seq.aligned)
write.fasta(seqs, names, "soy_rmbiom_d_aligned.fas",
            open = "w", as.string = TRUE)

近隣結合法

保存しておいたアライメントされた配列データを改めて読み込む。

rmbiom <- read.dna("soy_rmbiom_d_aligned.fas", format = "fasta")
rmbiom

## 100 DNA sequences in binary format stored in a matrix.
## 
## All sequences of same length: 257 
## 
## Labels:
## s001
## s002
## s003
## s004
## s005
## s006
## ...
## 
## Base composition:
##     a     c     g     t 
## 0.247 0.194 0.362 0.198 
## (Total: 25.7 kb)

塩基の組成も確認できる。

配列間の距離を計算する。まずは、Jukes and Cantor（1969）を使って距離を計算する。

dst.jc69 <- dist.dna(rmbiom, model = "JC69")

一部を確認してみる。

as.matrix(dst.jc69)[1:5, 1:5]

##             s001        s002       s003      s004       s005
## s001 0.000000000 0.003978789 0.28645096 0.2922876 0.29228756
## s002 0.003978789 0.000000000 0.28645096 0.2864510 0.29228756
## s003 0.286450958 0.286450958 0.00000000 0.2523542 0.04919796
## s004 0.292287563 0.286450958 0.25235418 0.0000000 0.24132457
## s005 0.292287563 0.292287563 0.04919796 0.2413246 0.00000000

次に、Kimura’s 2-parametersモデルで距離を計算する。

dst.k80 <- dist.dna(rmbiom, model = "K80")

一部を確認してみる。

as.matrix(dst.k80)[1:5, 1:5]

##             s001        s002       s003      s004       s005
## s001 0.000000000 0.003984085 0.28889981 0.2962456 0.29511107
## s002 0.003984085 0.000000000 0.28889981 0.2899562 0.29511107
## s003 0.288899805 0.288899805 0.00000000 0.2527193 0.04941107
## s004 0.296245634 0.289956243 0.25271933 0.0000000 0.24230115
## s005 0.295111065 0.295111065 0.04941107 0.2423012 0.00000000

体勢は変わらないが、少し異なる値になっていることがわかる。

JC69の距離をもとに、近隣結合法で系統樹を構築する。

tre.nj.jc69 <- nj(dst.jc69)

無根系統樹を描く。

plot(tre.nj.jc69, type = "u", lab4ut = "axial")

次に、K80の距離をもとに、近隣結合法で系統樹を構築する。

tre.nj.k80 <- nj(dst.k80)

無根系統樹を描く。

plot(tre.nj.k80, type = "u", lab4ut = "axial")

このデータの場合は、JC69もK80も大きく変わらないことが分かる。

最後にブートストラップ法を適用してみる。

まずは、K80で近隣結合法による系統樹を構築する自作関数を作る。

myfun <- function(x) nj(dist.dna(x, model = "K80"))

自作関数を使って系統樹を構築してみる。

tree <- myfun(rmbiom)

系統樹を描く。

plot(tree, type = "u", lab4ut = "axial")

問題なく構築できている。

では、作成しておいた関数と系統樹を使ってブートストラップ法を適用する。

1000回のブートストラップを行う。

nboot <- 1000
bstrees <- boot.phylo(tree, rmbiom, myfun, B = nboot, trees = T)$trees

## 
Running bootstraps:       100 / 1000
Running bootstraps:       200 / 1000
Running bootstraps:       300 / 1000
Running bootstraps:       400 / 1000
Running bootstraps:       500 / 1000
Running bootstraps:       600 / 1000
Running bootstraps:       700 / 1000
Running bootstraps:       800 / 1000
Running bootstraps:       900 / 1000
Running bootstraps:       1000 / 1000
## Calculating bootstrap values... done.

数え上げをする。

boot <- prop.clades(tree, bstrees)
boot

##  [1] 1000  570  253  385  391  161  231  230  173  334  895  966  576  252  134
## [16]  308  209  336  647  867  720  837  646  998  490  818  600  704  990  984
## [31] 1000  996  614  365  990  743  557  357  208  897  755  579  153  462  332
## [46]  644  313  307  514  915  501  603  999  959  812  994 1000  461  924  711
## [61]  777  485  315  254  491  781  523  563  695  365  997  694  776  603  519
## [76]  554  695  998 1000  874  984 1000  983  995 1000  603  594  508  368  411
## [91] 1000 1000  696  983  799  504  664  965

10で割って100分率（％）にする。小数点以下は丸める。

boot <- round(boot / (nboot / 100))

ブートストラップ値がついた系統樹を描く。

plot(tree, type = "u", lab4ut = "axial")
drawSupportOnEdges(boot)

各ブートストラップ値は、その位置における二分割が生じた確率を表す。

高い値は、その位置における分割が確からしいことを表す。

バイオメトリックス第9回　系統樹作成

岩田洋佳 hiroiwata@g.ecc.u-tokyo.ac.jp

2022/12/12

必要なパッケージ

はじめに

配列データの読み込み

近隣結合法

バイオメトリックス第9回 系統樹作成

岩田洋佳 hiroiwata@g.ecc.u-tokyo.ac.jp

2022/12/12

必要なパッケージ

はじめに

配列データの読み込み

近隣結合法

バイオメトリックス第9回　系統樹作成