search

liamrevell / phytools

GNU General Public License v3.0
207 stars 56 forks source link

Error in solve.default(C) : Lapack routine dgesv: system is exactly singular: U[393,393]=0 #89

Open vanessayang0927 opened 3 years ago

vanessayang0927 commented 3 years ago

Hey phytools! I'm trying to use the phylosig function to computes phylogenetic signal. It's going well using the lambda method. But it's wrong when using the K method. I have no idea what the issue is, but it seems a very common question. And there are some information about the error:

R version 4.0.5 phytools version 0.7-70

file files.zip

my code: library(phytools) my_tree <- read.tree("otus1.16s.pick725.Neighbor-Joinning.nwk") abun <- read.table("test.txt", header = TRUE, sep = "\t", row.names = 1, check.names = FALSE) C <- setNames(abun$C,rownames(abun)) phylosig(my_tree,C,method = "K",test = TRUE) phylosig(my_tree,C,method = "lambda",test = TRUE)

error: Error in solve.default(C) : Lapack routine dgesv: system is exactly singular: U[393,393]=0

I'm trying phylosig(my_tree,C*100,method = "K",test = TRUE) but still get the same error

liamrevell commented 3 years ago

Hi Vanessa. The problem is that some of the tips of your tree are separated by patristic distance of zero. (For instance, a pair of sister species with zero length terminal branches.) You can see this by running:

D<-cophenetic(my_tree) min(as.dist(D))

(The function as.dist converts the symmetric distance matrix to a lower/upper diagonal object of class "dist" so the diagonal is ignored.)

This makes the among species covariance matrix uninvertible ("exactly singular").

There's no single solution to this, but perhaps it would be OK to just prune one or the other tip that is affected in this way?

-- Liam

vanessayang0927 commented 3 years ago

Hi Liam. Thanks for your reply. I tried to use as.dist and get the same error. There are some information about the error:

R version 4.0.5 my code: D <- cophenetic(my_tree) tree <- ape::njs(D) phylosig(tree,C,method = "K",test = TRUE)

error: Error in solve.default(C) : Lapack routine dgesv: system is exactly singular: U[301,301] = 0

It seems to be correlated with negative or zero branch lengths produced by neighbor-joining. I'm trying to make negative or zero branch lengths a really small positive number(1e-10) and it's going well. my_tree$edge.length[which(my_tree$edge.length <=0)] <- 1e-10 phylosig(my_tree,C,method = "K",test = TRUE) phylosig(my_tree,C,method = "lambda",test = TRUE)