添加链接
link之家
链接快照平台
  • 输入网页链接,自动生成快照
  • 标签化管理网页链接
Collectives™ on Stack Overflow

Find centralized, trusted content and collaborate around the technologies you use most.

Learn more about Collectives

Teams

Q&A for work

Connect and share knowledge within a single location that is structured and easy to search.

Learn more about Teams

I am having a hard time understanding why my phylo tree (imported as a phyloseq object from QIIME2) appears to have no branch lengths when used in phylosig() . I am trying to compute the phylogenetic signal of my 16S dataset compared to a single continuous metadata variable. All example datasets are included at the bottom of this question.

    glacialpath <- metadata[,c(7)]
    phylosig(phylo, glacialpath, method = "K", test = TRUE)
Error in vcv.phylo(tree) : the tree has no branch lengths
    In addition: Warning message:
    In drop.tip(tree, setdiff(tree$tip.label, names(x))) :
      drop all tips of the tree: returning NULL

My tree phylo was separated from a complete phyloseq object phylo <- phy_tree(physeq = physeq) and has confirmed branch lengths:

Phylogenetic tree with 1290 tips and 1289 internal nodes.
Tip labels:
  e54924083c02bd088c69537d02406eb8, 3112899fc7a2adb7f74a081a82c7cde4, db5d745b02355b6fed513c4953b62326, 2faf2046aa9ab2f6598df79cd52e9c7b, bec8db81cc8ec453136c82ede8327a9f, 601e47b8adcbd21d159f74421710e1b5, ...
Node labels:
  0.942, 0.563, 0.711, 0.999, 0.000, 0.528, ...
Rooted; includes branch lengths.

A dput() of phylo is too large to include, so hopefully what I have provided is enough. Any help is appreciated.

> dput (metadata)
structure(list(sample = structure(1:48, .Label = c("sample-10", 
"sample-11", "sample-12", "sample-14", "sample-16", "sample-18", 
"sample-19", "sample-20", "sample-21", "sample-22", "sample-23", 
"sample-24", "sample-25", "sample-26", "sample-27", "sample-28", 
"sample-29", "sample-30", "sample-31", "sample-32", "sample-33", 
"sample-34", "sample-35", "sample-36", "sample-37", "sample-40", 
"sample-41", "sample-43", "sample-44", "sample-45", "sample-46", 
"sample-50", "sample-55", "sample-56", "sample-57", "sample-58", 
"sample-59", "sample-61", "sample-62", "sample-63", "sample-64", 
"sample-65", "sample-67", "sample-68", "sample-69", "sample-70", 
"sample-71", "sample-8"), class = "factor"), pond = structure(c(1L, 
1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 
3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 
5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 1L), .Label = c("Lilly", 
"RHM1", "RHM2", "SS", "TS"), class = "factor"), meh = structure(c(16L, 
17L, 19L, 4L, 1L, 16L, 18L, 20L, 3L, 5L, 7L, 10L, 11L, 22L, 1L, 
2L, 14L, 16L, 18L, 20L, 3L, 5L, 7L, 10L, 1L, 15L, 16L, 18L, 19L, 
20L, 21L, 8L, 1L, 2L, 14L, 15L, 16L, 18L, 19L, 20L, 21L, 3L, 
5L, 6L, 9L, 12L, 13L, 1L), .Label = c("0-1", "1-2.", "10-11.", 
"11-12.", "12-13.", "13-14.5", "14-15", "14-16.", "14.5-17.5", 
"16-17", "17-19", "17.5-18.5", "18.5-20", "2-3.", "3-4.", "4-5.", 
"5-6.", "6-7.", "7-8.", "8-9.", "9-10.", "9-11."), class = "factor"), 
    depth = c(4, 5, 7, 11, 0, 4, 6, 8, 10, 12, 14, 16, 17, 9, 
    0, 1, 2, 4, 6, 8, 10, 12, 14, 16, 0, 3, 4, 6, 7, 8, 9, 14, 
    0, 1, 2, 3, 4, 6, 7, 8, 9, 10, 12, 13, 14, 19, 18, 0), om = structure(c(29L, 
    13L, 9L, 9L, 27L, 27L, 26L, 26L, 23L, 23L, 12L, 12L, 20L, 
    34L, 24L, 24L, 22L, 22L, 25L, 25L, 21L, 21L, 11L, 11L, 8L, 
    7L, 5L, 3L, 3L, 2L, 1L, 4L, 33L, 30L, 31L, 32L, 28L, 18L, 
    19L, 17L, 15L, 17L, 16L, 14L, 18L, 19L, 10L, 6L), .Label = c("1.2", 
    "1.3", "1.4", "1.5", "1.7", "19.9", "2.1", "2.6", "2.8", 
    "2.9", "27.9", "29.8", "3.2", "3.3", "3.4", "3.5", "3.6", 
    "3.7", "3.8", "30.2", "30.6", "32.1", "32.5", "32.6", "32.9", 
    "35.7", "36.2", "4.2", "4.3", "5.4", "5.8", "6", "6.4", "na"
    ), class = "factor"), water = structure(c(19L, 16L, 12L, 
    13L, 28L, 28L, 35L, 35L, 30L, 30L, 31L, 31L, 28L, 36L, 27L, 
    27L, 33L, 33L, 34L, 34L, 31L, 31L, 29L, 29L, 20L, 21L, 15L, 
    5L, 3L, 2L, 1L, 4L, 26L, 23L, 24L, 25L, 22L, 13L, 14L, 7L, 
    8L, 9L, 11L, 10L, 17L, 18L, 6L, 32L), .Label = c("16.1", 
    "17", "17.5", "18.8", "21.6", "21.8", "22.3", "22.4", "22.6", 
    "22.7", "23.1", "23.6", "24.4", "24.5", "24.6", "25.5", "26.9", 
    "28", "29.1", "31.8", "32.8", "34.9", "40.5", "42.9", "43.2", 
    "50.9", "59", "60.1", "61.4", "61.5", "62.5", "62.9", "63.2", 
    "63.8", "66.4", "na"), class = "factor"), glacialpath = c(19.4, 
    19.4, 19.4, 19.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 
    6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 6.4, 
    18.7, 18.7, 18.7, 18.7, 18.7, 18.7, 18.7, 18.7, 16.8, 16.8, 
    16.8, 16.8, 16.8, 16.8, 16.8, 16.8, 16.8, 16.8, 16.8, 16.8, 
    16.8, 16.8, 16.8, 19.4), elevation = c(8.2, 8.2, 8.2, 8.2, 
    18.3, 18.3, 18.3, 18.3, 18.3, 18.3, 18.3, 18.3, 18.3, 18.3, 
    18.3, 18.3, 18.3, 18.3, 18.3, 18.3, 18.3, 18.3, 18.3, 18.3, 
    13.4, 13.4, 13.4, 13.4, 13.4, 13.4, 13.4, 13.4, 5.5, 5.5, 
    5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 5.5, 
    5.5, 8.2), area = c(20536L, 20536L, 20536L, 20536L, 3571L, 
    3571L, 3571L, 3571L, 3571L, 3571L, 3571L, 3571L, 3571L, 3571L, 
    3571L, 3571L, 3571L, 3571L, 3571L, 3571L, 3571L, 3571L, 3571L, 
    3571L, 10369L, 10369L, 10369L, 10369L, 10369L, 10369L, 10369L, 
    10369L, 57003L, 57003L, 57003L, 57003L, 57003L, 57003L, 57003L, 
    57003L, 57003L, 57003L, 57003L, 57003L, 57003L, 57003L, 57003L, 
    20536L), tectonics = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 
    1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L), .Label = c("outwash", 
    "uplift"), class = "factor"), lat = c(60.52, 60.52, 60.52, 
    60.52, 60.478, 60.478, 60.478, 60.478, 60.478, 60.478, 60.478, 
    60.478, 60.478, 60.478, 60.478, 60.478, 60.478, 60.478, 60.478, 
    60.478, 60.478, 60.478, 60.478, 60.478, 60.491, 60.491, 60.491, 
    60.491, 60.491, 60.491, 60.491, 60.491, 60.425, 60.425, 60.425, 
    60.425, 60.425, 60.425, 60.425, 60.425, 60.425, 60.425, 60.425, 
    60.425, 60.425, 60.425, 60.425, 60.52), long = c(145.601, 
    145.601, 145.601, 145.601, 145.382, 145.382, 145.382, 145.382, 
    145.382, 145.382, 145.382, 145.382, 145.382, 145.382, 145.382, 
    145.382, 145.382, 145.382, 145.382, 145.382, 145.382, 145.382, 
    145.382, 145.382, 145.547, 145.547, 145.547, 145.547, 145.547, 
    145.547, 145.547, 145.547, 145.473, 145.473, 145.473, 145.473, 
    145.473, 145.473, 145.473, 145.473, 145.473, 145.473, 145.473, 
    145.473, 145.473, 145.473, 145.473, 145.601), x.depth = c(1L, 
    1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 
    2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 
    3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 
    4L, 1L), y.depth = c(2L, 3L, 4L, 5L, 1L, 2L, 3L, 4L, 5L, 
    6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 
    19L, 20L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 1L, 2L, 3L, 4L, 
    5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 1L)), class = "data.frame", row.names = c(NA, 
-48L))

The problem is not the lack of branch lengths in the tree. What the error message is saying is that the tree does not have the tips that correspond to the names of values in the trait variable. glacialpath must be a named vector and the names must be present in the phylogeny.

In general, a good practice is to check for possible problems in these three areas.

  • phylo is not a correct phylo format.
  • phylo should contain tips with names corresponding to those in the metadata.
  • Related to #2, glacialpath lacks names.
  • phylo is not a correct phylo format

    phytools::phylosig requires that the tree is in correct phylo format. Try to explore the tree with str(phylo), and see whether all values in phylo$edge.length are numeric and there are no missing values.

    phylo should contain tips with names corresponding to those in the metadata

    What are the samples the phylogenetic signal should be calculated for? Assuming that the column sample in the metadata contains names, reduce the tree to the size of the available data:

    tr2 <- ape::keep.tip(phy = phylo, tip = metadata$sample)
    plot(tr2) # to check that the tree was trimmed correctly
    

    While phytools::phylosig can trim the tree during processing, preparing a subset of the tree beforehand is preferable, because the user can check that the trimmed tree plots as expected, including branch lengths. In the given example, the tree tips are not names present in the trait vector names.

    glacialpath lacks names

    The glacialpath is not a named vector. To set the names of a vector, use:

    glacialpath <- setNames(object = metadata[, 7], nm = metadata$sample)
    

    Lastly, a good practice is to name variables differently to names of functions or, in this case, classes.

    Thanks for contributing an answer to Stack Overflow!

    • Please be sure to answer the question. Provide details and share your research!

    But avoid

    • Asking for help, clarification, or responding to other answers.
    • Making statements based on opinion; back them up with references or personal experience.

    To learn more, see our tips on writing great answers.