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.