在R (ape库)中对基于简约的树进行根操作时出错

Question

我正在使用最大简约算法( ape库中的pratchet())在R中执行引导程序分析。当我对无根树(使用pratchet()函数生成)运行分析时，引导运行得很好。但是，当我想要在找到引导支持之前对每个引导树进行根操作时，我在随机地对100棵树进行根操作时遇到了错误。请注意，这发生在调用任何用于计算双分区或分支支持的代码之前。

如果我使用邻居连接算法(ape中的nj())，则在根或下游引导中根本没有问题，但显然在使用外组对基于简约的树进行根操作时(随机)会出现这种情况。我观察到的奇怪的事情是，如果我在对未根树进行根操作之前将它们写入文件(以防在根操作过程中发生错误)，然后想要对它们进行根操作，那么它可以很好地工作。

这是我用来分析的代码。

performBootstraping = function(charMatrix, bsIterations) {
    # charMatrix is a DNA alignment matrix
    # bsIteration are number of bootstrap iterations.

    library(phangorn)

    phySeq = phyDat(charMatrix)

    treeMPRooted = getRootedParsimonyTree(charMatrix)

    bValuesMP = boot.phylo(treeMPRooted, charMatrix, FUN=function(xx) {tt = getRootedParsimonyTree(xx); return(tt) }, 
    B = bsIterations, trees=T, rooted=T)

   # convert to percentage
   bValues = bValuesMP$BP/bsIterations * 100;

   plot(treeMPRooted, use.edge.length = F); 
   title('Max Parsimony tree with bootstrap percentage')
    nodelabels(bValues, frame = 'rect')

    # write the tree as newick 
    write.tree(treeMPRooted, paste0(outDir,'/rooted_MP.tree'))
    return(bValuesMP)
}

getRootedParsimonyTree = function(cMatrix) {
    phySeq = phyDat(cMatrix);
    treeMP = pratchet(phySeq)
    treeMPRooted = root(treeMP, outgroup='Germline', resolve.root=T)

    return(treeMPRooted)
}

下面是堆栈跟踪和错误

 Error in phy$edge[sndcol, 2] <- newNb[phy$edge[sndcol, 2]] <- n + 2:phy$Nnode : 
  number of items to replace is not a multiple of replacement length 
6 root(treeMP, outgroup = "Germline", resolve.root = T) at libaryFunctions.R#105
5 getRootedParsimonyTree(xx) at libaryFunctions.R#32
4 FUN(x[, boot.samp]) 
3 boot.phylo(treeMPRooted, t.vaf, FUN = function(xx) {
    tt = getRootedParsimonyTree(xx)
    return(tt)
}, B = bstrapCount, trees = T, rooted = T) at libaryFunctions.R#32
2 performBootstraping(vaf, outDir, i, bsIterations) at runAllSitesBootstrapForAllPatients.R#15
1 runAllSitesBootstrapForAllPatients(ccfFileDir = ccfDir, outDirPref = outDir) 
In addition: Warning message:
In newNb[phy$edge[sndcol, 2]] <- n + 2:phy$Nnode :
  number of items to replace is not a multiple of replacement length

Answer 1

我对“root”函数感到失望，所以编辑了ape源代码，创建了一个在不抛出错误的情况下将树作为根的函数。它作为RootTree()函数包含在'TreeTools‘包中，该函数还实现了一些用于设置根目录的其他选项。

使用以下命令将其安装到R：

install.packages('TreeTools')
RootTree(tree, outgroup)

恐怕它不支持边缘长度。

Answer 2

看看phangorn中的函数bootstrap.phyDat和plotBS的例子，应该会简化你的函数。此外，您不需要对单个引导程序树进行根操作，只需创建要绘制的树即可。这应该可以解决您的问题。

致敬克劳斯