深入学习ggtree: read.tree()是如何解析newick数据格式
最近在学习ggtree，为了更好的理解ggtree的，我去找了一些资料研究，Y叔在统计之都的文章《使用ggtree实现进化树的可视化和注释》提到
除了ggtree之外，我所了解到的其它画树软件在画树的时候都把树当成是线条的集合。很明显画出来的进化树就是在画一堆线条，但是线条表示的是父节点和子节点的关系，除此之外没有任何意义，而节点在进化树上代表物种，叶子节点是我们构建进化树的物种，内部节点是根据叶子节点推断的共同祖先。我们所有的进化分析、推断、实验都是针对节点，节点才是进化树上有意义的实体。这是ggtree设计的基础，ggtree只映射节点到坐标系统中，而线条在 geom_tree 图层中计算并画出来。这是与其它软件最根本的不同，也是ggtree能够简单地用图层加注释信息的基础。
那么问题来了，read.tree()是如何解析(B:6.0,(A:5.0,C:3.0,E:4.0):5.0,D:11.0);这种数据呢？于是我就去都源码了。
数据输入部分
if (!is.null(text)) {
if (!is.character(text))
stop("argument `text' must be of mode character")
tree <- text
} else {
tree <- scan(file = file, what = "", sep = "\n", quiet = TRUE,
skip = skip, comment.char = comment.char, ...)
}第一个条件语句，决定数据的读取方式。如果text不为null，则读取文本格式的数据，比如说
newick <- '(B,(A,C,E),D);'
read.tree(text=newick)如果text为null，则使用scan函数读取文件。根据scan函数的特点，当file=””，也就是说你没有指定数据的存放路径时候，他会让你从键盘上以parenthetic format读取数据，也就是说
> tree <- read.tree()1: (B,(A,C,E ),D);
2: (B, (A,C,E),D);
3: (B, (A,C,E),D) ;
> tree
3 phylogenetic trees于是就输入了3个进化发育树。scan的what=””表示读取的数据是字符型。sKip根据输入决定到底忽略多少行，sep=”\n”表示用换行符分割字符。如果没有这个参数，上述的结果就是7个字符串组成的向量。
if (identical(tree, character(0)))
{
warning("empty character string.")
return(NULL)
}第二个条件语句，利用identical判断scan输入的数据是否为空。这是因为如果scan()只要有一行输入为空，就结束输入。如果一行都没有输入，那么结果就是character(0)。
将输入拆分成tree名字和tree的字符串
tree <- gsub("[ \t]", "", tree)
tree <- unlist(strsplit(tree, NULL))后续，先用gsub把tree向量中的制表符和空格都去掉，我们刚才的输入里就掺杂了许多没有必要的空格。
然后用unlist(strsplit(tree,NULL))可以所有tree里面的字符串拆分成单个字符，如果有多行输入，会变成单行的向量.strsplit当第二个参数长度为0时会把字符串拆分当个单个字符串。如果输入时多个字符串的向量，返回的就是list，然后用unlist转换成一个向量。效果如下：
> unlist(strsplit(tree, NULL))
> tree
[1] "(" "B" "," "(" "A" "," "C" "," "E" ")" "," "D" ")" ";" "(" "B" "," "(" "A"
[20] "," "C" "," "E" ")" "," "D" ")" ";" "(" "B" "," "(" "A" "," "C" "," "E" ")"
[39] "," "D" ")" ";"也就是单个字符所组成的向量。
y <- which(tree == ";")
Ntree <- length(y)
x <- c(1, y[-Ntree] + 1)
if (is.na(y[1]))
return(NULL)后面就是用which找到”；“的位置赋值给y，用length判断scan输入了多少tree。然后把不同树起始位置传入给x。
c(1, y[-Ntree] + 1)
[1] 1 15 29后面的条件语句用于判断是否你是否真的输入了newick格式，毕竟你有可能输入错了。
STRING <- character(Ntree)
for (i in 1:Ntree) {
tmp <- paste(tree[x[i]:y[i]], sep = "", collapse = "")
STRING[i] <- gsub("\\[[^]]*\\]", "", tmp)
}character用于创建空字符串所组成的向量。接着用for循环，根据之前的位置信息，把数据读入到STRING中。R里面用paste组合多个字符。gsub("\\[[^]]*\\]", "", tmp)中的正则表达式直接看不好理解，可以分成三个部分
\\[： 先匹配到[[^]]*： 中间是任意多个除]（^])的字符\\]：最后是]结尾
所以功能就是删除形如[xxx]的内容。
注： 换做我们，可能会直接把STRING赋值为一个空字符串的向量，然后在循环中增加，如下：
mystring = c()for(i in 1:Ntree){
tmp <- paste(tree[x[i]:y[i]], sep="", collapse = "")
mystring <- c(gsub("\\[[^]]*\\]","",tmp ),mystring)
}由于每次循环都要重新在内存中开辟位置，因此会降低速度。而原代码一开始就定义长度，避免了动态创建，提高了计算机效率。
tmp <- unlist(lapply(STRING, unname))
tmpnames <- tmp[c(TRUE, FALSE)]
STRING <- tmp[c(FALSE, TRUE)]这部分首先用到了lapply对STRING循环套用函数unname。unname在前面定义了，功能就是把输入的结果拆分成两个个部分，如果是”ABC(A,(B,D),E)”，那么输出是(“ABC”,”(B,D),E)”)
unname <- function(treetext) {
# 输入字符的长度
nc <- nchar(treetext)
# 定义text的起始位置为1
tstart <- 1
# 找到(的字符串的起始位置
while (substr(treetext, tstart, tstart) != "(" && tstart <=
nc) tstart <- tstart + 1
# 返回一个向量，由”(“前的内容和"(”后的内容组成。
if (tstart > 1)
return(c(substr(treetext, 1, tstart - 1), substr(treetext,
tstart, nc)))
return(c("", treetext))
}
然后分别name和string赋值给tmpname,和STRING。
**注**：这里用到了循环补齐。
if (is.null(tree.names) && any(nzchar(tmpnames)))
tree.names <- tmpnames这个就是当函数输入中tree.names没有指定，且tmpnames也没找到，就把tree.name命名为tmpnames。这里的nzchar是快速判断向量元素是否为空的方式，any判断是否有一个为TRUE。
将tree的字符串拆分成分支和分支长度
colon <- grep(":", STRING)
if (!length(colon)) {
obj <- lapply(STRING, clado.build)
}else if (length(colon) == Ntree) {
obj <- lapply(STRING, tree.build)
}else {
obj <- vector("list", Ntree)
obj[colon] <- lapply(STRING[colon], tree.build)
nocolon <- (1:Ntree)[!1:Ntree %in% colon]
obj[nocolon] <- lapply(STRING[nocolon], clado.build)
}有一部分数据是(B:6.0,(A:5.0,C:3.0,E:4.0):5.0,D:11.0);。根据“：”的数目决定使用clado.build还是tree.build.这两个函数这次就不说，直接看下效果。
[[1]]
Phylogenetic tree with 5 tips and 2 internal nodes.
Tip labels:
[1] "A" "B" "A" "C" "NA"
Node labels:
[1] "D" "E"
Unrooted; no branch lengths.
[[2]]
Phylogenetic tree with 5 tips and 2 internal nodes.
Tip labels:
[1] "B" "A" "C" "E" "D"
Unrooted; no branch lengths.
[[3]]
Phylogenetic tree with 5 tips and 2 internal nodes.
Tip labels:
[1] "B" "A" "C" "E" "D"
Unrooted; no branch lengths.最后一部分大致有点难读懂，需要改天把clado.build和tree.build看懂才行。
总结
看到了如下函数的使用：character(), scan(), identical(), gsub(), lapply(), unlist(), nzchar(), grep().
学习了R语言里面如何使用正则表达式，知道[^]]就能用来匹配所有非]的字符，而不是[^\\]]。
了解到了一个小技巧，提前定义向量长度，提高效率。
阅读原文，可以跳转到Y叔的那篇文章