search

KlausVigo / phangorn

Phylogenetic analysis in R
http://klausvigo.github.io/phangorn/
203 stars 38 forks source link

Write ancestral reconstruction to file #93

Closed surh closed 4 years ago

surh commented 4 years ago

Is there any way to write the most ancestrally reconstructed sequences to a fasta file?

If I use wite.phyDat on the output of any ancestral recontruction I just get a sequences with all 'a'.

KlausVigo commented 4 years ago

Hi @surh ,

that should not happen. Which version of phangorn/R are you using? Here is some code which works on my machine. I start with simulating a tree and an alignment.

library(phangorn)
tree <- rcoal(6)
plot(tree)
nodelabels()
dat <- simSeq(tree, l=15)

It is important to use the argument return = "phyDat"

anc <- ancestral.pars(tree, dat, return = "phyDat")
as.character(anc)
write.phyDat(anc, file = "test.fas", format = "fasta")

If you have problems can you send me your script and the dataset/tree if possible. Kind regards, Klaus

surh commented 4 years ago

It was my bad. I was using the output of ancestral reconstruction with return="prob". Still, it would be nice if write.phyDat threw and error when this happens, instead of silently writing an incorrect file.

Sur

KlausVigo commented 4 years ago

No worries! I will add a check to this. Klaus

DiracZhu1998 commented 1 year ago

Hi Sur, Thank you for giving us such a wonderful tool! I would like to ask why the sequence after running ancestral.pars changed from original sequence? Appreciate your help! Code: # Please change NJtree.txt as NJtree.rds for usage NJtree <- readRDS("NJtree.rds") binary <- read.delim("test.txt") phyDat.binary <- phyDat(t(binary), type = "USER", levels = c("a","b")) anc.acctran <- ancestral.pars(NJtree, phyDat.binary, "ACCTRAN") as.character(anc.acctran) # sequence is different from binary sequence

test.txt NJtree.txt

DiracZhu1998 commented 1 year ago

solved already. it's prob matrix.

KlausVigo commented 1 year ago

Hi @DiracZhu1998, glad you figured it out. I am working on how to export ancestral states in a better way and also on some more tools to present them. In the development version I added an experimental function write.ancestral as a more or less tidy data.frame. The problem is that the output is not an alignment as one has probabilities for each state.