Howdy Biostars,
I have DNA sequence that I have successfully done a bunch of tree building for in BEAST2 to get a time calibrated tree (taxa=8) that answers all the questions I have about topology and branch length. However, I need to know the ancestral states for each bifurcation in the tree. There doesn't seem to be a ready-made way of doing this in BEAST2 (though it seems like I could manually define each possible clade and completely rerun BEAST2, but that gets computationally expensive fast). I have run through the phangorn tutorial for doing ASR and indeed have used my data to do it. However, there does not seem to be a way to actually get the DNA sequence back. Instead, I get a phyDat object. I can tell the ancestral states are there (9-14 in the names, below) but I have only 77 different site patterns from my 179 input characters (all of which are phylogenetically informative--this is a concatenated alignment). I'd like to get back to characters (that is, base calls) from the site patterns, which seem to be defined outside of the object somewhere, but I cannot find out where.
The contents of anc.ml are otherwise as I expect (6 ancestral states, 8 extant taxa) but I simply cannot figure out how to access the sequence. Any help understanding this data structure that leads to me being able to pull out the sequence information is appreciated. If you have another programmatic option to do this (must support GTR) I am also open to that as a solution. I can do it in FastML online, but I need to be able to have this as part of my own pipeline.
Below is the trace of calling the object itself and the contained attributes in R
> anc.ml 14 sequences with 179 character and 77 different site patterns. The states are a c g t
> attributesanc.ml)
$names
[1] "Met1" "Met2" "Met3" "Met4" "Met5" "Met6" "Normal" "Prim" "9" "10" "11" "12" "13" "14"
$class
[1] "phyDat"
$weight [1] 1 1 1 1 7 1 3 3 3 8 1 2 2 7 1 11 2 1 1 1 3 1 1 1 3 1 3 1 6 4 6 3 3 6 1 1 2 5 2 7 6 3 1 2 2 5 1 1 2 1 1 1 2 3 2 1 2 1 1 1 1 1 1 2 3 1 1 1 1 2 1 1 1 1 [75] 1 1 2
$nr [1] 77
$nc [1] 4
$index [1] 1 2 3 4 5 6 7 8 9 10 8 11 5 12 13 10 14 15 16 17 18 19 17 20 21 22 23 7 24 25 26 27 28 14 29 30 16 31 32 21 33 16 32 34 8 35 36 37 38 39 40 31 41 42 30 41 41 29 43 14 44 39 37 5 30 14 34 29 45 14 10 46 47 30 [75] 48 49 40 27 34 46 9 50 51 9 52 5 53 54 54 10 55 42 27 10 56 5 33 57 16 58 59 60 40 57 21 16 31 61 5 29 31 40 62 14 13 46 38 40 63 64 65 25 16 66 65 25 34 41 67 65 53 68 44 34 69 16 40 16 70 10 64 14 31 40 71 38 34 29 [149] 38 72 5 41 54 49 73 10 70 74 42 31 41 16 75 76 7 32 77 77 45 12 16 33 29 16 10 55 38 46 46
$levels [1] "a" "c" "g" "t"
$allLevels [1] "a" "c" "g" "t" "u" "m" "r" "w" "s" "y" "k" "v" "h" "d" "b" "n" "?" "-"
$type [1] "DNA"
$contrast a c g t [1,] 1 0 0 0 [2,] 0 1 0 0 [3,] 0 0 1 0 [4,] 0 0 0 1 [5,] 0 0 0 1 [6,] 1 1 0 0 [7,] 1 0 1 0 [8,] 1 0 0 1 [9,] 0 1 1 0 [10,] 0 1 0 1 [11,] 0 0 1 1 [12,] 1 1 1 0 [13,] 1 1 0 1 [14,] 1 0 1 1 [15,] 0 1 1 1 [16,] 1 1 1 1 [17,] 1 1 1 1 [18,] 1 1 1 1
$call ancestral.pml(object = fit2, type = "ml")
