A few days ago a phytools user wrote:
“I apologize if this has been discussed elsewhere, but I am curious if you
have an easy way to extract and compare homologous edges of a tree that all
contain the same taxa but have different edge lengths and are written
differently in there newick format? I would use grep to extract these normally
but I have a case where I dated a large topology and the tree output is written
differently than the tree I input in newick format. I want to check
correlation of dates to the actual branch lengths but they are not in the same
order in the newick format.”
I'm sure there are a number of ways to do this in R; however, one way is using
the phytools function matchNodes
.
library(phytools)
trees
## 9 phylogenetic trees
## first, the trees are all topologically identical:
densityTree(trees,type="cladogram",nodes="intermediate")
![plot of chunk unnamed-chunk-1](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAfgAAAH4CAMAAACR9g9NAAAA5FBMVEUAAAAAADoAAGYAOpAAZmYAZpAAZrYfH/8jI/8nJ/8sLP8xMf83N/86AAA6ADo6AGY6Ojo6OmY6OpA6ZmY6ZrY6kJA6kLY6kNs+Pv9GRv9PT/9ZWf9kZP9mAABmADpmAGZmOgBmOpBmZgBmZjpmZmZmZpBmZrZmkJBmkNtmtv9wcP9+fv+Ojv+QOgCQOjqQOmaQZgCQZjqQkDqQtpCQ29uQ2/+goP+0tP+2ZgC2Zjq2Zma2kDq2tma225C2/7a2///Kyv/bkDrbtmbb25Db/7bb///j4///tmb/25D//7b//9v///9cw7kTAAAACXBIWXMAAAsSAAALEgHS3X78AAAgAElEQVR4nO2dC4PbNpLnnWTPO4mdyW527tItMZmbu32ZpJO9hwnP7t0FvTdBt43v/32W9QLBhyQCfIgS8Z9Mu1utpij9gEIBKFS9skm71Ktr30DSdZTA71QJ/E6VwO9UCfxOlcDvVAn8TpXA71QJ/E6VwO9UCfxOlcDvVAn8TpXA71QJ/E6VwO9UCfxOlcDvVAn8TpXA71QJ/E6VwO9UCfxOlcDvVAn8TpXA71QJ/D3ozx+C/ySBvwN9+uOvwX+TwN+BXr7/LfhvEvg70PPr8L9J4G9O+tiIHnn6IfwqCfxtSB0HRb+M8O0S+A1LV/kQ7Lz1pLq3x/h2CfzWpKtiiHZx6vlPr+3zV+G+XQK/Belh2nmpL//tpx9ffRnR4RP4q+kkbbXKyyfwqwpo9wfufC3anhL4xaW12gptTwn8MkLafVMOtEcM3CsogZ9PujbkVdnr3HlebIW2pwR+os7RvqYpv6QEPkJ6mHZeG/Ky2l7nHlQCP1abpf30qlbwan0Cf04AWyHtvAP72rQ9fXwHyzjvAv8qge/J0S42TLsRrdQH78wm8Cg9TBtgF0h7a7hJpZUojOBYjD2D1zRsn6S9uc7N8u7VymZ88EbN7sAL7arwB+5867TN4A5tZWmIx026MO0CvD5Du9wwba2zIdyt53z+GaIwPv0UGotxv+B92kWfdqk2OnBrnQ/gzk7dKw7un9/vfDqnxUkj2rkPe+O0qwHceW4u/uUzzOK/2GVc/QjamzTlpr7pPOvhznN9Gfd03Sh4PUw7FyeNaG8QN9Au+7jrFroK7ka3BP6Waeuqhzurb1qbdXE32jp4ho1OGk64CXe+edoGaGPnzrq0q+vhbrRJ8F3axW3Rrrq4M3A2NoG70WbAa9+QN7QZNjppOOHeHm5Hu4U7A9rYubeEu9F1wQ/Rzvu0N9i5Dd83du6sTbvcLu5GVwA/TDv3Dfk2aZuGto/b0b4a7ueIyPqVwOvxtDeGWxuF7kZZZC3afN/Yua/cuz/99DfbWsAZoF3cBm1W3cXhtiHq5uHh6zdv3ta9++g692aM+dMPH0PDMOYH72CTSy60y9uhTSKzriqYfHM/P2TH+r2Um+FdC89PfvrjrxHnpGcB3+naZ2hvb+Dui5hXVZHX/yvwa1UprWBaXr+Toro2etVsOcGPT+9iMiNEg69nKqwyR9o59IdjWT5mnksutLePm0TM6xYLzCkgHqw6vs/6H6NVbfhhkrb6nWk9cD4DPtZPf/otJhdKLPjKgS+zvKJ/igq69rG4NdokWmlzzEtmzogBPrYAo/BdrmO6oBUWRXdhv2i99EfYngt366PBu1avazNo+J/6c6k9n8hLXk+0lg6TNGJe6YGpGWA3vBBbLbnDi2Nm0eWd58PnM2g//udgtz4WvHYdGu2gsWIOoZ9EXvMawt0yCJ+umcNK2yBzfio8juhre1bVs/k5XRZ2kcCnyHzaef0q5z7Pz+/RoQ936+PBW37T0A80dvZ6BETwt9LlYSnGMS+ry5M0+rXBNwnD2UT2Rje4Pd7oDqvzvKdrAnhDXZu6u4FYh7oroAe0/S7vM0effdyqm+GWwfNV7PmB7I2bAZU+b579wMXW+fQiwSNcNOqEue7s9QhkoJ1i/5/1HucVR8/XE/R6Uj6eOYvJG5nzgdkn9mdNhTfh5bluJqt/hHst3k5TwEtnx58rBQOg0kh+m+ANdvMCVl7zCOZyFbL3ln1CnuEPsG9o179G3Njc3OJGdQ3eTpPAW5rcwnfwFtDprXs+jPXz3eEs8phz9Hz8jooxxvkC2Ifxaghft6S4e9M6kHTv6rKBWEPTwFsBD6zJu4fARrMh/87gBgsyp26upm+ieeRpWlA3JvDzYZbnVi7rFyuFd0HdewO8nSaCr4EzeMMzOfha+04bIG885tlczJtrl5XXsXM24hn8B/8j3KXjvRngojjwTYNXtF5jcLTHDxZ9uyvP6WhFBplj1JOakTm/Qt3NK1Z9bVXmx8cDbN4dH4+PuBC02vLex6jA+kjwzmQZcOqsYfAGjbzW1/Psjax7kZEV5gu9kL87VdBmVN37C+ryhWxCLhug8fmXXyOOzk0Hj7aeFrWoq9M36w/zPMEC5rwisxTz4Vf22wD0d/ivaX68rrnA7eCabfj23HTwOLwTeF7Pxt6wJnn64GkBrOB1r3WYd2/EbwPk09PObiEzuHkbAGzLIvOPK52da8CzWad+bmkfQ8iv0dl4hxD6GDLHxZANRMe0mgDmUiH61AKqSoaACXdqqhwyWcelwIkD709mYIRX1M3xZ+z2OAIsO8zLBwsGFVdfhfmCrxkl07QBVcp0PicvoMJNoZgGYMC7hEAM2phbqcd3wOPqDe/ZaBrxLe/bxFz98strx7zEmVRebZN5W6ZpAtRSeQYIa8eVbAuOehNwHbAiEAZGKXBWcu66PZ7m7oYfYPIwyZmbhXxwisJ72GpeceEzTtIEFL8HWEOmVlDS7PBcC8C/hcUiauro24UnvZoKnvZpaDPWuPdEXxXO72KuP/iSbrwsq8pnPtMLXEeEWHELyFAcpVi/u4EGQM+HdWJN4PF8/KuVomx98PgjbMjy7jz3d1rImYW885Iq6A/4CRW4DHfbzNtiGwCLvgU1gFxagERzc3+SXQEzbdoSAx5vgb7jDg+7cw58Q54C1eJvrvGMwQKSb0zG8J6Yd8Rx3bCplPEYAEdsafWfNgLRH5g4X50FPGCuxNRbRx/9u7jht0GuKNRlF8xbksj+CvzXAx3iwf9g46GoCjVtvhwHvuPb4U0qws3PkODEcPJe0ALKY35bLtw8QuuOW351+88Ph8cjnevIzTT/Zhp4f/LG/rw8gyOSbQB54zOHTS8KbS/3ytyTOPq1lYfwXlOZXNtJUQ9TwXPvphUczxQweTOuz7eRt5hf4/DCJsXxH2BbrzTG94d4y3be24Yi8rhRe558a3XbMS95OhN+d/esK3v1bfBI2E3iG8RMXs6fDF7HESc/Bn24IjE/p948vtbLt69erbFy14CnId5aIQv8/FtkT61Hvr2PiczLxHykOit31j59+av9/D6U/ETw3OHFxBvfC3MenrYN+W435wB3Zp7G81Hy1+rty3e/uq8Bmg7e0t4s/a61VCeTOgq/9JHLHBVjFZD5hsIQb0Judy6qhriNAd9M103H0tu2f8dPxZh7sE3czY3bZSnoHHViHifaj+fDc8GaCN5N2+SXHfAalhnqno1nzUxiPrOOdl3wjW/nlui83xnx8GneQXM0g8vtQJuYa72btdeFRTEYwQZ/Gng25p7rhg6ef6JE6Qo2nR7LMi8Oifn8eobQ6qfgEvLTe7yR4BsrizmYNICfo+in4yEv1DGDNBrBL5h0XrAjHzyNDwffGeLJezP+L7GPN9s40P9NBYGv1NnTqL4JTQDvdXgfvOXx3Vu9xZ+UHCpMCzRb0HTwtrH0tgu7/RNFu+vrpf5MajQRPI3qffC8KdtBTzvsCf0WNBm8PQH+LHqT0F9boeBdaBWtwvLCXeu33nMH0Bs+NjzxGEnSRMWDRwedWoBp/bb17IR+o5oCvmfpe+D76JtoOi0hJYn9NEUksrVx4HmIVxcsvXuwh56CSJRJ6OdQ3PbcJPCXLL17eBi9VrLWn9DHKyKdKSgQvOfbSXKRy+C76BsH30ef2IcL9uPp1GSwosFjeAU+NAL8IHotc7uEfoyqXhZrylcfkbIcNAE8d/hx4C+iT+wbYeKmwQrife6Rvl08eAn182ld4GZMO/S+OTrRXGpv6A0ekiuyfnHhrsoT5yeeVgnEGADfCq+8+Ofn0d+1o4eVjTD34SXMmLhp5McQ6dsFgpeeiUM8owoBT7HBbfTKGXz3jNtnjwFHRPk85gyjiydk7Yn07aLBSxYKEwb+DtHzkcaG8inOWU6QN5GcaQL4hlLn15cvMoAeAraUt6+/hc+mJ6xDh5CJ8mnMGR8C2wrmviLAyxDPj4z27fzLjEN/vQ8MctUWmI5YGBdYfy47HnqMs5JyFyu3Gnm12w5QJHijqsYRb/127IXa6OWU7AbQ4wGlosgfHo/Hw8PD4+Obb96A/vqvf/fNN28fHw+YpIAag77dwNFI8JpLIhljo8CPQr8ye0Nn0mi8pkTUXGUSrfkj9vVeh99Fj/eG+EosfSz4LaH3kWPyvIo50skAt5W02zHega8tvTwS7Nu1rmd8soxet5OYLYnetHp5QZkYBLket3G8B6++AV+JpZ8EPgT9zB8VXFGpZrJdNMiNRzz2Vdebx/MKTvA6Thx4L7nZNPC2+xHTph0Z/O6zwq9+6gUROdCGzNKCXEaeacjPvPIiK3e8ghO8jhMF3gwO8VHgx6Kfzt4I8lIqQrFLJucDlkN+7q4mr9XT7lzwHl0I+GZlTerHzgH+FHoMwze9p0VfX9E4TNnTSpc32l4L+RlF7M4F79HFgNeNpZ8F/DB6Tej7zwu9LI+2JRc3xz6uvZ3gTTE/rVP78ZivPjybbTh4Y1yqG+ov7V9HqodeiZc3gH7UXpDRWkqBeeVp/FzLt4L8tNz5+PA9uijwSs/l27Wu3iLBUdgqEL2R2p8Nch7Itf/nN4/cKdK3iwGvtduZmxP8CfS2O9TbYfaSYoeRU5U/Sn5rete/C+Qk9u3+EPp3AeCbIX7Q0k8HH4fez3BdcfUfRm7uGjmJvLrwEOsI8Go58LwU6FF1C6e9yYxxxVwVdXKapSFy3/G4X+aTFA4eq8i6B+YGPwJ9U9sH5+Ru7a3CxLcJ+TiFg1cC3iwDnqF30OOSupdbp2ZeSVGnkjp5Qh6iCPBSUc7t0bZ/P4Pa6A2tdeICDGXMwwyoboFdD47lCfl5jQfvLP2JIX5G8LxRignTKC8ilKLJwFeXGIiyUn5Utk7IAxUOXp+w9LOClxUcyIxHyy4qPx7z/PD2m7cHnKi1kigZV713zlu4a4WCNy7aqGvpZwbPZrsZ1asCywtmeSEqye/DJytjdtvpn2Lqj4WC13r5IX7gpdtNAEvPFa4FgFFoFhN31wA+fvVbRMXBQPBcfKb5eRXw3mv4LcDVm2w3gO72y53zf8aklsvtzokj1Rri1wbvvZjXArBQY9E0AM8E3HkDgF2MjwufnWNL74FvvtrOt2vKawHUALCMkTQAirO4L/66tS0L+3JQm2SxJMZs6d2e7NK+XbB6JqCQAo0FTw1uuQHofCgwo7D200+4IRvc78PAe5WFNgfeyXTWdIvclejkHPn6RvircigiK289h/bjF9yWxVm7NlKOojfEbwd8I1d1nJb6oFJ34beALTYAo8p8AHd+Mvkz1CJCzz5MY8FLqBJszrmftw7eqcm/QXv10gKoHI7S11740wprSPcseTXmbp4WTVvOHd5g7ht7a+CdGhNQkQmgFWAqoLBuAzBMu4Mbcvqv8VEGgYcvWjWW/vbAOzUmoOItgEwaQGcQmPll61En7+HOi2pixdBgBYDHT8EN8fKo94TblORgcg0g5xZAjuAMDQCnGUUXd8a0r/S5TQJ/ux1+SNIAcJcfWwDWQSzafuDIixmh3cYNraraROrukeCF+236dmEyfgvIqAXgugCewTjXAAzPH+jPGkMOtNUmcDcKAG9piOeVUH7Uf8a9yZiOCaj7a4kxXo0J8Gnn3lHZTGhvC3ej0B4vvp086j/jbkVVVRTsBwnc7HA8PjZpEg7OkFOz2MyywEmNBs9vZZfgncgElPnh8Ph4ePvw8Obh8eEAqTOK0g0Da9/T5/dLHqFCvy6BZ4EBqGr+EONblRrDxK7WxRdNYuwcGoPnmczOwYMMrQKpivMaXG0wXzSJsevwhn07frT1lJiXv2lhx6dD9lwZ/RofwqJJjKVgeALfkaK1vxXRm9xf8qsf+IhD/EJHqIyW7TgAP2Tpdwqeur2iVR7eC5r7gzC6GNyMJ/DLJjHmutHs2wnjPQ/xnnCyBzHA0usn5zyEIyT5CdhZ5+KLJjE2WvZl2beT+/PvNerl70IGV3rgpG4Vi95AqpZiGHaWqzNWJNK3Gwdetyx9At+TQfaw04sndfUYiw+sFW0M9VBneTl2yHjGIf6LYHs/DrwkQGgN8Qm8J5zuQsQ/LOVwr+9+JgZZV7S420WdwVCxYo31keBNCzw9mny7lgi9opSoYvG1jALdbTpiXVCytWvc7hjwxq1NJPDnRLQxZi47HLP6f4cD/OehztAH3MIhvzHgG0uPRQYT+JPCw5uqyB4fvv76d9/Ugn/ePhwwa8N0b39OjQDfpJg2WnvpjhL4vvCzgmw8eYZ5z49+ZYNb6/HNnlMCHyhjbniMx1NTdG9QJsw6z7Plgq7oj96sbsyr11YmpVqXyis35t+l2dLwdRuaaR7Pi/XBOzVjerwrMaa13+OTqZ9P8St3C+7Hm3GDfMyrJw0qZK1+yf14fwp3knzq8kvp/O7ckvvxpiFPVt8xTl1+dQ3vx38ZbO7H7c75xl6lLr8lLbsf3zL2ifyWtGxRYTw75cgbn7z/rKg7SJqk8ITlpIADFUJe+8N86vJXVnjCctLoQ5NeqN1pYx93D0lXUNhpWf7+FPnU5W9HITlwvGFeDw/zCfzNKCCzpfWHeY986vK3qLCUpg7sKQcvgb8VhRUq8LbmhpduU5e/FQWVHztl7G3q8jen0LpzbuP9xDCfuvzaWqHEKFLVHvmhYT6BX18xCayDwCPVE8O8e0oiv7ailusDwfvG3iRjvw1FhWKEge8Z+0Hy4XeRNF6qsx8fGYoRCJ6MvSNvhob51OVnV348IfztU0yNimDwSF5+0O7gvD0Rl5MUKX0KdtF9ZlwoRij4trEfHuZTl4+UOgX7bPB6XChGMPi2sT9BPuJGdipzCnY++qRCXJhtBPj2nK41zPtPSTolXZ2AXax4LCUcfG9Op1OXvyijTsDO1oTtKwJ8z9j3yacub+Gk4Yl4+LzcwIGzKPBt8mZgmN8peDgNN3wMqtgCbF8x4NvDPByt6w3zO+ryWqnh4055cfbM25UVBR7x+oFYOzP2kOFqEDaceN8wbF+R4Im8/DRIfvKtbUuQ32AQdl4Ua55rn0tx4DtzOkj60yV/B10eklkMHl2GShWbgf0UtR0fC75j7NvkW//clsxJ2HXH3g7srp6CC03Gg++R7zp2N9PlMfP8EGyoNLJqepJYPYcflo0Gz+QbZ37A2MdeenlxSqI+7OxmYKNK/Prpx8UjcHxdNPYb6/JYCa7KB2BDCaHbgS0q4d4hp/Tn9zHH5yaAP01+K12e6ocNwMaOvW6WqVlVlrQAjBXFP8bEYUwAL8beI2/aw/xVujx2bDWQUgwrwZULFBJYV5AnuYD3Uyl64DncsQNNAG9b4IeG+ZU+Yi4M2IfNdWJvHraoLJk6lELhx15+v2RGjGF1jb3tDvMLdnnq2APJArFjF9tIGzqjFEEvcqh2rVTzi8iSNJPADxj7zjA/62dvTsGmjn1/sFlQ6QqhF7rqUJ+gSeDF2Lte3h/mJ13eWfEB2NixN5IQeDkR9Jq6gtTI8FHMtck3EXyLPBp7M9XYc8fuw852ApulBLrS8pHMRt1OBd8z9rZn7EdeprHiXdjFjmCLFA7ptbiw1dzU7WTw0q9Pkj9FSytdqbL+7/B4zPPj4fh4aHfs3cFmSUcvteJK5UR97jiOyeAvDfMnyEGbht58/OZ3b968ffO2/u9wIJf1TuZeEWqgM3PGvkRpi6ngLw7zJ/qsNtDla/zQwXMsx865+6GMh9a3UIJ9Rumjqm0cMafS1MBc5/ObeNFk8GzsrW/sWzVIh+l55ak1pOyuu3pZ5lLHgYp4cCGXPbQArGDCAR1o4/GTySFt/aU/pS3Zj6EbdDOA7xr7zjA/TE0K3ZiKipbKVbBSq7SADFtAXlDB5jtuAjrLqrziz4xY12+1rIoRb/fz+3cxG/LTwfeNvW0P86fAWwbv/7W7JDaBEptA5oxAzjW7tbmvJqBr6gosPPxA+cHxE6hGvcvnr36L2JCfAbwz9ieG+cGbPw++eRrX7qynswNGQPzeGd7DNaWr+t1YLkuK7wd3vMDkX1gIwUW8z+//KWK9fhbw7MT5xt43AEP3LmfuNIG/7MlTCU9uAplrAugI1MbB3G4TqNs1eHCAmhoyf9WtnHLdPypxtQN/+xxeWXYe8EPG3ic/xENzyUpdYZnigCmcZwQKagFkA9gImBtrAga4s1tktBJ7aaQUUPedaKhWL2tcuB//8t3/jIjEmAV839j75AeNPTuAJhy896roCdAGRqsFlL4R2HYTgELE2n0cUN+NdrowN7z258O4fFfiO62RF5U8/vL7DzEHpWcCL8Z+wMEbNPbGgVfx4L2Xxxg6bALCH9b/6gegBejNNgEYvfADkI8D2mvjuJB/jyt3WJ6+4Fmu8t4LBtw9vQ5+6XnA98hbvxLxgLF3vt0s4P37gBZQ+UYgx8gbbAJ6Y00AuStGTF1c+dNW2o9TuCVL6ztsyJpLvHwLZj4i3HIm8GLsvV7eHua7TyfwZm7wzQvA1kY96GPwArUAKO7Hn57eQhMgO4XcCbyhquN1h64bLy1oQD+vfX5ot25Jb55Xnw38gLH3hvnO7RoH3oilW0wUXYtNgKcD0BawBVTXbAHE3ZVyhBLt4NerMqctC7TtRN9vqXO9/FzgB4y9Pd3lebK6CngnmB9BC3AlfnnLF5tAtXITMMKdXq0mDdMSbKEFuG5o4/Mcdy3qO1Oz39ls4BtjT990hvn+yhwv168H3sm1AFkVwn3gnFvAGlbAdLjb+tVr3NjDIQUSNMnK2hKmqHjHavbi8zOCl7mH58z75NvPdX3+CuCb25AWUPHacJ57RmA5VwAviFM3WfOCrq6wv4Mvl5fs1ZXc4+2We/yFYd6/a7MR8M39wKfLLQCtQKcFzNoE6CJYkt3taHBNNxjnK/blK1iU9OhvdYy3YuxPDfMt8A37LYB3Qs9ajAC6Ajk1gYqHgRmagM/dyAMy9PEHxcFW5OHhjvVmvXqQGRzm2W/1Z2wOvCXw88/mpsr4RgBnAQW3gEqMgLcqEHD3/OwOd/5E6Mfmycy/OjOPlzR3T6Hr9XOC94b5YQfPPW3z4J0MNQEZBlC0JIRNoFJhTcDn7gWnD3NnSdDx8ModZasPP08zK/gmpLox9o68D94Na5sHLzJ+C0AjUDYtoJIWcKEJyMNgtN3HYJyva3tOsK9Ta/XQ5a+4ckcaMPbeMO8dtiLw9nbAOxluAqppAe0moE41AdPjzg+M4i7q7859/vlDxIHZecGfMPZumOfnGHKkbxK8U68F0JEXaQGVGAFpApqn4gadtBb3vlG8qNZ+fET5+LnBe8Z+gLxu7ABzv2HwTqbdBHgSiA5hSU2AZugUMwbrxxSM4FlCbwYcJDpGF5W3fHbwYuwHh3ntnuKDn2+OclUZ0zUCZdMEqA3Q3jE2AbcLZ9zyXeynsEJNmhEaMPYNWW7qFoHfG3gnrwVQE8gLGQhgcRDI82lnWquToSD2U1ilNMllnTf27mE+J3CP4J24AeBAT7su9XQMlmIwqhZVFSV7gTr65ERUCfkFwDfGvvlOyLuQ+12AFyn3/urxHX9okpWryv1u1Q9hfvCNsfdOSnvkeUqzJ/D+GoYM7Lr/q1sHP2TsHXnp6Zrb9z7Ae3idS9f8ZqJvF6dFwHeNve/gaZqCSEIPo7TpBejcoTzyPKQ3uGMnc5O0BHjf2HttgMPEFXX8nYFvyBsOm9fyuLV3A94f3JsvmgLtiL+4sLsB3yOvG9wn92cW1ELgPfLNd5pWcvD/+wPvkecV3KEPaTUtA35wTgdvWPES9R7BO/KyvDNkFyNEFUZDS4kvBH7Y2DNvHOF2CL4hb5su7y1zxX0ItFIfvGy7GPhh8rwhqXUDXsW/55uT24WzzWlYt6cR+SEQ8uByk0uBb5l497OmhNHauAxH+wLvyFPr9/p5nNmDfPWEPHi9fjHwXWPv9unwbKis3+wNPC9n0AZGK/NBKHiF+/FWVupDh/glwQ8YewQP//ecu32BF/LwHdi/Bvx48hyEAVmsLSMP35JfDnzbn7fizDJ5N43dG3gy9zy/9c7SjOryUO+OoUteJEIeXlF6QfD+7FSmc+LaS6afPYIn8viN1srtyV3YsYAkIAy9qPwnQrzd1YMt2+oZe5q/1O/WGokQluRHC97G9iRhV2D+KrdBeyoADRO/UHUVqHfXfdLn9zFl5xYF3zX2NHmneYzi0W2X4G0D3lRN9+2Dh3DukoplQaKvOUMTFwXfMvYuqlRrcmvsfsFbCS43kPZL3nzrIBmG8CP0HJnPHY+6MPgWeRraOIVbJ8/dorexQclepaF0dvQYx6Y45hSnu0yxrGXB+8ZevFlatlI8nd0reIsLGm3y2B/kxG5ZVAx9mc9mYfC+sfcSN6I3Sz/tFTx5PPhxaJ7gUjB2jpHYkqZlsZdfHLzn2VNuZvhea27mkuduh+A5ARTl7sYzmEXdzUs5erF0jYalwbc2aIm5Zc8eJ/Bqv+AtRl2XcPK6duGgn0N6K6XXKcyxOHgvXYLh/He8Q4HRd6rcL3iw7cUxKyDxQZYVqqgg11moJwfJDbG8LOYvH60VwDc+vaza4sqFxeg7Xe51jMdVy7rH0+E6/D8ctfAyf48TgH/++9ehtUaXBy/ktRvn+cQYBtzq3Tp3Ijpqx0dt4N9A9M+v7edf/s/r0A26FcDbFngtXg2eLeL1iz2Dt9ZLeiXsRx41yCyCf/ln+BK2T7MKeG+hXhw8svdwomhHoVenxRU1MQMnpDOlk5Tn/6SizAi1hf/XX2vwgTvya4C3MryTk4dLVLQpbSpM0Z/Ai/vD7PFE/eluDy5SiWW7LOzI/v8/1F/+8qewjdlVwFvx6WXxTksrqNv3Lo5QjZF0DU66WdX/9j8ZimioiuyY5SU+8vL9//hQf/m3jUTZttX49Fbunb8ry5vPiDGbmt1a7PZUeskr78NxLFClK8vziv/q5W/qwf3lb/9hK8GWbcnBKSM/0BG6+jRbPVUAABMESURBVB3mhXUH6nYv47o4bMJXlGS70lK7AGPWcPumbM5UYxDGy7fbiblriU5F89loDsHBB1SW6+bsbJJn3DGnCuRRyXG3hva2oCJbUc3wca0EnnMDcMAVOPJkwuoWDKO83WH14FMynm0H9JRJCyqUKVjTheTmcySAXQu8deM7/dCM+7I3n8izOB5VbDt18+JwOBzzY1bMVWt2NfB98l7wZSLfkqFkyfITrO/k2eHx4VBP3Kt5uK8IvkXe+4oH6bwA8yRrG29OhkaDM7x6vL+9Hn+CvDGJ/AkZnz6uc97iGA86Td6stAt9e/Lo36RXT2ozl4cS+Uty8Afm8XxoMrSO+LrgT5CnJCmJ/HmdWLnDY9KbOjs3KNPy5+WhRH6s+mv1yDw8q+na4Jtl+4a85bPikisg6YKMvzsH+7ER+atXB5/IzyTZj4cD8hFpjNcH70NvrUsbThCSyI8W9vjn1zEJ668BPpGfVS/f/3t4h78O+Ia8t1ArlVp0cIDxzvXpx9DSY6DrgG+6e5u8TuTDhUH1wboS+GHyio/ZJPLL61rgG0PvTd1NIr+argZ+iLxJ5FfT9cA78qZD3jD5hH5JXRH8KfImkV9B1wQ/TJ4zfaaF+2V1VfAe+eYhKsSXyI/Wpx9fRdQYvS74YfLGI5/QXxIei38K3Y6/Nnghb9rkTSI/Wrgf/+mnre/H99QnbzzyJpE/JzyA9DFmxfb64E+Q54iDRH5QirNXH/Gn582lNB0nR957pE0+oWep4thSwY+HH53bAvgh8i4nXiIPUmUxCLzWy3fg0D9tP/RqSBw73HHwHPndnrIxNfCsBbx/qhj35jaWtny8JHZ4gLzeIXkIpG0RL6rT58hhGn+bYzyoR97skbwC4HlDPJvl6MSwNgLeHazwfpaCfEL+ftFDmru6izvieVEtR5y1GfDnyKv7JG8IeEO8AOArbUhvB7xz5dzPukX+fsw9nIIG4HkDHNIerPv2NgT+JHmL6e1vf6DXCLwhnhHwlYmztgR+gLzpkr9F9A44G/UaONj060YZbQo8k/d+FvJ8lvamyINFVyV3caohBcCvTZy1LfBd8o23J6eob4C8B5yIO+DLGPWPuCX758Cdmo2B54m792NTqWfr5BF4VYpRz8CiFxV38eVu+vMv//3DhsqIR6tLXvvkXVWbLUn7wOsunjfA1wgcfPn+337YUhnxaN0KedMAJ+I5/FNWYtTXusfn17BP8xx6mmZ74O2pSV1zovKq5I3qAC984EvHDKmqtWtj4ZQ0BF9tqIx4vE5N6q5LnoArAQ6wSwa+JHFc3statNvga68u5qD0FsHbtrH3kltTjiS7JnnjA6e+Df8uBxyGkApYD8Iuys6rgVf36Y//LzBb/UbB21OuvdutW36g55IRPnAkLsDnI46vhOUJhmDj8t6ZV0Kv7umf/hD6qtsE707R8U/aI2+XJS/AKyZelvRtNSNx9gsFdpt2RhVlx74IZkF5+Tb4pPRGwdNZquYnz8GjuGs7t7nvAa/oe9UAjyeONVUJNbPO2qQxL/1o2HNoq+DXIi9DuAAvETjW82118cBXMoxaWBPszEddlAL7KlOUzYLHLjjs2htcwbWTzP0gcFUxdwIeQBw5V3Q9KilTZo/E+oCkHWpmffWViO2C75OX76aQF4suwBlxRdx9ox5yYQNOINZ7h2hY6trfvH0LOeYfHx8P9QNa69ZM5eraMHhcDG2fspBviXzIjN4BV23glQ9cRcZyUymJSmwGZJs9Zo+PyP+QenyE/MwYLdd+LHnTAy4d3o3AU4jL3RgvHhyzOJVUK7Ziw5/G+FCdI28c+f4n5wGvGLjuEFdRRr0v5u6mmAZLQuvWEZHk1YfKXCLfmdE7n03J7NsBFy9bzUVcbgRvQTwONPpSUubMu5ptHm8//Qhz+M/vAw9Kbxz8GfKW9+itN3lqAZeHOsRnHWEdd/5K3f089oG3OGHlzr783T9aODgZuISzdfBd8l7ftjCaOndcgOsOcNXq4jMvq9MFZWTXsPag9ERjErhWb5//6k+/2c+//LdNlhidolZ22ybWnkdpKMvGThsTd+QXJE73It2dLHz9OtDW5j7ZfXF37um//PKrff7PoZkRtg/eJ2+4h6OXzsHJ0Nnxq/bsfQN8MXfJ7+40SYfXNwtvyPf34z+++/MH+6//97tbD73qC+uuuWlYmeduCC+5CeA3VVWsQ5xvyhj5F2smUv7t9XO3fP75w58/vPzz1mvSRKmeFkOcC1pzowrCrqEwMUKuiePvcIBfaVrccIdXRHdOT9rGidanP/76VPf5pzsIveqp7kto4GmujKO4mzbT8IoNYb0cCu5lyJHTmqYVV8nP9vL9b89/94f7CL0aFM/Y3DROVr6JveEUWSsYea+7aymogjOy66Tle35tX37/4U5Cr07IsTeD7LXHfsm78LBr8kBqn1PdWmGNWwJv2+x1hz39uHS3d3M4GW7qcWYjp6KCdGPgLfU0Max99mrZbs85N42EiahKV1k5V6HfNXV74K2w1z32Vjy9pbq9xq7NM8v66lWR548Ph3zp9BUL6CbBW1meP8FeLcS+7udoVsr6+rQ2XE8z80OGCwszvs4aulXwPMwqMu22FeDCBmF29Dy602uCJ1/kZVHCrgxvE8z1QmvodsHbC+zZGs/H3s0m6Ko17RzCLWBsUbfH/qbB27ZL1zL51vmAM6GXmYRM5nC9kIYU7e8Hr84+eO0GdevgrT+V0x32RETP0u1lL4auAtGVlT975O3AYn32wSfjSXcA3rbZq46r523bTEAvr0A/0K6RzOjxyvRKwL5Yl33wyXjSfYD3d8UdexeHT1t2k7o9rwzSD7gVqL28JgKfN4bx4NXiGQpBUEw6opI06F7AW8ee6PvsyRI3hjni0rQHhN9jv+4VyTLeq1OwQLEUe6PbgRgxBUbvCrxts1dt9mzyo3bweC0Yv4XefqI4GhsdWkWqiiIv5stFW1+z6IbfHI9gguJ8uzsDb4W9ZfbKsW+8/1D2svtjcXBX1dlS5zzgoL8H7POijN+1gyuVPd5Z7j8n0re7P/D2HHte0A2J1xCOlrHriyXu2dfAO0D2eQh7+mM4i9OhnQ+vDUb6dncJvuVqazrf4A5fNPO+MexNUxfHyL7gBe70d25qAeyh45/NcSf31eFd084vrQU/vwKFj/P3Cd4OsjcSg++m/BfYE3Z6lpHtYDPy7KNxLoYh9vXEv+sP0q2UPu8MaNe4l54P3i1427C3zF65kze00CcW+cyfayPVDyWibix3vILzCqHISM0eYsJlwQFTnTraiHsF3k73DN7aZu5OB9grx14mdyfZm0aNlbftfFxjXp+mF/X/4Qxt9vBwOBwe2rTzakXeTncO3vrsjce+cb79xZ72HzXRnK67B3LHS8G8n3Lo5Pmbh2/evPnd129qQcrTaqb3GKH7B++ZfDmfJOyd8+0v8MsfyGE4L266XS9p9Ms3+0e6zGqDf3xMPX4t9dhXyu/2vKdrmudat0zL27707Mm7fM6jS2P8WjrF3vlfur2oa5q1XyM+4czxPMmrX0kNVVh7rYh90+15Y11iN2FB3mGfMaBj+MZi5/G0ZvvybehWza7A2yH2VePcNzEdeGDHLfC2zMXiNxi2cgeH5+BIRejr7A18h6Gw13z2TcvCbkUuoNfd7foHZcas1cNifcx6/f7A2w57I5kMdbPCo8qi9A5Bcqe/6i2f3J17/uovMfs0uwRvO0O2URwux0t6kvauWeHnL5tQdz/+PwVnwwDtFbz12RuEjflU8Nh9WbI/T3sy2+LuC4z+xy8+xOzQ7Ri8oFeyXGMqOGpfHgsAX9Eqval4XX+L3EF4TjYiGGPX4C1H5Onmh6J4hCLOGEvJv6UzFFe9yzPC3h7R5fcOvpbyJ2r1DwelCm1k/5xiaVcsLhQoirUM7/IJvO16enDavWrKEbiBYKPkI5XAo1rzOwiXkw1cfsBu1b2LVgIvaqHXjN79yt4b+QS+kY+eyevmN3dm7hN4Xw16zJ3ZhM/fH/kEvi0PvZB36O09mfsEvqvGz5Mjcm557546fQLfl0NPA738tFnyT69evUrz+HnEsPE8dBOVs01z//Gr3+zLt3eb2XJtEWzy7k0Tfrm9Tk/pi4MPSyfwJ8W79ngqWtBvjDylLYfvglfrE/gzMkbCcbTf3e2Vzb0XksG7cwn83DIUhof1hVwE1nXIQ3a1bgRO/fAnKk3xlEz93KLY+orj8Yx4eiuZ+/qli7xTmKQVdvn5PZj6Tz8m525+0TF75Tq9Xpw8Z9Ro4y6GgwKev/hQ0w/t8An8GBkhL9Xl9BLmHl8FjtdkLdzFpczYz+l8/GIyXfT9UofRl9bC2wMOGVQWjvpJ4EfJUNw91g61HHs/hTw7DoXPGysPr5b6PoEfJ+Ny6njobaC5Z8MBvB3wLBfe69Y6SOBHSg7XKumSrjTJJfSGc98xbwROuEu1Om+nBH6kZCJH6PEhOlk9TN4Ibsx1iLwzPAlbEO+r17BJ4MfKtNATNS3ZNZrnaO7eyJuAU/cW3hupXpPAj5ZpluvdaWqXyJK7t+NdYKID+L4S3tsALkrgx8sjrxVnz4Ms9RWlsMQcpgV/yyVw1+jgHyHRXWiF0QQ+RNy5MY0x5awkY47/z4o277U6OO7S0MJtiBL4EGFOcqpBU38HznkJCeqLqqzyslTn8+YtJNqXu98So5uQgaKSABfT51RYpAIehHTW1yhBBttzGILx/GXoEfkEPlQub4pz4q3buFlHYGu88/EwxAdzT+AjJNlw6+84l70Up1pSuiqPXVkXiBGsBD5KNIOj9CjUChby33XtVHTjL1oJMVO++nWFwzokR4IfJB/ebFfHWWIn501eDiZATfnq1xYkwhPyltPdTruixtp17fiLvLhQ0SqyFlECHy+e1XMMHkfgh1+GcLd558Xi1WoT+HhhPIaSIDyuWzL2b6U0ZeGFXxTIe52sKwn8FOFI77ZvVL8qWff5rhRpE2BFuFdPYJ3ATxKR9xISDwXQEG7iLR08K4j3KnUJh5TATxSTt0aKHLosyE2lYRdByev52MGvvFmXwE+Vy4WGU3moQFKURYs34Sbe1+rgPSXw0wVJsmD7Rjbjj4fj4UC4qXuvFkEZoAR+unSRZfV8O6//qWmrqh7HH2Hjbi3e6Xz8lWQq2KWtgSNqXXd9MuwrcU/n468rir1SVN8Og/CR/eLw0/n4DUhO2VEUDjr11PEX2rmD1ft0Pn4jMh58mdER+9ngGyXru+l8/KYkMfVa+NdGv/Ry6URfVuV5awunTOfjNyc5UYENoCKVkkMrhD7G9Oa5f6YyK7z13drUv3z/l3Q+fksyTVguHqnhKb3WI+DT4j9G52eyzFsD72/6w/n4p1fpfPzWZJrTNTTeV1z7Zpg+7/bQeQxEDidxCn0mwOfl2zSP36bk5JXbm6O1elfZ0j3DAa+J57gapEdGbIcHYCXw68jB12j0qePLuauqyrm+MB69guU/PRZ5rBL49STwIX8aBGDk+eFQ92sy51hVVoCvEKydwK8rDNegbBjHQ/bm67cPDw+H41GZiovMzxmzeU4J/NoyWkb0ShVcPpwOVOdF6vH7kOk5dWmM35cme/XhSuC3pYh5PJ6PD850l8BvVKNX7miX5uW7wIl8Ar91XVqrp7Ub3qsZrwT+djS4O8chGE+huVAS+BuUvx9P5+NTDpx9SSJwgvt7An/zisx9lMDfuijmKrzLJ/A3LvLtUoWKpJFK4HeqBH6nSuB3qgR+p0rgd6oEfqdK4HeqBH6nSuB3qjDwr5K2pdXAx79Qut62rpfA7/R6CfxOr5fA7/R6CfxOr5fA7/R6CfxOr5fA7/R6aeVup0rgd6oEfqdK4HeqBH6nSuB3qgR+p0rgd6oEfqcaBx4zrHz6EY/f8z8T9fItFruf52IWDgq/+uLDjNfjo8ezXQ8yDc95PTv1WqPAPwMk+CCeXss/EwUpW56++m2ei4E+voPiLPNdDzJNvJvpzVp6vy+//zDn/U281hjwH7/4X3WPhyQ7dc/nf6JfkASnuusrzXMx67I3z3a9+g7/9r++m+96eJj547sZ72/qvY039cjqpw/8T+zrsbjHz3MxuMHv/wVM/WzXs59//t/v3814PWvn+/BIE681HjyUuapfiP+JfT0RDVAzXcxi3bX6k5jtevbpBzCl810PTPMPs15v4rWu1OPr8c4+f/nrjD1+xpuj632et8d/+vGHyb20pdV6/KxjPDfX+cb4P817Pcoh9sN8YzKVgrzJMR5MFXr1P8zgmHKPn+diIKzF9Nt81yOvebbrcQnQee9v0rWuNY9/fjXvvLu+0KzrAjPP48mCvLu1eXzS/SmB36kS+J0qgd+pEvidKoHfqRL4nSqB36kS+J0qgd+pEvidKoHfqRL4nSqB36kS+J0qgd+pEvidKoHfqRL4nSqB36kS+J0qgd+pEvidKoHfqRL4nSqB36kS+J0qgd+pEvidKoHfqRL4nSqB36kS+J3qPwCBtjLPdpchUAAAAABJRU5ErkJggg==)
## now that they have different node rotations & numbering
par(mfrow=c(3,3))
plotTree(trees,lwd=1,node.numbers=T,fsize=0.8)
![plot of chunk unnamed-chunk-1](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAfgAAAH4CAMAAACR9g9NAAAAyVBMVEUAAAAAADoAAGYAMlEAOpAAZrYxAAA6AAA6OmY6OpA6ZpA6ZrY6kLY6kNtmAABmAGZmOpBmZgBmZpBmZrZmkJBmkNtmtrZmtttmtv+QOgCQOjqQOmaQZgCQZjqQZmaQZpCQkDqQkGaQkJCQkLaQkNuQtpCQtraQtv+Q2/+dkJC2ZgC2Zjq2Zma2kDq2kJC2tpC2tra225C2/7a2/9u2///bkDrbkGbbkJDbtmbb/7bb/9vb////tmb/tpD/25D/27b//7b//9v///+zyifPAAAACXBIWXMAAAsSAAALEgHS3X78AAAgAElEQVR4nO1dDZscN5EWvvjMcay5cAaSEFgSL7C3GLM4M54QT5aZ+f8/6vp7uqUq1YdKak123gfS6+lWqVpvq1T6KrnTFc8Sbm0FrlgHV+KfKa7EP1NciX+muBL/TFEn8Vvnbopl5gYokn68N9cmC4Dy1BHvMmP/s4Il6vTvZES8dfEFAMpTSXzyuxI4fvkydxYT3OIighXxJlJiCMuzTuLLmvq9e3n86v6nTLydqU/XpQh4hnDrXv77vYp4KzUvJs+LIZ750Mf7H6/EZ0k11S7e49ub0/tP0jygbFkPPbjPPumI321M9MxO/N65FxuLPOXEj1cedt/cpxco91tL8OqbtvOzdYiX6dkU5tNfk/PUpHLelcDuu9tv04lnZsn9QACsU+PlCr/bnPZfJOWpTTV4zmziN0//VYz4BKxEvDjV06vAMuUnfjSKredc1I1x25tt876VO6IKe2PyRctTO6GqbvhP6znnIR5r2d7c7Zvua+3ELy65kqBCFCk4Kc8k9J4zN7vGq9/yR3BgoQfnXrz5tNvkJB7wlqWQeiJdldveNOUpeLHGCU0fuRMRP16lXv2bH765E2vk/3y4k5t6poID3olZD2X0Xzm379H9+fTm7v2nnUDP3ebwv7g8ruqClCHfvDx2//jzX9KJj96FEsiUPLXe8t59QT8WKjX7oWGx+Z+I+NfS8tyGhklmsGefKJP4ziiJHazdh199SCdeSqOwz9ng8NoJ5xFD4jsWX4iIPx2k3ePQMknq7fCv/hMlUg6F3hulrO1s7B2EZl6VSoaQ+NPhrvmfiPjz16xXVUF8/4kSZLrx2Rf/fEU9mwo74nX2SZaH/wPTLnkk7N3LtL6qgvj+E6Vq/PRs04NfPMtpngDTxNZvcUeG0T6xW9DGW5Y18bCuChKaUvVMaUxnYGhZkSerVNzsv0sJHAffingZJvvE7oPI9ER11dS+08dfeMRHdN5tnv47OU9eKrfUAkhrU+PZlZMja7RP7BZ0vRrfNEnfL01prBMD6JmXeOhZTqlu2c6ypUMWmHgO8cIlYlbEh997jPjdpp0oScyTl2pGvI9kvwTKx8zc26u4lM/+MdTKIeSPkEyDFSA+uBP6JQmQjCrwhZqquBTN/hF4IPqck02DrUB86JcE2G0eCBM6feasUQUhOCp24DdJM9HMH4EH4sRHpsGUK3CMiR/9kpi3/DnlM00mnjWqIARHxV7PQl49r8b302AQlCtwpMRP+Snu9tj+7pZY33Bu2zmjCkJwVOywFa+8goknsb3Ze8MhqBDonnIFjph4gWwY9AocP48MzTEDRjWeTvUmGAcT4elV0J+rk3gawRhLxrwsodLz4Nom+qPpCvBLJz53PtZQ1vjDnfUkgpL4OE6MJikKemVLHcQbDeDQqRLNGrBPQUc88fwT1SRRvgzghUZ0Siee9q5gNG28cBGOWtdgAFSkZ7hPIQvxr6kmiWoMAC80opMB8cp0axIvkQfsU8hC/IlqkijiAS80olMa8VO1UKSV7+pdifiwl5SFeLAoF7bnlLrm3ZB4EynS3JITpu4ayEF8mMb/FFovIG3N+5X4xBIsRLx/fZ2y5t03KKnEp/ZBZLnZJDy4xF0DKxF/iq95j+9C9Zs3mScegOyDRPAgDdxhQ1SrtmTXALBBhUqodX0I4gkAS4VCaUZVdOqDKHC4k26pMCJe2p0LN6iQxCsV9j5QoTMXX9Ik30IUg98HkYg1WnolTijVFdigUoj4dsr8I3+VALBUaBLleml7etcRD2EzxMbx9tORmkZEchPDI76tSPzJaGCDCml8m8rwd0W3wWuAG3P63VuroZbOOD/urIgfMcshH/TEL/D05u0fXiU5pNTruqZqSeqqL7f/ozWnhzsz4nvjbMhQT7ysFilzspHy2v3s55uktWE08Y87TV31iPdqFQV88sOdZRszxK9F55cpOHLnSfF9+ngBA0vZSOKbqvUvRV09W6WxIom80HBb70xP4VfE1Zhdi8YPr9Xzb5I1d2aKh6yddQLLM1zKRhKfqqyqOQK29Y7ivKsdqLGFQIeuxv8W+UCJhIkIWYmLBpaycQoySVFVc4QvacpIPPsTPw+TC029aRvlj5R0SqGmHp2kyUc8uyIxxS2vtmC2RLK+6Vy8qofE0rXvMPNF5yfetkUu0+mK4uAYfVP4i9H1kDgYOswVEa8SgC+9qoD4dgXc0Dfdo3MKoH5O10Ni6dR3mFFTDyRYXmGhqVpJgS+9qoL44dLpGX/G+1HXQ5LpBGEt4sXAl15VQfzUeOFTNDDxkXtmOkG4GOKBEJwDaiD+jEZPokla/mjq7ggA7PGrk3gcdRGPAx8JrQRX4vPgSrwKD+g8d2XEH7BofHURHwmFQhAv9svTgJ/yM+lbBWokHtQzdJW5xOdREgVJfCWojPgz00sAC9muxKfgAV15tRLxcC7HL/HjxyojHkU1ihCoinjto1fi5ViL+HYZKmsVUTXEL1qnh5t9ZAVOPcCnZdcgfu7NLQGUZz3ED5de7y227VypSBafPnb48SrEo30eoDyrI7777/GPv0amP7TE65KRQA8/Xpl47w5QniziYfthiznxx6/useBHdRFfmalHvXqgPAUjd5mJ58WKr4t4eYa5iWdH3K+D+NaedLtjyEyuxMfyHPYYmXXnLJeKIVocXLc7ZnNZpp5eKVRIkbE5fj0rxTPUpn5/k/nkvtZ9HLYudljTuZN49fRKoQRFRHDDf+eleIbWuXvcZT65r2nglx7kj8ZDthLiBWLplUIp0iUYiMf88LA8KzD1QwN/vJ3lcCHE0yuFUqRLQPCoJT5rd65v4L+jg+xr9eCns3tTinhJg8JCY+G/ltRO5rB+Xme0b+DJIWZt74Kfzq7/QhKfnoUnd38jOluhfZDxnecm3ssFjjRRH/F4rJ7yxD/u0BP+mlb04Yun3wSKMF43N/EeYJepMWcv/qTYF59IfMS+4rF6ihPfBbCY701d6Pnh7VuvmV+ZeMTYwMXp2k83GikVg39UF5oO2noYKR48Bk5x4oOSnOf0dHvfdMwCRdYkHhQPrRhpnzp+dX+M7Q9D7rj9DRFTefpjOKCboWOHSKweIqscBQqq3V2Pv9/sgXBnDBNalnjs4UbPD7FTSDDiH3e8mMpj4MW549ObgSZnfhHEfKbcxC/VjpxH1td4woTWQjzlhGLvSPV0zsSfDnd793KZUTcAfrwVEB/RZkEOWyIz30DtyHlknSKUCc3Q7RzRdz/FhSq6yf9gvKIbLq+7UQZD4uPqqhGqvcHOIzub+qgJzYeh++mW1QDqfkwJotLkd4K7Qe+yHwAHTT08SRNrOgsS36797olvmvjD/wCTNJksDwdD99PNSqMF0P2YEkSlye8Ed0PiccCTNLGmsyDxM+y/CPrHOZ0NBobup/MUAbof89uoNMWd4G5XBTq/iC4TeJIm1nSO7VuOyQ88Ds7xbbCTf23iJ2OzcEKB7seYICpNcQe6O/hFdJnAJ2lQTWfbvuWY5wY7oz382NXVEO8c4YROCXQ3hcT3fpG2TKiuStu+5ZjnnnVGGQ9718I4Z084oX4C4c3IHcDDGfwiNfGEz5TP1PedUebDy2thzIg/O6GgF+olEN6k7yyI53q78dCraJaZvGlc7G6DnjtXAfHDX4gX6iUQ3pQRP+lBIR56VadqGjA96yN+jBK5qAaAFzoliEpT3OnsrjJoNR16Va5QDlRIfO+IBsUeeqFjgqg0/A6O0/6X6tDvdOhV7X1bNL0PzTbpjBAfR6UkPpbm8WNi6HeNKiuVd6jAaqa+c0T52dsTn8fFvhJPZA97okCbNCaIStNoAKXl+t3AsV5MVcqWNxBMamXiu6xZXuj4cFRSSu6QINKrD4/1YqpCSbYFEFqmAuIBFCQeTMctFeBYr0RVMuFSiAe80AGFiGeuUgeO9UpUJRMuhXgc+YnvbaNwlbqZKsVwJd5P5079NAq6Sl0uskZUSrw4cCB9k5Nu8oaCVeowGpde44uUR1ieM+LLOppRrEb8cGF79d/cc4gvUWRREMQriy0H1iZ+KrE4dt/dfsshXqmXHS6FeHmoUPomJ53z/k0BONYLUmX1kgXKs07iceQhfkI/fmtQEnURD+BK/FJAd0DYlfi1oFvZYkF8d0BY4NthzzOnZS0+JKbXASM+H5+HeJUX2iiKHPCEK5lUMqOM/ghlv48b0xOTxFOarxmlShTAPv78xGsSqYhPyG+Ssfx4yEGOFYhXpV6lxmsSyUOF2hDvCSBLGzjWC0ry7Ih3Xg2yEEndSSd+Bv0iDSnxpLVOWRoYUTADRafUzxQWGUG7U9F0JU3Clhcx8bQq+qWBkQzzDNkLpZ6/PuG57H3q4+33xvuSyC0vu80D1ftg1ij6fsrSQOCIn6qIH6+Nt/w3LMIhnvqr+7SZVEDkaOox7D4nRxiZLRB5n1IlCi7xOukA2riqW3YtPCuxdb9FNipEUmcw9SOwB7a/u6XW/1sRT6kSRYT4h3mxWRVe20a2cUQUxCtMfUrJKMEYq2cQn0YqB7EaP28h7Yh/bCicXFGOgUjyXCsae5y/FDeyVFH1R06O8xbSjPjW1B+6Gs/5qpM914qInzBUqWqJz2Pqz2xzXu7sue7RUKHx/DRKpgAPcDhhqFLR95YGIBYDP1R4WSFtu0QLkVQB9J4rtmWSl1cmAG0TFq1+nmpu6pEmrp0T5B0owtIqBDC0DAxQrkf8VBDIQD0zrwxAWioski2Qsv8H/MjpcQccKMLUzLuCiBwjfspG/NlxiweUnePpFbr/OJqXPIlQ8jIHJPQqmhojPsXUc4nHTpo85SJ+EkkFlDXJJLdkXQ4U8SnCZTUK0Ck78VRAWZNMckvOQzyjx4Ok6B0ETY0qRnx2z/Vkr/eseJPUJ4iP34k97fQ1qhjxsq/6QTxw12WiScSSB6sfOR8PlpODeM4nCawbCIiX2x02WD2PDoc7lV+fkfjxugR+Ph4sJwfxHHC6c6sMI/lgDIxAKEC89wR6TBogJ1qj+JWig2z267QNx8OqJP54++mITXvFkJF4pKRFxPtCUxQTzX5BawPnxI8f04vPc3thmZBj/GEkTVbSkF6mxC9mv3Qigh+Ot9/mPVA0GzKMOE6kpZS0PfHn2S+1iOCH1lVN2hpuAPyU5igyEg+X9MPNHll6tZQ1wMyULs0RDaA8Q+K7SQW7ITah19KjgrF6v8YP/1ritOXwPqa1NaWJ5RkSL/yYSP2QfKJ4t9lr3HrVR4ZhsbQaKRVmd27Qbj1TCpQnqEaGmiMTeXiN7lRg5GUCThgc7gDOINHWlEoAlGedxKflZSTL+jQBY1OaCCnxKpP5tXtxecQrZ73iQi2FpUFMvDyD/c14XB/SUoaoYOSOpSseiBERmqiUFvH5eOJH8haW4HE3HteH+MYh8F2oRF62OI+x43qWJ15udKFt5/mJn5l69nBwBTV+kher+sWJZ9vMJVYhfqaphHhBB/mclyINKS+mNB56NSIwWSO5IGCDioR45ec2S5hzNYZWuZjI6T8nm69qLeIjkugfDXJNWkdMC59f7ES2X+sP+HHccoEWAi6OeOY6YjzSREz4/GKCwYTsb+YnIc7sXuSYtJjAFChDNYQrmkTEp57iMZl6Chqv3mWaTu6+1tlJiPMPAD8mLRMW44n8zwVY0QQTjyA1mPdZEAFgxQgtPNN0cjcxBxF/ih2Tlgnn8URuUXaIb6gQ5JobqhqfaQ7E/1oXV/SYtEw4jydKWjZoRZOI+ByjmGZwGedA/Gb2a8swRCpFkl0aWY1PzCwrSs0GDWPQ661N6zN2iX0NEfHsdoX/JAwgggM7T2WW/Hz6Mei25Br6H754+k3uHH0Feiz7GrEnwQgj0nLiPc8elEegIZ6tXCImU99g9+HtW95CW3s1ln0N+JnhAgWTyky8UPqIqomfmZan2/vG3S2RK6DGoq8BPzNegGBSV+JTcPz9Zl/Yrx/BcWlmxKebeh4MhzklypXNblVwacAWweQpK3+YsxCeE/EczBfBhPfy5MhwPSKQznNP2RYGcFhvTZgWwQDBpDIRPw5zKgEE1mdlWxpaZ6QQRlMPBZPKRPxUMjooa3x5VE78yAJnQ0UVuBJvCyCYVJ3EX5EdV+KfKa7EP1NciX+muBL/THEl/pniSvwzxZX4Z4or8c8UV+KfKa7EP1NUSvzhbm0NRkTnm5Rh2cojLM/5+5Sc0IyjJuLHK4T3n05Pf42mrgMU8ZYFloYKiYduvtsQxxHVUqSXQvxDPQtb3N69PGLLB59eEUtwKilSoDzrJH49hGZy617qlw/WW6RX4pcILLtrV1uolw/WW6SVEq87msQAIfEPkZ0L5NrASoqUCGKcS0u5F6o8msRCV+8a9+p3G+LIwRxFKi9OKohxNuLF8pXhzgzQ+nLbeWCNqfCgp0k9sxAvl00EMa6GeO3RJAZofTl+1I/dpnX4Y+JMlIJlCmQTQYyrIX49tBX74y/son5UQnxMSsSgqWWP8M1nvRCegUvLY+Qohklxep+PLfHjlW8+z2+n2UmTjnP+rMejETh5ouQlPhUnX9XIgYNZieebz+GR92sZB7e9eeDEYxuw2xz8feczWd41/pRAR684WXpS2pkV+PxzFJjP4ZF3pbfZTnhzJwnoto1NzmUj3itOVo1HDxx0p+5btyN+dhXYpGFk/MmVne4cNTs0+f5TEDMtGpYtG/EjYrMIS8RqfHuc3v4mC/HdlYmkkfFkbU/ucGdzyn33nQ/tRvQhXoUAU5rMIrjXTR173GUjnp8sYWRcDX+czoD4qS5FhKXlZjOL0M7V96aeWzcJkJ87rE4/Mg7EbMmKYJyOmQ5fgeOmuhRpNxKJj80iLBEJaXp+Y5vK1n7uiri3w+M7F4TnyorFMJYA+AocN9WlmKlvasdW7VjZefVu8Usa2s9d4iaNyUYE4bmy4Wyj5MBX4DhPMExIUzv2TutYCRSP1fiTMfGnFDeppKlPGQF9eoUNNPlWHJTe1o43n+wcKxSRU6jMiU+oR2VhM/SNSCWI72qHXR9aglzE20rKiVWJz5AzE3USrzyFClCBgcbD2p8bJZGp2t5gsXp4xEutotaOAjGFshKvBnBOmlIFxiONhzXzRkQjTrtv7lOIZ2oYPiwtT2CJ2JL46UmJNhlQlPi2u336eO8TzxG/++7221WIlyQ6MWp8nvZODjOvnkP8srstI36DHSN+JX5dcIifWbnu0g2qpK4amRMftaMym60bD40oOBHfuCvbTCcCrgHmi8xKYRhUSS2CBfEyZWJSdeOhoKjZ1XVDSe3UwtrEm0W2VBDfD6rwyravJhGJ1sTrxkOhAJxLpRrBDRJPWJK6nACUQYwBXQSPjda0G1ThLcLZvfnhG3h7Zybix/FQIWji222Vh8Qab/Ceq9T4oQhGcNLt/vHnv8SJ70VZNaAi5WZgEL/8VYVaiBeU0rxiStTeffjVhyjx/d+rN6A08fGyYhoWG784GYIPsNe7VdvqIMEF8RYNqDXC0omox9PcyC9OhtTyuDdv//DKakX9gniLBtQaOYiX+MUZIWxaD8797OcbK7UXxCtb5qzgE8/WfvKL9Vqp4v8HTY6waXVztXlk7TYPyGSSn25d0on5eP8v8EkG8enQEe//W9i0Lqh23o8wdp9jm2XzEa8pT2DuIwPxBrAhXt60hqUQT7z9HbapNyPxijTPjXj4d5YIXteEnKRB/62HRlLkpMnZK2NYrFeoEEFZL2w3T8QZrXfwUb/APw/xildCRXnXyJNvjDaZ5AKkmVrbtmvy3Vu9gEzE20ljizq4xXqFvFCFCrUlvvEO2qiAlKnnKrMW8ZONiK3AIURQ2wPiGYugWn5jTDxHfwHxNjhJ29v+yfchbQLix7xF0LmEQLAedlbkbzxhrNflE28Db30gJ0X333dhn1ZkPBRvoyMeCNbDzor8zRBApImsGXvrAzkput7Jf/wnGueuLuJVWIF4HJmI99YHMoD1TmTE61qlrXtxJd5G6rIBYpCH9k4MOwhwxk/dnFcnnu8l4CY0lhXzN0Pgq4HzZezCP2PEd70T1gCOLZovrp3zcrM8GGYiGmIEzQr6LS/26PdZD/EtdhtHRL0yx3xfIL+9B4L1cLIqj+OXzNk5O8w/vBO5dmR61N92np34UcX2b/bKHLManxlrmPp5Jk/MtSMrmPpO9oi04W8ym0xyVShD/Gv92pESxA84uLThbwrPkPiEJS9tmrk5zgnHGf7uoApbnr+snVdSbUQ+YFq2VHkKJuvCgx/cqdwgS06vvkhhBwP4xz/+GnBCCw5aMfstwMEPRYlnkyOv8UXeIWwXf4S6c0WJZwEJflSlosKIGFfiI4AOfuiIf3Av/rRSDB4rOPfZ37Nv4mDuUr6I8uxrfPM9FFpikQq0Gdvf5OwqDnmTQUr75y6hPDviG+/0+LYi0xRp4zEl3eMuZ1exzyMIUvpwsweapCrLExrAaU3Th5V2+MDVFz+xFSXeffYvXldRDzfNik56biFXRFqeLMc8CcCcQsl+PPrCwR38NDeU+PhtE7hAYaQ7JyzPAgUfzikUHLkLgOZ9/FI8Bp7x4z3z3TbxTfN95h0ewInpmfqsDiuN1WPQ5E3OQNpjktstauAdCVEZ8ZFM1yF+e8M8W2OWRnMrCWfiuyaetxnvSjyS7Wg6//0eMNI6E1qAeD8jee8j9Vk/Ka+FA0KvrkV895+2d/QjNIIPO01TQvGtJJyJDyDvfaQ+CyelvHog9OqqxLd7WSfiF/fBodDwMe4tNfCORwd57yP1WTgpJQEIvbou8bOcKyU+KlvR+0h9Fk5KEh/u6l2V+KEytQPgL5gKVEU8I6H1s0OKsyWSFB6UaUnil+Zz6B1XSvyD++xrzeRPZuKnP0SFB4koSvzy0veOFwrAY+CzROCtHDjtb/Dgscffbw7/kzaAQzgRWKrpD6DwAAChVysgvu8dL4sbHAOfJRLeUsOdHnd+8NiZnvsvMLeeTbzwee9TgQoPABB6dR1TP587hz75arpzran/u7c9dVZkx7dYAM6MxM+vTHsBhF6NZqyzjjT2N/EAiPUM4EDFM/8NiV0tIJ5/0nZEJQpA6NU48QLZArjHnToA4lr9+Dk4I82xm0th+xvejvd59nYnKZQlPmGZ1Bojd+ci7//djTSnEL/4+3HH2/E+MWVzWMEqxCdIxxItuTGF88qoH2m2Ip7dX5y0sIkYu06Np0jCw5ajxMdvpyAgvhtpbok/3n56+OLpNxF9YIHA3wLim+yDODguXqbEocJRDW1BeH91Ez/T88Pbt/DoMpt4tqmatIBSBDeD8kTFSbVPg1+YHhTE51teHXGnnm7vGz8VTBQTKHgWfipIQRiOiyH+6ZV46RXVRUwA7k4df7/ZRw8jYt0zIT7a0wAijKAE8G2QBgTxWCJco5QuIpmvwp3KTHwATk8DFCc1HolwmtmPWDuecycN6E6Rafj35MSHNzk9DVBcYeJPwOzH2cYg3nKU+Mi9VODGT+6LQPcsiJ96GiCAAxVw4h+y7kYDZj9mninsLcf2pGVtmGbyU3sf0D0T4mNqQgcqRGp81t1owOzHTBXYWyb2pOVzRVGUIp50uIgHgAMVIsRn3Y0GaDrrNcHeMrEnbQXi5b0P6J5k8EYHydKrtkpm3o0GeKbRQWhiT9oKxOOoi3iJyLzdOVAVqtdEOHA/ZeK7htH0yMKVBnDg3MZeE7KkifgUVyD+QbWhIiCexmn/S1ZAO6S6Aoc0V0X8cFUuaWKUnzXwU5MkxNNoHK6PrIB2iHfftPH+DoCaiJ+K02BJUyGUIp7bt0boRIIfocqsV8rJS5oKoViNZ6ZEBkR3m3Ygny+yslJuUaFKIIyJPxtD8OdZ4xNbDg4pcSXeFrbEQ6nAbwFaDh4XVxfxOm95BRzukBto1wOsuhwgnUCPeGA5+Ak8zo3ux68Bou1cVbc5FMTHb0fAI368LtG4yk9/hcUBQ7ZqFQ1AEF9QkzjWJB5y9M+f5ALvNsFO/ivxKXhAQ+4WIF6w5ujpVdB2xsbq2znQzEFCxaiLeBxo66AOdxpU7bQ1R0ij4LrlEN/nWsWmx6UQj2Eo1nTiU9cc4TMiTZffYMOGDsRZL9UQLz/8eChWE+Kh1hwGEA0cJx429YluMBPEsV45iRfpKT/82MzUj4O4LOw2B+wUKjQfn3i5xhrEj/XKpcRUUNwE8sOPpTmEKUcBkgVS29Aw4cTDKpYhfiVTL+7JqA4/NlFftkAqsqECkS57vgDqIj4tmyQZqQukrsTPhYujFOiyMZRxScQzWlL8FCoF8QJnbYxSwG3sVcek2RAvchYiY/WwdNZPUuBDBzMvNJpYSLzgySFKAUNFQk8bhQg5bETG6rkaWhKPPxIhXvapE/n4T47LGLiGdL0aL0NkrB7EasTj58crdlvKn5QQLzwmTaaQFZ5eoYstQZgRvzQ8KWuFhcTL7IONinQ2WaQKISFebmVBQey1wqAk2eSR0vlNUpEhvQKIiOekoQUx1goDuzvHxLLJIzXxvOXMeJPEUaocgBVNBYj37QRjmBnY3TmlFk0eaVVmjoRrvXpKrjWA9Q1x4j00VvZz+Uyg83I6q4MC2N05pfZNffyVddGnGSr2egY+E0t6+TovJD549nj7rXym1ieeAWB3pydt9h3F897fMA4F1UJd460VoZBIfBtjlrV2d5msqRZbs7IPKiNB/ONucShoDViBeEgJwbOdlSV93cDc7m/acPqGL+vQfwRPag8ayInLI57X8vm33ePW5dniC2a3fJLZWBfFxRHPfD4gvqlxhzsRB8CJCmgOpDw6Q/XngS4GJhQqTTwQsqUI8YtSjXvhPYATFdAcLIgXqLaAhvisRggrz1WIH36SdKmBExXQHCyJJ5/0oCJem5lAuC97XeIFqxyAExXQHCyIf1B2PC6FeCBIUyHiO9PGPYtBlgNNPI2Teccjoo5uGIwpfHFhPCmVjN7GGzD+WQwyJRh59D4AABQTSURBVEzkmXc8YpmphsEYgjtwe6/WxOOPSbrUwOnHoBJGjlJn6g+KGo8HOIxkphsGYwhu/z8bqDyXDRQitiDx45UB4PRjUAmj9lLt1QPRpOjMeMNgurdws4HK0ew7MESsNfF4lB5J7QROPwaVsCJeazhUNd7GSgGCO+Fnqzpr74EQsWLiCbBjsEQRmaTplJj+yOkoMaAhvkU2x27+Uc3mSIAQscYqtB5ctgMDzrmcv7NMjlJmZPXoJ0S7KtbEl9mSMGWXy1HKjEwevXdAUbSrYl7jTcRytx9nc5S4eIjMKESQpX335bqpqwKFiLUmHvz0xAA2AMDZ5XgLgZ6HO/lGaXN9ZwKhAY7O6IchYjPUlUCkPA9gAwCcFfSdpUAoT7ehQlMVKPS9KbC+ncAQsXUSD2wA4GeXAv6QZ4vj7acjujiwKIbeFKp3OMmdm3jzOhnPzkjYpTmLY29KoHd24vPlk0NyazJ/+MpiPqEsht7URRAfMwWSodD0N5i3lY3J/Pd7kvhZAs36ansQdjUSttxSicU/+ioEET9cIACB9XnZqeDmfz7uNj/SxI/XKlgfEClQwFXOQbyX5w189mWsOZV4y7bEt/N0AuJVnblMcLhHBVSj3K3ZUIWExAOB9SM5pMArKaYrOj3+ShrnLiMc2rAevwyDSWUnvq9CQz5LY2Cz4j2R+FAI2Ph46NYWVDZUjBOPPJ0VZ1V8M+SMNjeZE89K1WjedOGrIx71qKCn82oz0e1/jS6yuUkSKpRTQeN11x/yYuTZRplrCxiNx7cC5m+1vNNOZWEnTWbVaJHTuTkd4s0A4I7VW+gWH/JCUo2mfvM31ZaKXPBdlRHHP/66wJAtoo7fLsIfZw/uWL2FbtIhry7VhN/6MWJXBdrGA0dzZ1dmKiZBuwgE1s8E+ZDXHHjo1XVQGfFjTkO7WCxHDmJ25/Ig6G0UJL7GLlCH+jTSQmBVS9Z4fo66+HFa6MugrJ40MKv6cLNnH1RgDkkfCwisX4lmvp4n3RKcTJisqofTNvg+6zRzQGD9KlEd8YjLslJ3Tg5dUKHyqM2rR7DWAM4VFeJK/DPFlfhniivxzxRX4p8prsQ/U1yJf6a4Ev9McSX+meJK/DPFlfhnijqJv5AxcOislyphtYVKP4/JBH5+/IpKAVCGr86NUE+jLVT5zQR6fjyOVWxXpcR7V7MtVNnLWGPqr8RP8Ik320JVpWNQpVLrgLPK7Ur8TwpDm753L7fEvsQ6idcsYixD/NlfOsHBgdeFtzd1+NVsBY62jNnesnSxpf+m+bDcBwQEB14X3m61EVZr7tTEc9NLF1tKlm6nYZkTEBx4XcDbFM120mQnXrrYsiDxTfN5nFauA8GB18VA/ELJFpdCfDWCZ1kM2LqXcGiXKjAQz1BS/ALa5pTrbargTtubbcbtWed3dm3fHQztsj4mLTlKyolXJmyfb79EL5o9/B1JT35wT2/e/uFVvmDGbrosQrtAwYFXxKTlMv7MyWYLVRLxzZf4C4/48ep59fxwZ13q1+5nP9/ki1s/I378pbdgYXDgFXEmfvxhqlUGW6g0xE8KhNHsYXHS4MDudLjL0oZMuo9N1NJCAcGB+SLt0bZ3rZZLJU+aLVSA9OHgVxHxvrTZrcABbSEJd4YItoLr5ftN1IDYCbhxkRng3ty1kweheMUATni7aU7vpIfzzat18JVaecmZirMnPmiikkXa4+Dci+/ecsXLiX/d0iVzoCPNg52XnJV40wNXstX4w93hLh/xTXN6Ig7ng5qHB+RQMt8B7SAJd4armozwPTzoTqHK3+n0IT9UGCCe0Zz699rmATmUDPbqNeHOcrXvJ7iJ6lAX8afJPjH0FBPPSBQS/9q9+Cd8WBT8HWnCnWUmHkSFxIPluRrxXW9L0EZqwp2tQTxwPLdAbAZMffhlHvJjxI2IZzQPyQBtcSJOY8/YXks7kQvxMPHgk4rbUuJZmqQiR41/wnrGepHe1RgN8cxvtU7ia/HqXxM9Y/nxY7mJd+7N3d5142tn46I4cNCQeAmq8epnPWOwKVC4oHnnEQ/OvRgO9511mOQHDsqJN2nDavHqp5dB3ko2ldSJzDyP2I7ijNMK4y+aAwcVxLPkEpAcOOhlbA03F+7lAa1Xp8RlnkecaTpvVcQHDq5EvAYrEK8Sl20ecchg0SdBxktPWuJjaLLbq+OAJyAv8Vb9uszO3SIrfLyUVkCuX5vdNpX4fuZXkKmJZxETPfnKCwArW3jiSsRJP4+XQuv/7Yl/3fbEPkYG+RnYvfnhmztJpmp1eWh9Zdi6aYJXF7KG03gpuP7fnvgmuwfodB9JM7n7x5//UhPxra8M5QGsbOEIs9CIkc3ZuADr/+2Jx0yZiPgPv/pQF/FgHtDKFp6wwgDW/5sTD6cavL6t+KAvSBCA/hCxvL5y1EkWCTNQKBmliO+9vsZBSiUe+bX5oG0H1cM8ok6ySJKBNskoSPzrbjAxE/GPO/5yM2XO+KICqSQDbZKRZfR1mWIS0w1dsKQ1fTlkgQNCfGOD/8VebqYDsKjgYtp4IEhTlkE4B/1L1J375l5GvPob5QPSU+nVFwcQsqUc8RI5u+9uv5UQLxRvB9Xm+BU0FRKvdZV94sUzkbsNZkKvxKuAEw9aCK2r7ImxnImsjHgVqtA0ZoedjatsOhN5Jd4IUeJtXGVwJvLsd0DAI2LURbwu9GqueUS8RIGlbFHi8VsihaAsok2MMyY+2SvGoAm9mo/44QIAWMoWr/E9UhWC5MRdfWviI9qlQRF6NTvxwC1gKRvR17LSMZATX7K+dVhVqov4ukx9twIGlA5sUClFfADlkvW6iNfBSJ2gSN9ItisVIj4UTC1Zx9Ip7vxUiV9eD65bAXNfO/Gn2Gbu3eYBaTvrIl6+8QN0eHTwqXPdChiw8QSCSa1H/HAFsfscW1WfkXieI095y1QmIpVEoiJePRBMajXio1/+9ne3yORHTuJZTy0g3/iRTvyZTi8aUaREI0eTlCY+DuOx+mzEyzd+pK8VmlIKohEBwaTqJB6HmHheo2rY9HoCfSStFVqoqYlGNBPlXZHb1UBOPJmQ/5QAmKC0tUILNTnRiJii8NuFgR/rtQLx8bYTMfUo8UlrhZbE860UcPADj3jKu7UGfsiPnHhWoxpbQivzlgk900xLoCazPIGDH5jEq9RMgIZ47CPa30C7n7x0rZe0h5+KlAEeGAElflJLCFLNGICDH346xKMJHndvYquiR0fpxeYRfipSBnjo1ZieqtIk1YwBmPSqlHgUcuIpUz85SqMNDWwGPvmh0zOB+BM+OqeRhyvzEyCeSugWTwGWWDb5QWVH3IsnMux81Ek8HjhQ2TZGGlXQQ549fXD45Ac+Vp+JeKWHIFyBcypAPOaFmhFPJnLzP6DiGCc/AOBj9XmI7/5QQLgCp/9Z+5Ux3gWVvCLxYTT+yPvjY/WmxCdToFqBY9iuhMIRyXio0JzEt5BF48dDr9oSr004AphTqJN4Oo0qI/pm2vg3M8vyxMtlOkZ/KCHvdhxKd+aBMCP0jmdD08a/WVleCPGOHvrS590a1eNt2IJFDuu1Jt57QFjE+OHHGuJx56ype18DdY/d9gPBpOg3JYe+1OjW3X019ZNmL4of1pubeBnww49VxOMJ9jdQTAag5BA9w2BSDOIzmvq2o/QCUAE/rDcz8ULghx8bE/+4g2IysO0TEEyKU+OZ0sXwSn2WAX5Yr474CEAbyoVm44eG+K7uhYOHfOLDYFIM4rkNiRiTYEFUIWMtHGxDDQQr7uHEhxTMSk7XB+F9NPYFM5Msmmu0Jh62oQaCFff4vc7hn4N3rDJYNRAfzDVGDuu1Jh62oVwUMvXAzZ743jum1QdiCq1J/GTAgsFw/LBec+KjzRjpLa9C/Dn/3jveUGpCMYVWJX64hKWPH9a7jlHGnlmJ+PHKdr+AmEIVEA9JR2NX10W8ZlcvfI/kECI+kscSwD4FJvF5gAxIRTURPKuSuNRQ7TaLiafSeVolR8/P10vn5C3vTOUnfvEPtdtsTrynVXIndF3i5Z2pwsRz3WZKKnlPMj46BIe9YOK7l4VaN/zkh6zEB6pMbjOSFt7OTTbXEARTYa1Wi+WfUBF6egJClteSWGq7bMSw0DJ5ifd/opwUeDu386TwtOBPhQHfJ6FuZcT3+UIawDFiqe9aq0HwJzsHeDu3injhVJiD/nEpNb7Pt/VRj8t2FIkKrSlQjgazP4XuMrydW1fjZYk84oMiXALodlZAfOOj/vs9x4HKRPwMNmcRNLX3xZ+EZ3GQzTT2/Kg3rwhnArxrabh+qP7H9YifZ2DhLrcKNj1AdhSieTpthuwiXGaWp+3k5t8tc/uRZeoVNUmozOQu82ogfKCCa/Q/vpXzqCc+LMIlgN6HOxWoSRHMnLsZEOdOWZMEykxX0lHuvWXwQIX2A/2gWLCZQDyhLND7WJ34M+Y/gwvutDVJrgy3CYQPVNBaUD3xVI5A72Nt4pEBSJh4ZU2SqTNeCEe5gyZIE5l1BiCTNNnbzij4HlURXyTBUV6IKZRGjb7GZ247CQ28AUji6TLQOMqL9IXSqNERn7vtjGsAVGR5iBFrkI5yB/zU67qIB6Jsj6Y+e9sZhe+F4vEiyxEPa+Z5y+ip1yriswE4UGHtfjwMLJZUQeJZhYKfel1RYbYID1RYe+QOxPFL9ACAqvSMnHpdl56IqZ9fa8dVTyNcic+D6vWsk3hdUKFiODtN4pCmK+EQNEhliJd6oe9U69WLYSo0bPX/XM9Eh9wEqxB/zp2Pd5u9IrZMMUyFxvlAq9B4FeIVORxeYzsVqijG6ZWeXtFNUg0aP4TVqE7iaWmZQZjOboFcdF14XcQDKEB8U0rfJ+xI9aSZSCHyoNomd9r/kjgf+ydKvMyz2N+MsWC57b1uT5oVyGlqd3r82J2PjZ80uTrxi5IGeh9K4kXPPu7GWLC8EXDt9mMr0MSPa6F3m7+RH+hqk579pQO0TMi7iqTynp1tjeTmtkUPdypj6om2afpAt+63ZNSrlYnv/gv0PkrU+POVm9sKNX5ucM5tE4GaTf38AvQ+FMQz22n/8Q60O0xJU6eMCV2+kju3TXqRwF/2iLaX8ZLWEC9OMUtKusO8zG3hFUJiDLSFyAI2CrlLlHRp4gd3mHiu8Aocn3h+Vg905I58xEflu3lJ20zLphA/usN1efVBjeficIcN2rqlpIxaY8TPSxrofVDEw62H1hJyC0JQ4/ksRdAG3doq4jRTemYdHot3PhYlDfQ+SOKBn/Y3mrPPJBBEhbZQox1jEhxPMOF4++mIHHudi/j558rufGhMPUR84/Qqzj6zQSbiH5uiMQ1wmYN4z1ymdT40xBs4vWrkIb419RbHE8xEjm1Iu1XFqKy8Rj2NBxXxrHY6BehBk5mI174S3iQNgpvW4HuzQx584pf/jCLsfSiI52WW5GwVJl6tb2Nr8bUYrdDjV/d20bGDGs8G0PvIRzz/0RArEK8CSfzS1Ks+rhn6zofijYHex2UTP5aIJg8L4GP1PdxSwUQ9h86HXArU+4jxgxcrh3ht35gQC/1r1VpPwYF/6iR1nQ/TdgMmHr3FsUvKvjFTX0rDemBIfN/5WJN4lmDzvnGoTJPB56v1LJkwJN6yYctHfErfmDxbdiyD4+23Bk6z0xYp7ovMhIN/aqETAXQ7MeIHL7KtUEriU7CldqiMl6/unwyOlFF0kHqsQLwKQO8DJf7UDz80FWoFQ0rX+OHS9pbSF/CqeyDliddBSvyxq1A1E7+4qKE5MaHDpRAvmKTpiO+GH15W5Tr5Nb4H8iwbp4wTjv7gTY48FIgSX5eqHYI6HtGOrXi3AjzXhCM8wr46YsRDv5cBOR9vTLzW1ONNUqBFhgLlGxHBkK2lniZe6FIfkniRtVJ79XjvA9XYkni22oJ49abEi1MYEC/JeiogtoIDWDV+wIP5YBOfo7VqvDgFuVHBlng1OMSfRk1sBpvmItktlGQAB/ldrJy2MmHyFpcTTny/nrcWV2rqG9sQP5ZpQmckO/FGcjx5M+IR9GfaVUT8g81g08x9SeqMQMTMa2lh4s/E4UO2LCMynOqUmXg8SFMAO9M3ffsSU886k2b+m8WgmOB0peGx9+jjjvktDeEYxOrTHvIc0eU3oPIGOBPP1pnp1S+IT4XsGNHhsXfpxDOeCdJMOnOBB2mKKYVnzERTz1tnga8vsO2cIj4VkoP0zlk+uUhQIdaaZTGHJ81740Gaohlwf0eebvw54Xm3sRo/h+Fq8CkCDhMP4+eMyWvJZK5ZluqfYYwFzoD7O/L040533i2Zp+Vq8EVzwXmc/pzZa5bFxDdf3Q/pRRrJAPpNbJscUTcSdDFcDT5/Jxbx3eccXb3KNkjSprP96rrjSNhU4EGakLfDf5MQL27FgN4HqIvpxh/Ruw2m/v8cFlRIJk6Grl98PpWC8aXg27mRDPDfdG4oE5EYOKBYiS4xNfkmdMoaCyqUQb2z4HaZ4Ix4OgUepAnOIPJb1hEHoPdh42/EMZlQdgpqo0KeYvK+KEYe0hoPQXiMtApA76MM8QsTaiMyE3oyLLwnboaygQ7DfGW/q/JYmlATkXaiAOHyjnJCZrKBDsN8Zb+n5CGQiQYVkouSg99RxkOh8DNjhvpNcmyARaFFiBcDDyqUQT1A397U0xC28XBmPEb5XQ0AKxEvB1mTMhPPouJkUuO7DPnPXFSNFyMSVGjAut/lBCoiBhNX4tmoXD0hePbadjj5SvyFwEnHQih5wt8LITw8Z4laiKf0tINLGQsBFoVarQ2wBU18HShJ/DAWogKw/r+WquOhXIGmoSTxw7emgaDGrwvguKwqcSl6XgzxV2THlfhniivxzxRX4p8prsQ/U/w/woY52rHIyLQAAAAASUVORK5CYII=)
The function matchNodes
matches node numbers between different
trees based on (by default) commonality in the identity of the descendants
of each node. In the following code I will match each of the trees in the set
to the first - but this is arbitrary. Note that we must do this separately
for labels (using sister function matchLabels
) and then combine
the two matrices.
M1<-matrix(NA,Ntip(trees[[1]]),length(trees),
dimnames=list(trees[[1]]$tip.label,
paste("t[[",1:length(trees),"]]",sep="")))
M2<-matrix(NA,trees[[1]]$Nnode,length(trees),
dimnames=list(1:trees[[1]]$Nnode+Ntip(trees[[1]]),
paste("t[[",1:length(trees),"]]",sep="")))
for(i in 1:length(trees)){
M1[,i]<-matchLabels(trees[[1]],trees[[i]])[,2]
M2[,i]<-matchNodes(trees[[1]],trees[[i]])[,2]
}
M<-rbind(M1,M2)
M
## t[[1]] t[[2]] t[[3]] t[[4]] t[[5]] t[[6]] t[[7]] t[[8]] t[[9]]
## H 1 25 9 21 5 22 6 5 1
## I 2 26 8 20 6 21 7 6 2
## F 3 22 6 18 7 24 9 3 4
## E 4 23 7 17 8 23 8 2 5
## G 5 24 5 19 9 25 10 4 3
## A 6 20 3 25 4 17 2 8 7
## C 7 18 1 24 3 19 4 10 9
## B 8 19 2 23 2 18 3 9 8
## D 9 21 4 22 1 20 5 7 6
## J 10 17 10 26 10 26 1 1 10
## U 11 4 12 13 12 6 16 23 23
## V 12 3 13 11 14 4 14 24 22
## W 13 2 14 12 13 5 15 25 21
## Z 14 7 17 16 15 3 11 22 26
## Y 15 6 15 15 17 2 12 20 25
## X 16 5 16 14 16 1 13 21 24
## T 17 1 11 10 11 7 17 26 20
## O 18 15 25 4 23 8 22 18 16
## P 19 12 24 1 24 9 25 15 17
## R 20 14 23 3 25 11 24 17 18
## Q 21 13 22 2 26 10 23 16 19
## K 22 8 19 7 20 12 18 13 12
## L 23 9 18 8 19 13 19 14 13
## M 24 11 20 6 21 14 21 11 15
## N 25 10 21 5 22 15 20 12 14
## S 26 16 26 9 18 16 26 19 11
## 27 27 27 27 27 27 27 27 27 27
## 28 28 43 28 43 28 43 28 28 28
## 29 29 44 29 44 29 44 29 29 29
## 30 30 48 33 45 33 48 33 30 30
## 31 31 51 36 48 34 49 34 33 31
## 32 32 49 34 46 35 50 35 31 32
## 33 33 50 35 47 36 51 36 32 33
## 34 34 45 30 49 30 45 30 34 34
## 35 35 46 31 50 31 46 31 35 35
## 36 36 47 32 51 32 47 32 36 36
## 37 37 28 37 28 37 28 37 37 37
## 38 38 29 38 37 38 29 38 46 46
## 39 39 30 39 38 39 30 39 47 47
## 40 40 31 40 39 40 33 42 50 48
## 41 41 32 41 40 41 34 43 51 49
## 42 42 33 42 41 42 31 40 48 50
## 43 43 34 43 42 43 32 41 49 51
## 44 44 35 44 29 44 35 44 38 38
## 45 45 36 45 30 45 36 45 39 39
## 46 46 40 49 31 49 37 49 43 43
## 47 47 41 50 32 50 38 50 44 44
## 48 48 42 51 33 51 39 51 45 45
## 49 49 37 46 34 46 40 46 40 40
## 50 50 38 47 36 47 41 47 42 41
## 51 51 39 48 35 48 42 48 41 42
The way to interpret this is that the i,jth element of
M
matches the ith taxon label or node number in the
row names to the nodes in the column jth tree.
Next we want to get the edge lengths of all the corresponding edges in all
the trees. This should be fairly easy now. Note that we have one few edge
length in each tree than the number of rows in M
because
(generally) the root node does not have an edge length.
M<-M[-Ntip(trees[[1]])-1,] ## trim root node from M
E<-matrix(NA,nrow(M),ncol(M),dimnames=dimnames(M))
for(i in 1:ncol(M)) for(j in 1:nrow(M))
E[j,i]<-trees[[i]]$edge.length[which(trees[[i]]$edge[,2]==M[j,i])]
print(E,digits=3)
## t[[1]] t[[2]] t[[3]] t[[4]] t[[5]] t[[6]] t[[7]] t[[8]] t[[9]]
## H 41.507 38.43 40.89 40.93 35.380 36.741 40.698 36.70 39.46
## I 41.507 38.43 40.89 40.93 35.380 36.741 40.698 36.70 39.46
## F 25.372 25.51 22.98 24.74 23.830 27.529 25.377 26.37 25.94
## E 25.372 25.51 22.98 24.74 23.830 27.529 25.377 26.37 25.94
## G 34.043 31.93 31.06 32.62 29.796 33.251 33.857 31.85 32.47
## A 28.602 33.14 29.13 29.46 28.879 25.782 32.828 26.38 26.89
## C 3.499 3.74 2.98 1.27 1.862 0.522 1.192 2.49 5.84
## B 3.499 3.74 2.98 1.27 1.862 0.522 1.192 2.49 5.84
## D 36.666 38.67 36.18 37.60 39.167 37.088 35.332 36.74 37.25
## J 97.374 94.51 92.68 90.95 93.051 97.177 91.018 94.77 91.71
## U 12.275 16.01 15.30 16.34 14.012 13.259 15.002 13.60 11.46
## V 5.248 4.90 4.79 3.60 3.848 4.543 9.768 6.24 4.03
## W 5.248 4.90 4.79 3.60 3.848 4.543 9.768 6.24 4.03
## Z 7.952 12.70 13.57 13.49 8.555 13.479 10.517 12.95 8.01
## Y 5.860 12.70 13.57 10.00 7.919 13.479 10.517 10.32 8.01
## X 5.860 12.70 13.57 10.00 7.919 13.479 10.517 10.32 8.01
## T 77.345 78.97 78.68 79.36 77.708 80.346 80.827 81.22 78.95
## O 33.773 30.63 32.62 32.31 29.050 29.524 31.983 35.08 35.77
## P 29.987 25.67 30.19 25.47 28.086 28.537 25.924 29.73 26.43
## R 18.898 15.71 16.82 15.43 16.459 13.036 12.350 15.81 15.68
## Q 18.898 15.71 16.82 15.43 16.459 13.036 12.350 15.81 15.68
## K 21.646 19.60 18.16 19.51 22.487 18.644 18.263 19.11 19.64
## L 21.646 19.60 18.16 19.51 22.487 18.644 18.263 19.11 19.64
## M 29.576 24.77 29.13 25.65 27.042 24.809 28.559 22.98 24.85
## N 29.576 24.77 29.13 25.65 27.042 24.809 28.559 22.98 24.85
## S 66.926 60.54 69.51 66.75 64.339 70.155 69.758 62.07 62.50
## 28 2.466 5.65 9.39 2.68 3.798 5.466 12.050 5.17 7.41
## 29 6.992 8.90 7.17 4.12 8.859 6.314 5.507 7.45 7.40
## 30 19.058 18.65 18.28 18.35 22.574 18.007 16.631 18.48 16.88
## 31 29.817 28.53 26.34 27.56 26.239 36.115 28.183 32.14 27.96
## 32 37.280 35.04 36.17 35.86 31.822 39.605 35.023 36.99 34.95
## 33 8.672 6.42 8.08 7.88 5.966 5.722 8.480 5.47 6.53
## 34 53.716 46.94 49.33 49.24 45.026 53.775 50.179 50.58 47.05
## 35 8.064 5.53 7.05 8.13 10.288 11.306 2.504 10.35 10.36
## 36 25.102 29.41 26.15 28.20 27.016 25.259 31.636 23.89 21.05
## 37 19.982 21.20 17.92 14.27 19.142 20.337 17.424 18.72 20.17
## 38 2.513 0.00 5.46 0.00 0.000 1.960 4.817 0.00 0.00
## 39 64.567 56.63 60.85 59.94 60.023 61.396 65.656 61.40 64.04
## 40 0.502 6.32 2.54 3.07 3.672 5.691 0.169 6.22 3.46
## 41 7.027 11.11 10.51 12.74 10.165 8.716 5.234 7.36 7.43
## 42 4.825 9.63 4.27 5.93 9.130 5.471 4.654 6.87 6.90
## 43 2.093 0.00 0.00 3.49 0.636 0.000 0.000 2.64 0.00
## 44 12.931 18.43 14.63 12.61 13.368 12.151 15.886 19.15 16.45
## 45 24.173 21.19 26.63 23.24 26.666 27.403 22.921 20.84 20.47
## 46 8.980 8.72 10.26 11.20 8.623 13.228 14.854 6.15 6.27
## 47 3.786 4.96 2.43 6.84 0.964 0.987 6.059 5.36 9.34
## 48 11.089 9.96 13.37 10.04 11.627 15.501 13.574 13.92 10.75
## 49 6.915 8.75 11.83 12.77 5.908 8.559 13.627 10.35 10.51
## 50 14.191 11.00 12.89 11.23 9.278 15.549 14.946 11.77 11.88
## 51 6.261 5.83 1.92 5.09 4.723 9.384 4.651 7.91 6.67
Neat.
Of course, if we want we could get the average edge lengths for corresponding
edges across the set of trees as follows:
edge.length<-rowMeans(E)
ii<-sapply(trees[[1]]$edge[,2],function(x,y) which(y==x),y=M[,1])
tree<-trees[[1]]
tree$edge.length<-edge.length[ii]
plotTree(tree)
![plot of chunk unnamed-chunk-4](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAfgAAAH4CAMAAACR9g9NAAAAkFBMVEUAAAAAADoAAGYAOpAAZrY6AAA6ADo6AGY6OgA6OmY6OpA6ZpA6ZrY6kNtmAABmADpmAGZmOjpmOpBmZjpmZmZmZrZmtv+QOgCQOjqQOmaQZpCQkGaQtpCQ27aQ2/+2ZgC2Zjq2tma225C2/9u2///bkDrbkGbb/7bb/9vb////tmb/tpD/25D//7b//9v///9UN9/YAAAACXBIWXMAAAsSAAALEgHS3X78AAAOHUlEQVR4nO2dYXvbthGAFcdOtqx22i6zk62z0jbqYo0S//+/m0k5ycw65hEEcHe4933y5ENE4XD3hqRIgMSmh5BstDsAOiA+KIgPCuKDgvigID4oiA8K4oOC+KAgPiiIDwrig4L4oCA+KIhvisPVZnP2SbIl4lvicHXd97uXd4JNEd8S3ev7vf3w9lawqRvxG5ij7483or19rGdRWxnRrqoDhjLtN5sXkh3ek3jtHnihe9XYoV67Bw7Ynw9/b68Fm7opJ+IFjL/qm/txp90DDwzX8Zzj4RnclBPxeXFTTsTnxU05EZ8XN+VEfF7clBPxeXFTTsTnxU05EZ8XN+VEfF7clBPxeXFTTsQL6C6G8ZlxOsYcbsqJeAHdxTDhDvHh6F7/eo54RbRm4Nwr314jXg9F8Yef7xCvhlZnB+W7S8SroSn++P4XxGuhKf7hp/0cbsqJeAGng/wO8Vo46Kz9Hj7goJbfcNBZ+z18wEEtv+Ggs/Z7+IDFWj57TW0c+z18wGItEV8Bi7X8bp8sdnaC/R4+YLGWiK+AxVoivgIWa4n4ClisJeIrYLGW5sR3F8MD8ltu2RbGnvhXf7vrDz8xOlcYe+Jf//22794hvjDPiNdgnHN33f/BeHxpDIr/7Yfjh98QXxhzfepe//7P//yDqVelMdene+Uf/8Wcu+KY69O98v2LW8SXxlyf7pWf/sxvaq3r38VckXubfZJSp+uZfrRaw2KfpCB+BRb7JKWWeAtNCOP4/s8oBfF/ioP4jFE8iS+zrTUQvyIO4mejIN4aiF8RB/GzURBfBdE9uxOIXxEH8bNREF8FxNeJk9inZbcKltxUQHydOIifjYL49V8TgPg6cRA/GwXx678mAPF14tgTfzGc689FnSjVh8dREL/+a3lB/Io4iJ+Ngvj1X8sL4lfEWXHFrQ7iV8RB/GyUNsXbDTEP4uvHQXzlJqzEQXzlJqzEQXzlJqzEQXzlJqzEKXyv/ngjuWeL+Ppxyoo/3lyLOlGqD4+jIL5KiHvxQu+IV4iTfN9n9rZQ9/p3oXfEK8QpKP7VX2SDsohXi1OE7uLycHUp2hTxSnGKMPy4617dSjZFvFKcIoyXc7uXd4JNEa8UpwgP1/ES84hXiqMN4pXiaIN4pTjaIF4pjjaIV4qjDeKV4miDeKU42iBeKY42iFeKo40j8bXIkK8eu/sEWpuIgfh5ti/v+sOVxLwb8bVw1NU/c3o+ft/UvfpaqHQ116FqP87COLwVDMwifoJr8btR/PE94pejJD5PO+zx6bgWP5zj/3vbveEcvxzX4ofZN4cr0fUc4if4Fj9cx5/9uhHMt0T8BOfi5SHrREG8taCIn4D4rFEQby0o4ies7Gr6zZfKIH4C4rNGCSReIWhSyDpREF8yaFLIOlEQXzJoUsg6URBfMmhSyDpREF8y6FdO76s/k6xWgPgJvsWPM3BEz0kjfkID4kULlCB+gpL4HLDHr8G3+PEc39YbMWqhIj4Xwx6/f9HWO3Bq4V58v+NXfQoLuvq9A64eLb7Lthb+xQ9vu5vfFPETFolf9XVdED8B8VmjeCkH4jNH8VIOxGeO4qUciM8cxUs5QovPc//Q0EXOIhCP+JRN3WRap+9uyoH4zFG8lMO5+NPonGiJCsRP8C1+uGUrei8C4qf4Fy9bgAzxE1aL1+NBvGxRGsRP8C1+PMenvuAQ8fk3rYRomuUJxE9YuI/ZAvHpID4n9kqUA3tZIb4KrrNCfDqus0J8Oq6zQnw6rrNCfDqus0J8Oq6zQnw6rrNCfDqus0J8Oq6zQnw6rrNCfDoWs1qz7hzihRjMajeuOyeZdIf4dOxldfhxzcuPEC9k0QhuBb4sNSiadIf4dLRFT/my/Fi/RXwwTuLZ48PBOT4qwxuv+FUfkf2G63h4DsQHBfFBQXxQEB8UxAcF8UFBfFAQHxTEBwXxDXG8GYdnU9ekQbxjxA9KI74tbIrXnaDSPEOJER+QocRWxWdvGR6D+KAgPiiIDwri4XkQHxTEBwXxQUF8UBAfFMQHBfFBQXxQEN8Swytw+q29decQX5rtdd+9kcy8QnxbdG8+v5csPoZ40yRMxNi9Ex3oEW+alBk4F6LXIiDeNMvrdXz/i+wUj3jLLK/X/lz4ox7xlllcr+GtV8KpGIg3zNJ6HW+GvX1n70kaxC+jZL0QbxjEBwXxAfl2aV6o/adC5o+C+KUgPihDrRAfEMQHBfFBQXxQWhG/Kf0rtTVSxB9vhqVJtuei9p8OmRnELyZpj+9e3doapCnecnukHeq3xoZli7fcHifxS2fgHN7amohRvOX2SBPf72Qr0iDeLIm1MvZGjOIttwfiG2b+yL0YxHsA8UF5rh6la4V4RRAfFMQHBfFBQXxQEB8UxAcF8UFBfFCyi9+Pt/wMvfyoeMs+KbLHD7NwJLEzhpwPhvhHlBB/elRaEDtfSEEwxD/iefGLGb8nezoe8aoUEL8/Y3TOPvnrcbgSzrxCvCb56yGbUz/GfuKfEF+J7PXYC0/wPeJVyV2Ph2XEDT1JU7xln2jWA/GKID4oiA8K4oOC+KAgPiiID0og8fCYUpUWuHjinxBfjVKVFrh44p84IFcie6W318P0q+txpYq52BW6A98he6V3lw/rUO1mp+EgXpHsle7e3I3LER3nF59DvCJplX7mB8Phx0+Hnz9/+DSsTTPXSq7uwHKyi7/f07sf7s/093/NtpKrO7CcVPHf/2x3vbvu95fzp3jEa5Jf/P78423f/fBxfmo94hXJL/7w05u7+wP+X+en2iJekfzi5a8xRrwi+cUvaKVUwzAP4oOC+KBMKl11aAfxiiA+KH8Sn/Kt1NilGoZ5EB8UxAcF8UHJLf54Mzwd34lWpUG8Itn3+HFUTjA01yNeleziu3GIxt5br+AxqeKfm4ghPNIjXpPs4ofDvOxIj3hN0sQ/R/fms+xIj3hN8os/fvjN4kqT8Jj84k+z6kWxZ7sDxSgw+iJ8ky3iVSky7CaN/VR3KsaHbyA+KIgPCuKDgvigID4oiA8K4oOC+KAgPijrC9+9/vT177lgJeJDEogPCuKDIn6E6rszcBDvEsRDIogPyuHtMA1DNM8W8U0xvP1Gtqww4ttiu7G4fjzYAfFBQXxQEB8UxAcF8UFBfFAQHxTEBwXxjSEaoekR3xyID0oZ8WunCkBJEB+UouJFjYIeiA8K4gMyzLzavyzx1ivE22a/2bwo8vIjxDcD4oOC+KAgPiiIDwrig4L4oCA+KIgPCuKDgviW6C6G4dlryaaIb4lxaO5wJTGP+JZ4eCNG/hcjhBWvPblGxBfxp/dizGX0VJLP5J9eO9doOxXxRbxorUnEi3CSOHt8bpwkfhIvmoSDeBFOEudXfW6cJM51fG7aSxzxItpLHPEi2ksc8SLaSxzxItpLHPEi2ksc8SLaSxzxItpLHPEi2ksc8SLaSxzxIpwkzi3b3DhJfByk6V5JnpRGvAgniY/iRfMwEC/DauJPzcARTblDvAyriU/Fj+d40ctQsoqvM7VMB0EtFZj063Sov6k+EUNbTkkEtVTgKfH9VkG8IKJLrGb2lHjRXMunxaf+x7danvVYzWwqfuV1POKnWM0svV/Lvoh4YyC+MFYzQ3xhJJnlu4JYQnJGy7ZG/HPbIL45ZOLL9yNfTMSLQHyhblgH8YW6YR3EF+qGddyI32/SR+eWxFn0sWO8iN+dfer7beLz8QviLPvYMU7En8ZnROOyiBdhV/yjS3rpgjQ94oV4EX8u/+bCOGs+doxV8RNk0+1GEC/CiXjO8blxIp5f9bkpO2CSEa7j8+JGvBjEZ8NX9ojPhq/sEZ8NX9kjPhu+skd8Nnxlj/hs+Moe8dnwlT3is+Ere8Rnw1f2iM+Giexlj8r2iM+Iiez37y5lGyI+GxayP77/98+yIXnEZ0M2kFOQfpyJIXofBuIzou19KP7uUjr9CvEtcby59z/MxZgH8S0xzrmTHesR3xK7wbnsWI/4hjh+EK8wivioID4oiA8K4oOC+KAsFT9/8whcgPigZFWFeD8gPiiIb4vhockCt2znGkO8MvuzT/3xRjIJB/EtcXolgui9GIhPZ/3Uibx8WZNG1vm8lcjZmnm0RU8p+A6c2UrkbM089tId9vjuYvNCMC6L+HTspXs6x4sG5BGfjsF098POvmePL4vFdIcFyESzLRGfjut0EZ+O63QRn47rdBGfjut0EZ+O63QRn47rdBGfjut0EZ+O63Rti9ce9pgjc7o1QfwaMqdbE+viMzeYFXu9G27YbkSDc4hfgb3eNTMRw15p/x97vUN8FWYfKKv+iwPxVTAofjzHF1iTZq4xxMs/LgF7fBUQ/7UxxMs/LgHiq2BQ/HiO30jed4b4dOyJXwDi00H818YQL/9YGcSnI7iwtgvi00H818ZiiZ/BducRXwzbnUd8MWx3HvHFsN15xBdD7c6daB0qxBeDe/WGGqwJ4g01WJNVkypSQLwNNMQP5/gSy4/NppqztRINNo7eHl/gfzHIQXxQtMRnB/HLQDzMYbuyiC+G7coivhi2K4v4YtiuLOKLYbuyiC+G7coivhi2K4v4YtiuLOKLYbuyiC+G7coifhncsg0K4oOC+HwdcATic3bAEYhvqAPLQHwzHVgG4pvpwDIQ30wHloH4ZjpQDNuJqdddvQPFsJ2Yet3VO1AM24lpXxlvEK+DtvUN4oOC+KAgPiiIDwrig+JM/OnOnej+na/EqoP4oCA+KNp3ERaB+Hxou1zE16VJzhAfDPb4oCA+KIgPCuJhDsQHBfFBQXxQEB8UxAcF8UFBfFAQHxTEBwXxQUF8UBAfFMQHBfFBQXxQ/gde0ggLcFRAwgAAAABJRU5ErkJggg==)
We can see that when we average the edge lengths from a set of topologically
identical ultrametric trees, the resulting tree is also ultrametric.
And, to see that this is indeed the average, let's overlay it on our plot
made with densityTree
, this time using a square phylogram
plotting style:
densityTree(trees,nodes="centered",compute.consensus=FALSE)
par(fg="transparent")
plotTree(tree,nodes="centered",mar=par()$mar,
direction="leftwards",xlim=get("last_plot.phylo",
envir=.PlotPhyloEnv)$x.lim,add=TRUE,
color=make.transparent("darkgrey",0.5),lwd=6,ftype="i")
![plot of chunk unnamed-chunk-5](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAfgAAAH4CAMAAACR9g9NAAABO1BMVEUAAAAAADoAAGYAOpAAZmYAZpAAZrY6AAA6ADo6AGY6Ojo6OmY6OpA6ZmY6ZrY6kJA6kLY6kNtkZP9mAABmADpmAGZmOgBmOpBmZgBmZjpmZmZmZpBmZrZmkJBmkNtmtv9ra9RtbdRwcNRwcP90dNR4eNR8fNR+fv+BgdSHh9SNjdSOjv+Pj76QOgCQOjqQOmaQZgCQZjqQkDqQtpCQ29uQ2/+Rkb6Tk76UlNSVlb6YmL6bm76cnNSfn76goP+jo76kpK6lpdSmpq6np76oqLOqqrOrq66rq76srK6trbOurq6urtSwsLOxsb6zs7O0tP+2ZgC2Zjq2Zma2kDq2tma225C2/7a2//+3t765udS+vr7GxtTKyv/U1NTbkDrbtmbb25Db/7bb///j4///tmb/25D//7b//9v///8eDW4xAAAACXBIWXMAAAsSAAALEgHS3X78AAAaa0lEQVR4nO2dj3/iRnqHvdl2N3uXS7vZpOmvu73rxthIcdJfvrUBobq0pVdK5Ny1XbBk2ixg6///CzrvjCTABswMM9KI+T6fTyB412LQsxo0877zzlEKnOSo6gaAaoB4R4F4R4F4R4F4R4F4R4F4R4F4R4F4R4F4R4F4R4F4R4F4R4F4R4F4R4F4R4F4R4F4R4F4R4F4R4F4R4F4R4F4R4F4R4F4R4F4R4H42nP35ujokw+yv1VH8bdgQXp/9jZNJ88/Sp7Eeoq/ubmp+oTbQnr3Fbva7748lzyJdRA/nSZRQdf3OyHR6VwIPK/JOXUQJj69fiYrnaiD+CSJ/G4u3n93kov3PC6elJP4RqNqDeVD4tPZ0dFb6ZNaB/FTdslPixfJ9IZzexvHvLOLY/E/2UvH4Cdl/lLaPMTXnPkruqGfvJA9qXUULz4x/Tz/46l4rqh9lXJ/9ppGdE5c8RC/DA3jHfmOh3gNQLyjQLyjQLyjQLyjQLyjQLyjQHztueYh2R8lIzUQrxszM7MbSe+//wvmnAdnZYB43ZScLZDOP/vD65Q9SmZiQLwWWAuznIHuVY8RdjrNMrIE2GmYvaA4zey1ZIshXgukvctzBvwOiQ8uOsfHTfNZAuw0TF5T8tW17Gw9xGthOi1aMaKY8Si+HQx4M00Hi+mubvbi/rvDTL2C+I3QXd3dV//z9SEmW9ZLvEgKmaZRlJbRKn5XN/nl57K/B/FaqE78jHJv5i9l7+0gXpmVDjeOiwQw0dWXJl4ViFdlZbwO8SawR7zo0ZOERm+Dfr/PxusXfKROg7Zmk2f5tyn3OzjpvnuX/WXjrVIC4mVIEnqzKEqSbpfEh+0L7/ibBsSbwR7xizea8iTv0TCOB4MYXb0ZIN4AEC/XkvxxylsRj6njn1Y5nFMF4uVakj/aI/7uDY3h788kF0pDvFxL8sfN4sueq5//4q9TWjgpOYUD8XItyR/tET/7o68/pvff/7lkeA7i5VqSPz4QvxSWDcsOy/7p9x/S2Z/IVkaAeLmW5I/2iL9+++N5+vv/foXUK6MtyR8fiF/iNlvEXRL3353/eD7/1Uy2CA7Ey7Ukf9wivtycu7uvPkzYNT9B6pXZluSPm8WXzPyzj7NffI7UK9MtyR+tET97kc4/PUfqlQa29KuLmitxLKZsKxevCsQ/YsuXNMSXilnx02WSJPE7YRAEbU+wmhzfaIhX9POA/lbrxH93YmngdTsQL2LsGVGUfNEM2+2g1RTiV0voQXypmL7iV19E/ZvRiOZcH3Tu6OpLB+INAPGPxA/IbDxd/tPpo1fWDOdUgfjai5eeu+FAfHnit0wQqCO/Ml4A8WWKNzCLL78yXgDxxsXnqfjJeDgMw+DC0xuWnUnXL+ZAvHHxNFFAqfgkvtVqe81vjnWKl47LCSC+hCs++0XW0ff7ozgejHV29Wr3dhBfd/GK93YQX6J49guDQcx+Q+ccr+K9HcTXXfzsiJD/nod4VfErS6i2IWZ4b/lgjnX1msWrAvEQby32iF9eJt1fXia9DYrlUiTXa7fb79+Hfvekq9RQvUA8xFuLPeKX2Tl/Hl29KpaK33XmHeJVsVO81BuYG86pAvG1F09ztvOXsqGaisTv0kUu9ZXxuBgysf8VW4zWT7yZuXpaSzH/VHpD6crES4SmIX6b+LuvPqjM15ckvshaFzXiogGtJe5ddToXOxR1bzROG82LfLf45gX/3V636/uiQnw3isTRFZu2+kK7+CIsOzQSlp09/0klTlOS+MUAOKExcDQIaRn5FfN43HhyETnEbxM/+WPpahhEaVf84pGeKIOZ9fVx1u9trxdS866+wMjK+fT62blKhA7iSxRvIueOr5NVSMaoSLxQTc742Gb7ia35cM4s/GpXuOQhvubiRa6l/CVfR/E32Xdl9gVRFBNVKyy28mvsxUAsoaqLeFUgHuJNspd4NvL3u1GG/+4k5HQ6nscHeDQNsEdhMbH0udglzmu874Vhv9fN5xz2GCnaDMRDvEn2Ev8QveNhdPUm0Ste63gY4k2iVbyJpi1e1Gw4pwrE11/85Ojo6CDH8WU0bfGiRPE6vqnS9Pr5x3T+0tLKlhC/Dh03K6koXyy9WBriyxZflNSLokHQZhQF9RTGosVqWenZeogvW3xeVy/yu4PL37Zy8aenKhMRWXQO4pWbtnhh/orPnpPpTb83yvt6ek+VOec7sTXFBF29WtMWL2omnrr6+Wc/vcHNnVrTFi9KFH87uIrpvjz7sdq7zJ6dswGddCEciK+7eDaUwzhevWmLF+rid+ubi/Ko41h09cUb7XEW5BOsIV6n+F2G5UbEywPx+4unITmleV+JcDGlgO8QBj5tNIPL973hsIj/lhr8hfj9xSddni7gd4KAxHtes3G8eUgO8QcjnhmjP5+Oipzx8XjzyAxdPcRD/GGIj4sPlSSbf0H7cE4NiK+9+GulSncQryp+qfMWPXssunoJ8VqmbEWUZm7ppsIHKT4ftVcrXszd3Fm6jfiBiM/X+Ed8iT/tT8cXbp82PRq9t4OAlm/7PhvebR6ZLYVl+xrCslkKxsTS3aQPRDyTxle605L8QdhutUk8T+u/oDp2lYif8K94We8QL3nFi/+lh5vRcDi6ET096+sZu3X1S+hYIMAzcKSvd4ivVryGnDt+b3d/huicYtMWL3YVL+7pUj6aE933bsM5zYicK/lLHuJrLl7c290hA0exaYsXT2wxyv+XHkRXn1YtXhWIh3iTWCz+0Tbira3biPMC5Gz4dhGw0VzYpzi832XDO3YYiH8MxNsGxD9iy+gaXb0stRK/eXQN8bLUSfwWLBzOZWFZjOONsl38bhk4mkkVp+4gXgYbxYupO+kKhxAvw3rxEQXkul0/3CnLVit5WHb2CRZUmMRG8TwsK+0d4qVYL37KMzSmZsqSP0W2Pl4aiJfhCfEmypI/heXbiDshvhIs30Yc4k0hXfUoQ0b8Hh3SSpCLPd0UySoQXw1y4pW/wyDeNp4UT/OPWeHoK8onzpJCeX2u5arP2xFLRHlki8e2gnabjXz8kxP/nT+tT3Vol8ST9q4oGt6hfOI8XpkFLDOeygyGeNvY4YovPhQPR+WpKY+qPqOrrxUQL4Oj4vknzdNPV2U99bExnLMNiJfBTvFmypZD/BLTtStp4sVOWOVjrGw5xC9ho3hTZcshfonl1bL9YrUsG+DyLa8rwGDZcohfwkLxxsqWQ/wGKom+P8ZY2XKI30Al0ffHiDRLA4smId5yqGz5/Zn+suUQbzszhWJnEO8sEO8oEO8oEO8oEO8oEO8oEO8oEF93RElT6YB8teLjOB6Mb5/O3AIboRMqvUY6rV78qD+8hfg9SFU2nUurE0/1XwfRVS+8vGx7XhURzYOAi1daPleV+AHjqnvVC779tuV5jUY14ezaw8UrLZ+r7opn3+/T0WgYhKNb/qrqPrOmsPM5UfiKh/i6k6ZqpRGqFR/H44tOLKpLyrcdcNRKI0B87VErjQDxjgLxjgLxjgLxjgLxjgLxjgLxjgLxdUdsOXctW80W4m1DdsqWr6EysOEgxJeM5Jq8lBbIy1cth3g7SIpigovtyWkB9lNlBHmQ5uyXn8qHaSDeCrrdJCsmuNie/OJCVAfcVkZQxOOfKYTnIN4KkqQ4hUWJFYpWZ/35lqB1ms5f/b38mkmIt4M9xM8/PVcJzEK8FSyJL8qo8Z8KtpwdXhJBuh4GxFuCsvj5S+rm5QtiQLwdZOKp7y52sNztilcF4q2gEH9z8wPEO0QUJVRHbTwcDn/3r73ekO9ZO+12xdk3UdYd4q3g5CTyu5n4fwmF+CQ5OYH4A8f3+epRqpz3w+/62S7V7IrP/hhd/aEC8Y6Sied39cNBXg4d4g+ePcSLtPqJ7Hw9xJfPmonXTocXCqCu/j/+rZ+XQ7+62mWunlbHz6XjcxBfPmsi7vuIZ5c8Zu6sZFrAt3KjxeF5xP0i277v179uftPwvDbjn//p29++vzylQui/+c0uYdn7786vMVdvI/TpBGyw7vtMfKuVi282TjPxf3OsKD6dKWwfD/ElsPhsU/onwAZtw2He18dj3osPBkVX/8O/j8Q3ATs54409/FJXr7ZOGuLNIyn+P0fiD3cWX7MaOO6KFxkWxU/oKYqK4Vxxhneeqa1hKRSI1yF+opB5BfElYFq8EhBvHoiHeIiH+OIn9ATxhwjEQzzEQ3zxE3qC+EPErPiJUrn6qsTfDIej0SgWJU0hPhfPt6sdZTHYnadsVZbHQ3wJGBevsDy+LPEiMsm34OabcPfDsNfrD/r9XqvdF38g3fTasIP4kxP/i3e+HxJtL4vYNo+PdypbrpCGAfFlYFj8/ZnFc/UPuR0RMevX+v2ReBO5A9SJHcRnf0thT/o0VUi/SasTH5N3iH8oXn5PeiqBo9KoisRnh2V3MINBDPGrf0sK+fxaAcSbx6j4az6Mt3b/+A2HFV09jVwOeWuSxWeLadrqVqReZedhX/GqQLx5IH6JhKeY9/vh5fuW51Eu8VMjl/pSVC1rnjaantdhI7Z2vy+Km3X9SJxhoxPza4B480D8IxTGrTUEXf0jFMatNQTiHcXocE4ViDePWfFU3JBP1/P65TsD8ebZX/yW7xEufvaXL1LJuVuIN48O8RvvhWgB1f33//VCFEjYHYg3j7z46QpUBq3fC4ML71HheiF+/it6kNtpFOLNIy9+saKeiKLxsBe0W16zcdx4LP75x99/YOLlLniILwGVK37192+owO0w5gIedvXzz/7vc/bw09dy0VmIN49p8X93zh7+IJl+BfHm2V88Tf2Mx/Rt//g8z3/GvtznP/8ryXQMiDePWfE81XL+UjbfEuLN84T4RxO7D8ub8RQ11tWvF68IxJvnSfF8kA7xB8cm8VNKNT/xr3phGLS95Q3mFoFc8cN2u3XZavtdv6stdgvx5tkkPulGif/FCdU7bF9ywXlBO4g/CDZf8VMmP2H9+Gg0jHn6WdbHo6s/CCDeUZ4Qz3NNx+I0Z39VZjinCMSbx6x4EZyRXk8D8eZ5KD4vaUrlTMdXY9HVbxW/bcqWFzSVL2cL8eYxK547l69qCvHmeVi2vJeVLfeaF1dew+uxcXwv5MnWeY71w7DscEtYluKxCvWrId48hsVPXquUMYb48slXE1BXz/7Lunp+Ejecxy3rD6ibVylYD/Hlk2fQSYjfknPHbuv/V6E0AsRXBw3neJCGD+e2id/K3RvZrccIiK8OTeJrVQMHpNrEqwHx1QHxjgLxjgLxjgLxjgLxjgLxjqJJ/ISn1EsunYP4CsnEZ6lXPOFKvupbFpmT3mYU4qtDj3ihXHKRNMRXCaVXd0/8nkiv9rxFevXO8EWTPPVKNk4D8dWhR7zYWVb2Kx7iq2e/Yn9COXLuash+xf6EcumveIivPbRBhcKuNBBfd+7PVLadg3hXgXhHgXhHgXhHgXhHgXhHgXhHgfhD4MdqllDdPqzdAcqDBEgH47WJvxmNqj4BrkIC5Gfq1cXT9mF8T/Ak6l7RQt+Qdr72PE8urAj2Q4hXWCWtLL7bzcV3TzokPmh6tOW59/Se50AfQvxEYfGc+hWfd/XTZMSTh8YieYjq+oDyIAUK93YQX3fY1a5yb6dFPDceJyJBuOSdMp1n8kK+1BkB8TXn7s3RJwoXPMS7CsQ7CsQ7CsQ7CsQ7CsQ7CsTXnskRQ3rSFuJ3oerpuS2IZXPyKyogfhf2W+VklGzCtrTVsq6Ip10hkmTc74Xt9tL+YLZQLJMubX28K+JF3sG4F7TfX56eNqXXMRumWCZd2tYkrogXH542iOgP49jCHLNsZfwEXb1erBeflcD5sqz18W6Jp02gBuPp1MJ0Yv7lrlDAGuK3Y734GY3iFQrWQ/x2rBevCsRvB+JTiIf4FOJrD8RvB+JTiIf4FOJrD8Rvpwbir5UG8hC/HevEP5yyTe+//6AwVQ/xT2Ch+NXUADFnW9424hBfdiuIJIquaEl6cJGlBjDx3Pk15uo1Y4f4fFvxyO92wiAI216jwVMDKB6vlHIH8U9gh/jFuZ6OGDcj1p4sQizCsrjidWO9eJFyh5s73Vgnnp1pksBOs/ghv7crsWw5xJfbiuz5sXgejy+xbPmu4ktKQDIGP8+3IvVqi3jjrciex7Ho6hfiVTEv3tqM9J3YWbzRj1kL8dE0WWLQ7zGuPK+mldBozOR57dZlqxVEEWXZ5+clefAx2TDrQpR80599zwZvnNNGs90Ogl7o+/50vy9V/eK7RQU8guYbmPhGg4mvYyW0jeJp77CVj9lutbl40q47+74W4qPpch/Ee6ab28GAdVR1LIi1saunbSJXPuZwOCp+Q3cSdk26eoiH+FSchDSK+O/s09Jq2DicWxVP+0GP4+I3dA/6tgznVIH47dRB/Pzl0ZEFM3eHKH7NEqp4vNydi66++A2D4h9O2abp5JMP6f1Z9WvnIL5c8fNXfBvxVyXtH++KeLFMevh4mXTj+NRb0Gqx0VwYEaKut+5WFGHZ8FFYVqF0dQrxT2G9+Psz+Xl6AuJ3Yc1O36V29VuaAvEmWTMVX5n4hzl3IgdDusOHeEVKHc5tY0ap1RPpyuUQr4g14nlEXn5TGohXxB7xakC8IhAP8SnEQzzEZ0C8pUC8IhAP8SnEQ3wl4mcKG5BBvCKlpl5tIaV6pj+zocAhxJctfvL62oaVNI6IX02v7plNr94Iib/76oNCTB7iFbFH/OStygbyEK+FNQH7srr6u68/yu9PAfGaqG6JIC96pbCvMMTXHFGv/jsLyp1BfJlkmVfyt/Uy4mnEcjMSIxX2b6BY1jO6uWHDGYivFRDvKE+Kp/xiyhhm0IglaDWzte6edyHWgzeDXi8Ignc+/cWk60N8HYB4R3lS/BLLg9XFkmfe1bPOHl19rZASvzRYhfiaIyN+mYVEPpxj/iG+VkB8/ZmoFLrTIV509WOxiHg8WAlKilIohyjewPSrCml6/fyjif3jNwDxllTwy/aR1r9//AbYwC2jexJSKbug6VFI0jv+hhcIW1o/3A5D/4sv6G92Fd/MNqY0xh2EYdBue81Kd5UvdpOWBuIVSLpRtztot9uXLY8WzVe3q7woW05FcLSvlt2B1WB0fPhd/ZTqFNwMh6PhSBQkqXBz8WwDcenrXov4lW87iC8VEZ4TVetl0CH+wVk5/OEcF09xqnEs6hJVmUdPRa/4nb0cEK+AVeL5MN5Azp00EF8LIF4BiF8HxNcCiFcA4tcB8bUA4hWA+HVAfC2AeAWsEn9/prTzXH3EVzYp+hgmPI7FlO168SW2RWG2llMn8VbEv4kdxJfWVpUFk4S94sVKjqUdvkLaOF3Evx9RbjC0cXzabAatVhj0EtqTbGk3Ot7maNALKFbvGW9Zvn+8PPaKFzXaF3vYXbSCMGyfNk7XiS83IL5ZPG905F+Fl61Wy6NgvdmW5atl5TeQt1e8OMriWMM+LdUaj9dGQMuNi27u6sWrZDTqD0cj3izTLVNYKMuBeAVsEq94b1cj8eMBZXMmydqRU7njqc3DuUx8HA/Yvw3+0nTLFO/tIF4Fm8SL/eOf6S+MII0r4temXolXY9HVlyJeFYhXAOLX4YB4nl7dX5dezff7Pm202pftXo+y7/fc7NscEK8AxK/DAfGctaXt0NXnHLD4dfPxEJ9zuOLXUv5wThWI10oV4vlkvf7946VZFS+WUN0OBrfL5VN2QQyRFjOeYsq2HuJLnLK1NR4P8abFWxOPFxtvL9XzDnu9q0bD87zjY5mIIx8ZNRdLri/fB+xgJ340XYMtg6YsLNsrMSxrTTwe4ssVb6p69Z6oF3KvZVefUWL5envi8auoZ5/VWnx5OXf2xON1UcvhXPnMP1f7PYivOQr7EHEg3lEg3lEg3lEg3lEg3lEg3lEg3lEgvvaoRWkgfivlTLvug+oUDsRvxaJF+RtQ2V6UgPgNb550iQ6tyQ+y7cErrUy/HpsycHRRqfgkid6dMPEXtCZf7Asv8gMqq0y/HpsycHRR8RU/jSi3Y8S3T+V9qkgcq7BA+XrsycDRBcTvRDopd2sS81ghPuY7LPEfiewu68aNtmbgqAPxO3H4GTillkKJaUnELe2sQl09f39LxasC8RBvGQ/LnQVlljtjxzs+7XQ6AaV0h3wTdd+3KY97byAe4i2nxFx1dPU2Ueq0OcQ7Sn2Gc6pA/FpqJF7M4EjP40D8Wmo0ZStmcKTncSB+LTUSL6Jz0jE6iF9LHpYN7A/LiuicdIwO4tdSI/GTI0v2pDkoSp08UEJsIy4fo4P4rdifc6d4bwfxdSe7t5NeJQ/xNUfc1cmnWEO8o0C8o0C8o0C8o0C8o0C8o0C8o0C8o0C8o0C8o8iJPwJ2UZp49TfC8ew6HsQ7ejyId/R4EO/o8SDe0eNBvKPHg3hHjwfxjh4PM3eOAvGOAvGOAvGOAvGOAvGOAvGOAvGOAvGOspv4+asPaXr35uj5x+JpT+Yv+Wb3eg7GuD87enau8Xh0xLcaj8c+75HO46X7Hmsn8TOSRCdi8iJ/2pO7L8/TyfOPeg5GXL9NZzqPl6YTJkrb8ejzzj8919m+PY+1i/jrZ//Arnhae8+u/OxJ+Q0FtKqbHUnPwdKiMIC247EW/vzP3uo7Hl/MfP1WY/v2bdvuXT139eV59qT6fhnZFa/nYNTAz/6Wunptx0vvv/tHdkXpO16a6jt5gj2Ptbt41pPSG2VPqu+XI76gNB2MvkK5JW3HSyevqSvVdzzqml9rPd6ex6roimffd+nskw8ar3iNjRPHu9d7xd+9eb33VbpCaVe81u/47J+rvu/4r/UeT9QQe63vO5m6JK33IOV9x1NXxe/qX2u4Mc2ueD0HI67fiqtU1/HEXbO24wnv+o63/7GqGsfPjvSOu9mBtM4LaB7Hix7kbd3G8eDwgHhHgXhHgXhHgXhHgXhHgXhHgXhHgXhHgXhHgXhHgXhHgXhHgXhHgXhHgXhHgXhHgXhHgXhHgXhHgXhHgXhHgXhHgXhHgXhHgXhHgXhHgXhHgXhHgXhHgXhHgXhHgXhH+X8cVrDyFItl8AAAAABJRU5ErkJggg==)
par(fg="black")
The 'average' tree is show in grey, whereas our sample of trees are given
in blue.
Note that it is also possible to do this automatically (that is, obtaining
the average values of corresponding edges from a set of trees) using
the phytools function consensus.edges
.
I hope this addresses the user's question.
For this exercise I simulated my set of trees as follows:
tree<-pbtree(n=26,tip.label=LETTERS,scale=100)
trees<-list()
for(i in 1:9){
nodes<-sample(1:tree$Nnode+Ntip(tree),10)
trees[[i]]<-tree
trees[[i]]$edge.length<-abs(tree$edge.length+
rnorm(n=nrow(tree$edge),sd=3))
trees[[i]]<-force.ultrametric(trees[[i]])
for(j in 1:length(nodes))
trees[[i]]<-untangle(rotate(trees[[i]],nodes[j]),"read.tree")
}
class(trees)<-"multiPhylo"
although I assume that in the empirical case they would come from my
inference method.
That's it.