用phylotools建立supermatrix

phylotools是用来建立DNA supermatrix矩阵的R程序包，是为了将比对好的基因文件（Fasta文件或Phylip文件）合并成超级矩阵（supermatrix）而编写的。

Supermatrix与Relaxed Phylip格式

supermatrix，就是将不同基因的片段分别比对后，再按照物种名合并形成的大型数据集，以便建立进化树用。

supermatrix一般为Relaxed Phylip格式。这是因为能够建立大型进化树的软件，如RAxML、Q-Tree等，都能识别这种格式。

Relaxed Phylip格式是Phylip格式的扩展，即第一行的两个数字分别声明序列数和位点数，两个数字中间有一个空格。第二行以序列名开始，序列名称和序列数据之间，有一个空格。传统的Phylip文件，序列名称只允许10个字符，如果超过10个字符，Phylip等软件在建树时会将名称截断，只取前10个字符。Relaxed Phylip格式则没有这种限制。

建立supermatrix的几个常用片段

为什么建立进化树要用多个DNA片段呢？这主要是因为DNA条形码的单个基因一般都太短，难以全面反映物种之间的关系，所以建立更为可信的进化树，一般需要更长的序列，即所谓多基因策略。

不同的序列，由于编码功能的差别，有些突变的速度快，有些速度慢。有些基因片段就较为稳定，在不同科属之间差异较小，而有些则变异很快，在同一个种不同个体之间都有很大的区别。在DNA条形码研究中，常用比较稳定的rbcLa基因确定植物物种的大致类群，一般可确定到科的水平；用matK等变异稍快的类群鉴定到属甚至到种的水平；用trnH-psbA等变异速率非常快的非编码基因确定更精确的种或种下的水平等。

上述三个DNA条形码序列，都是属于叶绿体基因组。当然，除上述基因外，还有很多别的基因，如ITS等也可用来建树。

DNA条形码序列一般是对PCR产物用传统的桑格测序法获得的，对样品和反应条件的要求较高，有一定的失败率。理想情况下， 一份样品的几个条形码序列都应该是完整的。但是由于各种各样的原因，有些种的某些序列总是不能正常扩增和测序，造成数据缺失。

supermatrix的基本内容

在建树开始之前，每个基因的序列需要单独比对，然后再拼接在一起。当然，这时候需要保证同一个种的不同序列要接在一起，但是位点的顺序不能错位。与此同时，若某一个种的某一个基因的序列缺失，则其对应的全部位点都需要用?补齐。

类似软件

用比对好的序列文件建立supermatrix，除了phylotools，还有以下软件：

  1. Biopython (https://biopython.org/wiki/Concatenate_nexus)
  2. Utensils (https://github.com/ballesterus/Utensils)
  3. Geneious (https://assets.geneious.com/manual/11.1/GeneiousManualsu53.html)
  4. catfasta2phyml (https://github.com/nylander/catfasta2phyml)
  5. SEDA (http://www.sing-group.org/seda/ , https://www.biostars.org/p/332853/)
  6. SuperMatrix: (http://coleoguy.blogspot.com/2015/04/r-concatenate-alignments-into.html)
  7. concat (https://figshare.com/articles/concatenation_R/3839538)
  8. AMAS (https://www.ncbi.nlm.nih.gov/pmc/articl::es/PMC4734057/)

其中Genious可以通过图形界面实现。

Phylotools的特点

  • phylotools全部用R写成，操作简单。
  • phylotools假设同一个种，在不同的fasta或者phylip文件中有相同的名称。
  • phylotools的supermatrix可以识别比对好的fasta文件或者phylip文件，生成的supermatrix以Relaxed Phylip格式保存。

phylotools的安装

phylotools的最新版保存在github上，而并未上传到CRAN。要安装phylotools，可使用如下命令：

1
devtools::install_github("helixcn/phylotools")

若未安装devtools，可使用 install.packages("devtools")安装

phylotools程序包函数列表

  • clean.fasta.name 用于清除fasta序列名中的特殊字符，以便建树的时候能够正常输入和显示
  • dat2fasta 将读取为data.frame的DNA序列转换为fasta文件
  • dat2phylip 将读取为data.frame的DNA序列转换为phylip格式
  • get.fasta.name 获取fasta文件中所有序列的名称
  • get.phylip.name 获取phylip文件中所有序列的名称
  • read.fasta 读取fasta文件
  • read.phylip 读取phylip文件
  • rename.fasta 重命名fasta文件
  • rm.sequence.fasta 从fasta文件中删除某一条特定的序列
  • split_dat 根据序列所在的分组，将数据分成独立的fasta文件，例如，可用于trnH-psbA序列按照目(orders)分组
  • sub.taxa.label 替换进化树中的物种名
  • supermat 用已经比对好的fasta文件或者phylip文件建立supermatrix超级矩阵以及RAxML分隔文件。

用phylip文件生成Relaxed Phylip格式的supermatrix

假设现有三个Phylip文件，分别编码不同的基因，建立者三个基因的supermatrix，并保存为Relaxed Phylip格式，同时生成RAxML Partition File，记录每个基因的起讫位点。

相应代码如下：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
cat("6 22",
"seq_1 --TTACAAATTGACTTATTATA",
"seq_2 GATTACAAATTGACTTATTATA",
"seq_3 GATTACAAATTGACTTATTATA",
"seq_5 GATTACAAATTGACTTATTATA",
"seq_8 GATTACAAATTGACTTATTATA",
"seq_10---TACAAATTGAATTATTATA",file = "matk.phy", sep = "\n")

cat("5 15",
"seq_1  GATTACAAATTGACT",
"seq_3  GATTACAAATTGACT",
"seq_4  GATTACAAATTGACT",
"seq_5  GATTACAAATTGACT",
"seq_8  GATTACAAATTGACT",
file = "rbcla.phy", sep = "\n")

cat("5 50",
"seq_2 GTCTTATAAGAAAGAATAAGAAAG--AAATACAAA-------AAAAAAGA",
"seq_3 GTCTTATAAGAAAGAAATAGAAAAGTAAAAAAAAA-------AAAAAAAG",
"seq_5 GACATAAGACATAAAATAGAATACTCAATCAGAAACCAACCCATAAAAAC",
"seq_8 ATTCCAAAATAAAATACAAAAAGAAAAAACTAGAAAGTTTTTTTTCTTTG",
"seq_9 ATTCTTTGTTCTTTTTTTTCTTTAATCTTTAAATAAACCTTTTTTTTTTA",
file = "trn1.phy", sep = "\n")

supermat(infiles = c("matk.phy", "rbcla.phy", "trn1.phy"))

用fasta文件生成supermatrix

要使用比对好的fasta文件生成supermatrix， 则命令改为对应的fasta文件名即可：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
cat(">seq_1",  "--TTACAAATTGACTTATTATA",
">seq_2",  "GATTACAAATTGACTTATTATA",
">seq_3",  "GATTACAAATTGACTTATTATA",
">seq_5",  "GATTACAAATTGACTTATTATA",
">seq_8",  "GATTACAAATTGACTTATTATA",
">seq_10", "---TACAAATTGAATTATTATA",
file = "matk.fasta", sep = "\n")

cat(
">seq_1", "GATTACAAATTGACT",
">seq_3", "GATTACAAATTGACT",
">seq_4", "GATTACAAATTGACT",
">seq_5", "GATTACAAATTGACT",
">seq_8", "GATTACAAATTGACT",
file = "rbcla.fasta", sep = "\n")

cat(
">seq_2", "GTCTTATAAGAAAGAATAAGAAAG--AAATACAAA-------AAAAAAGA",
">seq_3", "GTCTTATAAGAAAGAAATAGAAAAGTAAAAAAAAA-------AAAAAAAG",
">seq_5", "GACATAAGACATAAAATAGAATACTCAATCAGAAACCAACCCATAAAAAC",
">seq_8", "ATTCCAAAATAAAATACAAAAAGAAAAAACTAGAAAGTTTTTTTTCTTTG",
">seq_9", "ATTCTTTGTTCTTTTTTTTCTTTAATCTTTAAATAAACCTTTTTTTTTTA",
file = "trn1.fasta", sep = "\n")

supermat(infiles = c("matk.fasta", "rbcla.fasta", "trn1.fasta"))

参考文献

  • Kress, W. J., Erickson, D. L., Jones, F. A., Swenson, N. G., Perez, R., Sanjur, O., & Bermingham, E. (2009). Plant DNA barcodes and a community phylogeny of a tropical forest dynamics plot in Panama. Proceedings of the National Academy of Sciences, 106(44), 18621-18626.
  • 裴男才, 张金龙, 米湘成, 葛学军 (2011) 植物DNA条形码促进系统发育群落生态学发展, 生物多样性, 19(3):284–294
  • Roquet, C., Thuiller, W., & Lavergne, S. (2013). Building megaphylogenies for macroecology: taking up the challenge. Ecography, 36(1), 13-26.
  • https://cme.h-its.org/exelixis/web/software/raxml/index.html