[med-svn] [r-cran-ape] 01/03: Imported Upstream version 4.0
Dylan Aïssi
bob.dybian-guest at moszumanska.debian.org
Sat Dec 3 21:59:29 UTC 2016
This is an automated email from the git hooks/post-receive script.
bob.dybian-guest pushed a commit to branch master
in repository r-cran-ape.
commit 3b7b29efce9e96af34b16a47b69f9a3bb0632dda
Author: Dylan Aïssi <bob.dybian at gmail.com>
Date: Sat Dec 3 22:57:06 2016 +0100
Imported Upstream version 4.0
---
DESCRIPTION | 14 +-
MD5 | 126 +++++++-------
NAMESPACE | 39 ++++-
R/CDF.birth.death.R | 6 +-
R/DNA.R | 202 ++++++++++++++++------
R/MPR.R | 2 +-
R/SDM.R | 2 +-
R/SlowinskiGuyer.R | 25 +--
R/ace.R | 6 +-
R/as.phylo.R | 2 +-
R/checkValidPhylo.R | 9 +-
R/chronoMPL.R | 2 +-
R/clustal.R | 35 +++-
R/coalescent.intervals.R | 2 +-
R/dist.topo.R | 124 +++++++------
R/is.binary.tree.R | 37 ++--
R/is.ultrametric.R | 45 +++--
R/ltt.plot.R | 6 +-
R/makeLabel.R | 6 +-
R/multi2di.R | 54 +++++-
R/node.dating.R | 277 +++++++++++++++++++++++++++++
R/pic.R | 2 +-
R/plot.phylo.R | 29 +++-
R/plot.phyloExtra.R | 32 +++-
R/read.GenBank.R | 80 ++++-----
R/read.dna.R | 16 +-
R/read.gff.R | 22 +++
R/reconstruct.R | 421 ++++++++++++++++++++++++++++++++++++++-------
R/reorder.phylo.R | 37 +++-
R/root.R | 103 ++++++++---
R/rtt.R | 16 +-
R/summary.phylo.R | 116 ++++++++++---
R/write.dna.R | 2 +-
R/yule.R | 2 +-
build/vignette.rds | Bin 192 -> 194 bytes
data/hivtree.newick.rda | Bin 2189 -> 2190 bytes
data/landplants.newick.rda | Bin 575 -> 0 bytes
data/opsin.newick.rda | Bin 445 -> 0 bytes
inst/doc/MoranI.pdf | Bin 238263 -> 238589 bytes
man/AAbin.Rd | 3 +-
man/all.equal.DNAbin.Rd | 65 +++++++
man/alview.Rd | 2 +-
man/ape-internal.Rd | 6 +
man/ape-package.Rd | 4 +
man/boot.phylo.Rd | 9 +-
man/c.phylo.Rd | 23 +--
man/checkAlignment.Rd | 2 +-
man/clustal.Rd | 31 +++-
man/dist.topo.Rd | 21 ++-
man/image.DNAbin.Rd | 5 +-
man/is.binary.tree.Rd | 40 +++--
man/is.ultrametric.Rd | 39 +++--
man/landplants.Rd | 36 ----
man/mixedFontLabel.Rd | 2 +-
man/multi2di.Rd | 25 ++-
man/node.dating.Rd | 107 ++++++++++++
man/opsin.Rd | 35 ----
man/plot.phyloExtra.Rd | 2 +-
man/plotTreeTime.Rd | 44 +++++
man/read.GenBank.Rd | 62 +++----
man/read.dna.Rd | 4 +-
man/read.gff.Rd | 46 +++++
man/reconstruct.Rd | 85 ++++++---
man/reorder.phylo.Rd | 8 +-
man/root.Rd | 64 +++++--
man/summary.phylo.Rd | 20 ++-
src/ape.c | 8 +-
src/dist_dna.c | 126 +++++++++-----
src/read_dna.c | 26 +--
69 files changed, 2121 insertions(+), 728 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index c23536e..c8452d6 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
Package: ape
-Version: 3.5
-Date: 2016-05-23
+Version: 4.0
+Date: 2016-11-29
Title: Analyses of Phylogenetics and Evolution
Authors at R: c(person("Emmanuel", "Paradis", role = c("aut", "cre", "cph"), email = "Emmanuel.Paradis at ird.fr"),
person("Simon", "Blomberg", role = c("aut", "cph")),
@@ -14,6 +14,7 @@ Authors at R: c(person("Emmanuel", "Paradis", role = c("aut", "cre", "cph"), email
person("Olivier", "Gascuel", role = c("aut", "cph")),
person("Christoph", "Heibl", role = c("aut", "cph")),
person("Anthony", "Ives", role = c("aut", "cph")),
+ person("Bradley", "Jones", role = c("aut", "cph")),
person("Daniel", "Lawson", role = c("aut", "cph")),
person("Vincent", "Lefort", role = c("aut", "cph")),
person("Pierre", "Legendre", role = c("aut", "cph")),
@@ -26,16 +27,16 @@ Authors at R: c(person("Emmanuel", "Paradis", role = c("aut", "cre", "cph"), email
person("Klaus", "Schliep", role = c("aut", "cph")),
person("Korbinian", "Strimmer", role = c("aut", "cph")),
person("Damien", "de Vienne", role = c("aut", "cph")))
-Depends: R (>= 3.0.0)
+Depends: R (>= 3.2.0)
Suggests: gee, expm
Imports: nlme, lattice, graphics, methods, stats, tools, utils,
parallel
ZipData: no
-Description: Functions for reading, writing, plotting, and manipulating phylogenetic trees, analyses of comparative data in a phylogenetic framework, ancestral character analyses, analyses of diversification and macroevolution, computing distances from allelic and nucleotide data, reading and writing nucleotide sequences as well as importing from BioConductor, and several tools such as Mantel's test, generalized skyline plots, graphical exploration of phylogenetic data (alex, trex, krono [...]
+Description: Functions for reading, writing, plotting, and manipulating phylogenetic trees, analyses of comparative data in a phylogenetic framework, ancestral character analyses, analyses of diversification and macroevolution, computing distances from DNA sequences, reading and writing nucleotide sequences as well as importing from BioConductor, and several tools such as Mantel's test, generalized skyline plots, graphical exploration of phylogenetic data (alex, trex, kronoviz), estimati [...]
License: GPL (>= 2)
URL: http://ape-package.ird.fr/
NeedsCompilation: yes
-Packaged: 2016-05-23 07:32:24 UTC; paradis
+Packaged: 2016-11-30 09:33:29 UTC; paradis
Author: Emmanuel Paradis [aut, cre, cph],
Simon Blomberg [aut, cph],
Ben Bolker [aut, cph],
@@ -48,6 +49,7 @@ Author: Emmanuel Paradis [aut, cre, cph],
Olivier Gascuel [aut, cph],
Christoph Heibl [aut, cph],
Anthony Ives [aut, cph],
+ Bradley Jones [aut, cph],
Daniel Lawson [aut, cph],
Vincent Lefort [aut, cph],
Pierre Legendre [aut, cph],
@@ -62,4 +64,4 @@ Author: Emmanuel Paradis [aut, cre, cph],
Damien de Vienne [aut, cph]
Maintainer: Emmanuel Paradis <Emmanuel.Paradis at ird.fr>
Repository: CRAN
-Date/Publication: 2016-05-24 05:38:54
+Date/Publication: 2016-12-01 17:40:32
diff --git a/MD5 b/MD5
index 3a67d27..6a8e3b2 100644
--- a/MD5
+++ b/MD5
@@ -1,24 +1,24 @@
eb723b61539feef013de476e68b5c50a *COPYING
-d13712c284641e2ad23e16dc3ea68ec6 *DESCRIPTION
-030bd4cdd235a95eeb7a9463793b1865 *NAMESPACE
+e28f0a5656abc88ffa75eb720a68e0ff *DESCRIPTION
+6936565af4821e47088568525c3bda01 *NAMESPACE
854a025cb7e5da3e4fe230c0be950d08 *NEWS
0c7bc9101516fd26fb3ddbedbe25b6a3 *R/CADM.global.R
e042efca1d8a00d4178204f5f481b64d *R/CADM.post.R
-df7599318fae6e660d57461413148967 *R/CDF.birth.death.R
+a10ebb0e9ecf861eea483eb41da1bd23 *R/CDF.birth.death.R
fdb9cfa0cbda82bda982b290693e44e3 *R/Cheverud.R
-43dc0cac925dff273611fe3530075d2a *R/DNA.R
-0fbc7715dfdc8d54ded16d78ca3100f8 *R/MPR.R
+32741db1a4135fdd540d71a103d78d6f *R/DNA.R
+2187a1289b767066d1efe1ebbe7c3b0c *R/MPR.R
bb95af56d882b6162aa517a77140175e *R/MoranI.R
218b8cf3c13a700c757b7b2303dc897e *R/PGLS.R
-a683ed4e25ad199edc7d30e7b55181a6 *R/SDM.R
-6b54088bf2514e3aa70edf05fe88127d *R/SlowinskiGuyer.R
-eaad016718d9ee8c79b15d252879072c *R/ace.R
+6be0924b9f043abaee0968de5cf62aa6 *R/SDM.R
+315ae41ee77b323b07e63ff861359ae2 *R/SlowinskiGuyer.R
+8240ac5d3ec4c439e9a33a33f66e8428 *R/ace.R
4ce79cf3f3ff49bef989454d86d0c891 *R/additive.R
9fca35b1005cce7b2a260b553dd971a3 *R/alex.R
9fe874382f024a98f62a0ccfcd6d09ac *R/all.equal.phylo.R
06b049b681ab2e594124a2b010c5cddd *R/as.bitsplits.R
c94018d5e792c72e20ce84085b2df9e7 *R/as.matching.R
-ea82e116c735fbf545d888f3aeb36fd6 *R/as.phylo.R
+86d55d04b20c1ecad73ea817866d632d *R/as.phylo.R
1bb4dcfb44a7485710ec7f599d4b2ee8 *R/as.phylo.formula.R
844a38494dd34adbe916bb9116cd0bc2 *R/balance.R
607d7d6bc0ec5559ce1de5b1d1aab0c1 *R/binaryPGLMM.R
@@ -26,13 +26,13 @@ ea82e116c735fbf545d888f3aeb36fd6 *R/as.phylo.R
a28930cdae07a271403e988e86a0c324 *R/biplot.pcoa.R
92d9db6d7c1150d1a7c9443fd3d5cb04 *R/birthdeath.R
4ee816e0dc69630d104fb5da4cc542c4 *R/branching.times.R
-cc55d092fc68c6433da930832ba03177 *R/checkValidPhylo.R
+b4557270916972166e13e76391bb0a69 *R/checkValidPhylo.R
e43b5dec7eae6d4bf9371e50117bf6ed *R/cherry.R
-02bdacc928572e8ab2d88380a075d7a8 *R/chronoMPL.R
+5fd65efbdae6e0934e8bb42191142d4b *R/chronoMPL.R
74e1019810b06458e808a447bb099a91 *R/chronopl.R
8ae1af4a30c40d15676d70e180f3d593 *R/chronos.R
-3121541e330547af3b5eb7e2ea78b6b3 *R/clustal.R
-4f52fb37d822be586f0c5971295167dc *R/coalescent.intervals.R
+2a9983307f686facf6d9ec843d214f38 *R/clustal.R
+dedf66f0595977a11a7dbc2b0d2f56bb *R/coalescent.intervals.R
c3d44d960238a711764aaeebe9570559 *R/collapse.singles.R
338accc0accb8e67f31db5a90c033d2f *R/collapsed.intervals.R
01e979242ba4667e48a80fd32c03f254 *R/compar.gee.R
@@ -46,7 +46,7 @@ fd39268020a494980293c92d0971caca *R/corphylo.R
3822f0bb0a9ed4c8c19654e86ef34359 *R/def.R
93480a5b64e0d37f5450891629557615 *R/delta.plot.R
dfd5bb35f1cb1fd9154d023e0e4cfc2b *R/dist.gene.R
-9ed23c3cf1f576f93de0dfb725289ab1 *R/dist.topo.R
+c59b82a465d826b821896e0da82f1017 *R/dist.topo.R
b28ced504fedeb7f991f7eba10ad06df *R/diversi.gof.R
8b2ec4004022afdc7e2cb42f2657b628 *R/diversi.time.R
ea5fcbd0f207ab1d874c0d032392967d *R/drop.tip.R
@@ -56,14 +56,14 @@ aa09abeb90ef891384128f978ffce843 *R/extract.popsize.R
5c29d3ee785da587f4ad5288ec36b76a *R/gammaStat.R
499b5f8596f40c32592fc85ee0e226ce *R/howmanytrees.R
68d848281455c4c91e9b91f16170e2f7 *R/identify.phylo.R
-670971f4507c1138ca5d55cee936e4bc *R/is.binary.tree.R
+7d6ba4bcc70903878b1bba0565ac4722 *R/is.binary.tree.R
be075b35bb032eaca5946abd14ded1e3 *R/is.compatible.R
35921387c705612d8f7c5baa06f9ab79 *R/is.monophyletic.R
-135e6538fd2db96b8bd39ad42d161506 *R/is.ultrametric.R
+a7bd37de10eb3f5c02096dfb9a492307 *R/is.ultrametric.R
944879f7a9a134866f31d8476a07f9c1 *R/ladderize.R
65b2494ebd918c6f9a31c80e25370035 *R/lmorigin.R
-b4c7ccf40ec0c69a21da38965747d350 *R/ltt.plot.R
-2adfbd3f101466e8b04c4f9aaf3d5377 *R/makeLabel.R
+1553fd068a844a4274c62b936e242131 *R/ltt.plot.R
+0e3a341c99497686db7583eb0894ba77 *R/makeLabel.R
34069210fd7b12dda0979c45822e4d3a *R/makeNodeLabel.R
a84000c0d7d2f53aec6a7105477f886a *R/mantel.test.R
d2c16632492bfafd2ee18f2fe3d3d64a *R/matexpo.R
@@ -71,43 +71,45 @@ d2c16632492bfafd2ee18f2fe3d3d64a *R/matexpo.R
61021f7af1175c46a743c7fee4cdc87e *R/me.R
e1327a592199bac7ac8f670c3a068d46 *R/mrca.R
a078728fb5907565f85b54b30e5bf83f *R/mst.R
-eb5dd3843951b8bde1b8d53ab2d1353a *R/multi2di.R
+93f3fe244219d565bb2e111ff249b540 *R/multi2di.R
0850fdd19c01d37ac632fc308632a463 *R/mvr.R
282308c27ac1e78e330711d635d2e572 *R/nj.R
e3f22d0f260c43be87a17b0ab091e2bb *R/njs.R
+5e43665c088ced872dd2a03f7adc2add *R/node.dating.R
e4a0712fa45754cd02f5d3da6ef23204 *R/nodelabels.R
ae2aeb0e8aef7f8d4b19939ca61b3482 *R/nodepath.R
d9fd8e402e6fce6939a05332823a9390 *R/parafit.R
fc0260f9e8e17fc1446117580bbc93cc *R/pcoa.R
431d25db557d480985ac16bfc3d6425c *R/phydataplot.R
e71db002c66c277bfb57f6914ca143d4 *R/phymltest.R
-61d1360195bf4d5e6e2d8c30de335dc8 *R/pic.R
-7feec819ed70a21c0119d137da22c7d9 *R/plot.phylo.R
-377ed3a292afcad0773924b52444f6b4 *R/plot.phyloExtra.R
+d04e04530104c16108684db049d715da *R/pic.R
+198dfc699186f842a0cce2bedf4eabc9 *R/plot.phylo.R
+b31f836442e0a6dd739d446faffda284 *R/plot.phyloExtra.R
e579ec65c808c29e1ecaae1501784b37 *R/plot.popsize.R
736506e2f67e90cf9326abead437d298 *R/plotPhyloCoor.R
1e2485437566ca9af99d93b4580cbbc2 *R/print.lmorigin.R
d0e8bd41d5acc217fdee3578adcf635b *R/print.parafit.R
8c401518738b9cda403fa9f0eb382757 *R/rTrait.R
-1895483cd4f45095860092900e0672c8 *R/read.GenBank.R
+9c261f833a2f3325949b00f49a3d9d61 *R/read.GenBank.R
b13dfb8f455b1c9e74a364085f72dbce *R/read.caic.R
-c1fc8ab4715f1d58145f533fbaf776e4 *R/read.dna.R
+740c43ccf599f0ea029fddf4361f5738 *R/read.dna.R
+b25674d229dfa4d2bf76a50e1874c690 *R/read.gff.R
12cc225d7dbe722e9e61cd3911de0b5a *R/read.nexus.R
13ce7f5c7d1bcb7101469d12651e99c8 *R/read.nexus.data.R
88e633c80886f5a6c9f2b34918452d15 *R/read.tree.R
-86361682a3d477e8e63fff53dc67279e *R/reconstruct.R
-4a39f9ba8cc833bdabe7bf32d9e40651 *R/reorder.phylo.R
-623e12a556b5d22e6615b1bb1205ce04 *R/root.R
+df2ac5d3de7185c7e26fc95b35192a40 *R/reconstruct.R
+1770cbeb0312916ebe52647fc29a83fc *R/reorder.phylo.R
+4400af9cbac48963c7cb287c51aad62c *R/root.R
f584366b32e7414c669714ba5b84951b *R/rotate.R
b7158b84c7ee7d9dcb2e0abeb6005cb0 *R/rtree.R
-6d359208c90e5aff93426f3c6ce39f46 *R/rtt.R
+be352a2ab62b61eaf8186a7f9b682798 *R/rtt.R
d099c8987470c4506511be54e29a5ddd *R/scales.R
d2e06f8288af941a00c46248b586225a *R/skyline.R
1f82059f740388b7430b2359e54a147f *R/skylineplot.R
9c7b02a4625099f715700fb868226b0f *R/speciesTree.R
28fce999b31a0c7d53816cec6916981f *R/subtreeplot.R
bcc8f1fc8363728caba82129412d9e31 *R/subtrees.R
-fbbe3d3c1eec119690a296bfaff4b252 *R/summary.phylo.R
+e9daa8dbb2eef717c51a73e92669ae31 *R/summary.phylo.R
8fbd1589f5d98d76b1154cffb8d4d1f5 *R/treePop.R
b5081fca8758fe4458183c3e25e3e661 *R/triangMtd.R
6e92716e8004feb088d5c093bad3828f *R/unique.multiPhylo.R
@@ -115,15 +117,15 @@ b5081fca8758fe4458183c3e25e3e661 *R/triangMtd.R
a40ae9ad30c221d4ed14b90e8b406f93 *R/vcv.phylo.R
31b3bb1feed474692f07fcebe3a61ac7 *R/vcv2phylo.R
16d5511c871e41b9b2673052e03c8b33 *R/which.edge.R
-aa484437bf64a39b737665b0769e320b *R/write.dna.R
+43b881eea3763d0079e88012e0b84d37 *R/write.dna.R
4e5d7a7904fe6496a44c4b69f674290b *R/write.nexus.R
89496c0df70205654456d3b7bb7ba41f *R/write.nexus.data.R
17d72a136a8131ea46ef9a20fcfd4d36 *R/write.tree.R
-a918c086a449bcca5ccbb04f7e6d50a9 *R/yule.R
+774ce72875903259aade5344f9a70aa4 *R/yule.R
c8d3aa3fe64e75e61af07a1b11c74f3f *R/yule.time.R
1eb44ff9e5a036eb845faa1598ce5009 *R/zoom.R
3387c0d0c1f913f8471e1bb34bd2e516 *R/zzz.R
-e5443aa31e5f8c448d6d3458b580b2e6 *build/vignette.rds
+b77a7eef47f2f8d2874a638562f49ea4 *build/vignette.rds
db9083e8750aff839d5ebf3ed982f1f1 *data/HP.links.rda
9d9f9232839665422709ded1e541d038 *data/bird.families.rda
a14a6df0f3a735ebc056065077788c90 *data/bird.orders.rda
@@ -131,22 +133,20 @@ f74f9ed80c04756021cc093d40ca9ff9 *data/carnivora.csv.gz
4eaf8cbaefa2e8f8d395a9b482ee9967 *data/chiroptera.rda
1c74c3b99d08b0e17eea3ec1065c12d2 *data/cynipids.rda
7fe760c2f3b4deba0554aae6138cb602 *data/gopher.D.rda
-3a0678a086d0d8ff8c1a722fd7baac19 *data/hivtree.newick.rda
+1c8e6ca4c99ce95b2f2692c15b6d768d *data/hivtree.newick.rda
8d14f95319d0a5cdc8faa60a1d0085ce *data/hivtree.table.txt.gz
-41e2b0f8336384bf030aca421c987d0f *data/landplants.newick.rda
31be81fe3faca11f98d3e74c090bc59e *data/lice.D.rda
38edbd84a0a067322c40db8d71fb1289 *data/lmorigin.ex1.rda
e3ce9e3444182fea2e65df2e150ea0db *data/lmorigin.ex2.rda
ce7a56faebdf286fdf5ba6c8c3699a79 *data/mat3.RData
e2d1339025ed901009bfed58dc6505ff *data/mat5M3ID.RData
101d0ab2e981b0987cde704a2dee1d8d *data/mat5Mrand.RData
-68e738ac2bee2dd81f65bc1d138e25f4 *data/opsin.newick.rda
39e4fece2bdc527d7a9d4d83d023a947 *data/woodmouse.rda
828290996b613493f96e0fab024beebb *inst/CITATION
3f54f3775bcf382e25df2a12228894f6 *inst/doc/MoranI.R
2277b0efdb2f70dfdc5df8278702c6b6 *inst/doc/MoranI.Rnw
-a5c7aca72eeb81a034cebbc363457a5f *inst/doc/MoranI.pdf
-e6b3b19323b56c8e10cffc858b3af117 *man/AAbin.Rd
+31ca3eafda127eab7c4799764f8fa68c *inst/doc/MoranI.pdf
+db2ed5302bddcfac232dcd6fabe4a651 *man/AAbin.Rd
e6876b193a0df06697c788a8e48cf4bc *man/CADM.global.Rd
dfa15a3e3bb57c9b21a30b8d5d790c62 *man/DNAbin.Rd
d94f358593695b1713840df5a8c000ba *man/DNAbin2indel.Rd
@@ -159,10 +159,11 @@ b4b2d3c708bf1ebca6f1f3b4cb752243 *man/ace.Rd
60adac0f324c7b08a756d8f45c02b024 *man/add.scale.bar.Rd
0d68dd79c42bad2f3a68635ba75a59b0 *man/additive.Rd
25a2859708b9a281a0d682e4376f3f53 *man/alex.Rd
+2ff5d30c6fb1c5458643f1b3c09e76da *man/all.equal.DNAbin.Rd
f0b85d459edb5ead2e6d39a86cbb350c *man/all.equal.phylo.Rd
-c987c7b1cf09a272a6ce77eb3e9554dc *man/alview.Rd
-51a4ced7c3c2533b34887117e566c46f *man/ape-internal.Rd
-587cbd9d9a326b3990ed20bfe5016165 *man/ape-package.Rd
+c3bd9e7b02831d6086fabd5de6db917c *man/alview.Rd
+42489c977f59031dc17a0ba7e55663b6 *man/ape-internal.Rd
+f84e59c805efa00000d775b81bf9f8d0 *man/ape-package.Rd
5bba4ae4bfc66b613855cfc182d9b1bc *man/as.alignment.Rd
2771f272ad46b279022892d9f5d21eb2 *man/as.bitsplits.Rd
4f014cf2923e2eab6188acd48e8096fa *man/as.matching.Rd
@@ -179,18 +180,18 @@ f929bc1b6391c57a6b0099c4561fd7be *man/binaryPGLMM.Rd
71a008cfe65c4f524a5b66e68bbf81ab *man/bird.families.Rd
0e41770e1e6d0b8d90c4cf51049213cb *man/bird.orders.Rd
ef1c15d5d93410c21179997431112209 *man/birthdeath.Rd
-22fbad77d4728f10b2e117e2c6a53efe *man/boot.phylo.Rd
+87fafa0226908ea87c6d11aecc9bf72c *man/boot.phylo.Rd
5a64b90d3a6c7a8204946b00f45f4cfc *man/branching.times.Rd
-7dad9d632d914de8122f2455f0762546 *man/c.phylo.Rd
+c808768a942998429ecc9fc925f3f3a3 *man/c.phylo.Rd
9d3f9614732671a610cc8a37454885e2 *man/carnivora.Rd
-7c4994dab4cfa672ad54e945011fa560 *man/checkAlignment.Rd
+953b3c21c60a38e8e53ec6894170dd36 *man/checkAlignment.Rd
5ff8c7e8fad519d978f166948c03059c *man/checkValidPhylo.Rd
64c3996ca6bcc97d0d2e2cf3361f8f71 *man/cherry.Rd
f82a24dc4bb17b5675ec7b422d575587 *man/chiroptera.Rd
b40dd77fd2529626e983ee0300be6781 *man/chronoMPL.Rd
c1f01c6200b2f1e2901d45d40daae404 *man/chronopl.Rd
506d0332fb092ab87ca7674faef63ab7 *man/chronos.Rd
-61a3ac5c65d72509d2ef3fcb56507f30 *man/clustal.Rd
+1edeca1697a517d946b2062d20b90cbb *man/clustal.Rd
866af6e8d769b3d6972ef8e1ac849a12 *man/coalescent.intervals.Rd
18c77b4104aa22d3e6cb88339f2eff90 *man/collapse.singles.Rd
bff5a7826f5a39767601e32ceb776247 *man/collapsed.intervals.Rd
@@ -218,7 +219,7 @@ c0763a70c4965a6b03df3e5be68e450d *man/def.Rd
fbcd1d4bcf74e21fc93c195c7af3db98 *man/delta.plot.Rd
d4c1853e139f84101de2b9309823f366 *man/dist.dna.Rd
38011e81d28a120d88eead09e62c154a *man/dist.gene.Rd
-e813d96ce8d458a535fffbb072701300 *man/dist.topo.Rd
+4042a6c6113e8a0eb3a1ffaf51138f0f *man/dist.topo.Rd
c7cc398115be066740ca4fb037394727 *man/diversi.gof.Rd
d646ea0343999bd0e38e86dcf6c12018 *man/diversi.time.Rd
da8898476bb15b627b34ee1093b9aeb4 *man/diversity.contrast.test.Rd
@@ -231,15 +232,14 @@ eea313e8ee32597b4cec120d23113642 *man/gammaStat.Rd
3991fa7864e326579f1ab8b671095e4b *man/hivtree.Rd
18012643a904f657fc5f5896f3d14054 *man/howmanytrees.Rd
86c49d080fdffd614d8056021e91cc55 *man/identify.phylo.Rd
-5edaf778ad69fe8c9984bf3817a86a88 *man/image.DNAbin.Rd
-9957425ef61231ba72538e771c6565b5 *man/is.binary.tree.Rd
+dc00b2b43a44f4ba7ea76ea23a12c615 *man/image.DNAbin.Rd
+0ede757d7bed82216980d3c4fd592cb6 *man/is.binary.tree.Rd
80a5a228a43679edf75903590253a875 *man/is.compatible.Rd
d2de8fd9549ef01a1dddeb726dd77fcf *man/is.monophyletic.Rd
-9d09e72a00da6b7b3f967f76166077cd *man/is.ultrametric.Rd
+ad9e7316219c3238b44b02d55c44b4d3 *man/is.ultrametric.Rd
3f6ff8340b6c9770a1f4d2fed72a295d *man/kronoviz.Rd
8edde42ffc43737ba82e18594801421f *man/label2table.Rd
2afa6305e48e4c47a7d94d46401f85a3 *man/ladderize.Rd
-23db31affc0dc60f6c772e850bb640c3 *man/landplants.Rd
4b86460e4e5d993abf13c99b9744dbd6 *man/lmorigin.Rd
21bb31db5bcc8d8650467ef73fe0c0d3 *man/ltt.plot.Rd
7ca1aeffd0cd8106e00d0ee1952b5ccb *man/makeLabel.Rd
@@ -251,18 +251,18 @@ f56f6f49c89c6bc850a7aee1bce5e0bd *man/mat5M3ID.Rd
69ae0cb181240bb8ec168e69f1ba44bb *man/matexpo.Rd
a8b9d6b04d35d75f43d1b283361a1642 *man/mcconwaysims.test.Rd
8c3dfdb6c70caee1051b1e576e3eaa74 *man/mcmc.popsize.Rd
-c1d52335ab9d9ccac265c06c9fae019c *man/mixedFontLabel.Rd
+d974e980cb56568e337ef4e926f43617 *man/mixedFontLabel.Rd
dbcdc231c2f1f221818307bb33486159 *man/mrca.Rd
5c88230ad597ea9fe41536a8321c326b *man/mst.Rd
-114f7834d51151cb6f92f8ed77becbcd *man/multi2di.Rd
+b74a379a0a9b80941118922ca7843b02 *man/multi2di.Rd
029b04adeafe89ea5edf9a1ab00cd154 *man/multiphylo.Rd
00fb7ade93c2dd887be27e3dad8e2115 *man/mvr.Rd
3df9e16b8a09df3f1dba5c4327a635fc *man/nj.Rd
9ea7d5899a190171de4165e3240de62e *man/njs.Rd
+8a9189a62be3d0e871deb20be5bd6e8a *man/node.dating.Rd
a589b4cc04505185dc9ef1817c7ae304 *man/node.depth.Rd
ba754eba67181a1c39d417ffd9d23b76 *man/nodelabels.Rd
447ae03684ff56a4a24932aec182acf3 *man/nodepath.Rd
-4f8b4cf2a725a8001e49390056181bef *man/opsin.Rd
c2e2f35f4e233265c86b7967ec2d0630 *man/parafit.Rd
89db814a5fc2e9fa5f0b9f4ad2f0d070 *man/pcoa.Rd
583200254fe87b2a36fc18a88825298d *man/phydataplot.Rd
@@ -271,23 +271,25 @@ c2e2f35f4e233265c86b7967ec2d0630 *man/parafit.Rd
0363aa3fa3e86160934e2359c8ac5323 *man/pic.ortho.Rd
265527313d479d3625df7680c7479cd1 *man/plot.correlogram.Rd
d23d21863b2dffa5388f3fc128368af4 *man/plot.phylo.Rd
-a66c2f038b7657e5d619dee5ee9ed9ba *man/plot.phyloExtra.Rd
+896198fa3a0ce916346c8b507bf707bf *man/plot.phyloExtra.Rd
6bb1cf751772b9b6acf560b5f3c6f3c1 *man/plot.varcomp.Rd
+e0e2f60baddaaaa8f1d9796bd50a8ba2 *man/plotTreeTime.Rd
b24438c42cea969302ec6ba61002426e *man/print.phylo.Rd
80bf68fd21b5e665329cf30b10a1655f *man/rTraitCont.Rd
59e81eaae91dc77b4275c9a5b2910dda *man/rTraitDisc.Rd
81f756fdf2ec4c968095595dece8338f *man/rTraitMult.Rd
-38877ac887281607cb0a1bcd6325f983 *man/read.GenBank.Rd
+028b11582b3493cdefea60ccb78ad429 *man/read.GenBank.Rd
9e31bccfa08c1d2da2180cd8f50ab250 *man/read.caic.Rd
-0c91321b05d51c501e0f5c51de78538c *man/read.dna.Rd
+2b32811be5bed71ffeca3801c4408e5b *man/read.dna.Rd
+51670875393d70cb7ef84b4b0388e8e4 *man/read.gff.Rd
5c9aed0b6f0a7e0f0fac9b5a59ac732f *man/read.nexus.Rd
bc02e36c51d67074e661468993ed359b *man/read.nexus.data.Rd
3c8b53f5af4fdd36e7c0b449412d5a77 *man/read.tree.Rd
-48d08588177df45a1e53dede1dc56bae *man/reconstruct.Rd
-61e0bc86ae5a1f04ee6e88cdc92a2089 *man/reorder.phylo.Rd
+9224f164351fea4be9f7ce9d827f8e58 *man/reconstruct.Rd
+5c2f0e05ef560aa4ec066c898c430aec *man/reorder.phylo.Rd
af19f262b468dfd7c12448d9d8085331 *man/richness.yule.test.Rd
3afe25cebe4628e4844e41aad512d6e4 *man/rlineage.Rd
-0f62fb19db3810dfc66da173fec46104 *man/root.Rd
+c2503e4a9511cafefa7bc1e058fdde2a *man/root.Rd
6b97ea2fd96bf948fca764deab3a6e76 *man/rotate.Rd
1a7314ee0aa76712737e76162df6cc74 *man/rtree.Rd
e53b4b3f7e009289c628e16143b3e8b4 *man/rtt.Rd
@@ -299,7 +301,7 @@ a0db170675b1c265f1f4e7d56e88d7cb *man/speciesTree.Rd
8fc51ff26614e7c713be2570d9af47b6 *man/stree.Rd
1f1309e2ec6327974952d546241a6024 *man/subtreeplot.Rd
ef8aa9721275f20ddd38612cddf8300c *man/subtrees.Rd
-0091fcdcb4aa48aa364da2512545f6c1 *man/summary.phylo.Rd
+d684cc640acf041c6f1124e08426c8c8 *man/summary.phylo.Rd
18eb8ab87509076b49c50cd86b231cf7 *man/trans.Rd
3340a4f09c55010c15c0149a20b41c99 *man/treePop.Rd
7749863f505076f2288c106eb465805e *man/trex.Rd
@@ -327,14 +329,14 @@ a00006ae345bb9379312e81694de3885 *man/zoom.Rd
122ba51b574f7f7acd5ea4b6256cea5f *src/SPR.c
d11fb19958f2953583a9c5876da82df9 *src/TBR.c
9733c82cd67c4bd41aea44981f4ac8b8 *src/additive.c
-91ee394e252f7eec643d899f60ac4b76 *src/ape.c
+f795f77892d496e05e247b114eacf83a *src/ape.c
7eaffd3d85df9a1e4460905e9ca5eced *src/ape.h
19b564aab62c464c092cdf2f0d5cf447 *src/bNNI.c
79ca5c473decf13964192524caebf9f1 *src/bionjs.c
f782a97466ed8008d31afd7a3e9d4ea9 *src/bipartition.c
a5d4f692aca36f6b50e123b443792fa1 *src/bitsplits.c
81e4ad60b418966cdd09e608e074bb28 *src/delta_plot.c
-58108a4ba25655f5b5101ebdfea78d1e *src/dist_dna.c
+bcbc45836832bf51f35e920283cc37d6 *src/dist_dna.c
2a6a59a3a220eb855a3c48bc73dc54ba *src/dist_nodes.c
005ab69356525fbbc1b69950341308c2 *src/ewLasso.c
afa27403972145d083d7f7f8b2536b98 *src/heap.c
@@ -350,7 +352,7 @@ d3364013ad597425d01ac50ae4ee0b34 *src/nj.c
7de13de6d828d870915c2efb6eb47ab6 *src/pic.c
28ef72b2ccede817c0431afcb04e9c8c *src/plot_phylo.c
aa8d9966da3b7e970879a49a40c05a07 *src/rTrait.c
-ab68729ab24b62eda61e1a2bd4a4481e *src/read_dna.c
+23d15137e1f3d9b0d365d04c1f4d9b05 *src/read_dna.c
b5505d4f7c732c82bc54687d6fc37c0d *src/reorder_phylo.c
d7d48424f600f5dad372a8c3ccfbbcad *src/treePop.c
40d4ed144aba4b1acbc64d2807f33466 *src/tree_build.c
diff --git a/NAMESPACE b/NAMESPACE
index 1b635ff..5f76a37 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -14,20 +14,23 @@ importFrom(lattice, xyplot, panel.lines, panel.points)
importFrom(methods, as, show)
importFrom(nlme, corMatrix, Dim, getCovariate, getGroups,
getGroupsFormula, gls, Initialize)
-importFrom(stats, as.dist, as.hclust, biplot, coef, complete.cases,
+importFrom(stats, AIC, as.dist, as.hclust, biplot, coef, complete.cases,
cophenetic, cor, cor.test, cov, cov2cor, density, dgamma,
dpois, drop1, formula, gaussian, glm, hclust, integrate,
lm, mahalanobis, median, model.frame, model.matrix,
- model.response, na.fail, nlm, nlminb, optim, optimize, p.adjust,
- pchisq, pf, pgamma, pnorm, ppois, printCoefmat, pt, qbinom, qnorm,
- qt, quantile, quasibinomial, rbinom, reorder, resid, rexp, rgamma,
- rnorm, runif, sd, terms, uniroot, var, wilcox.test)
-importFrom(utils, download.file, read.table, setTxtProgressBar, str,
- txtProgressBar)
+ model.response, na.fail, na.omit, nlm, nlminb, optim, optimize,
+ p.adjust, pchisq, pf, pgamma, pnorm, ppois, printCoefmat, pt,
+ qbinom, qnorm, qt, quantile, quasibinomial, rbinom, reorder, resid,
+ rexp, rgamma, rnorm, runif, sd, setNames, terms, uniroot, var,
+ wilcox.test)
+importFrom(utils, download.file, read.table, str)
+importFrom(parallel, mclapply)
## Methods for the classes defined in ape, including for the generics
## defined in ape (see also below). Some methods are exported above.
+S3method(is.binary, tree) # to delete when removing the function
+
S3method("[", AAbin)
S3method(as.character, AAbin)
S3method(image, AAbin)
@@ -54,6 +57,7 @@ S3method(print, compar.gee)
S3method(print, corphylo)
S3method("[", DNAbin)
+S3method(all.equal, DNAbin)
S3method(as.character, DNAbin)
S3method(as.list, DNAbin)
S3method(as.matrix, DNAbin)
@@ -76,11 +80,22 @@ S3method("[[<-", multiPhylo)
S3method("$", multiPhylo)
S3method("$<-", multiPhylo)
S3method(c, multiPhylo)
+S3method(di2multi, multiPhylo)
+S3method(is.binary, multiPhylo)
+S3method(is.rooted, multiPhylo)
+S3method(is.ultrametric, multiPhylo)
S3method(makeLabel, multiPhylo)
+S3method(multi2di, multiPhylo)
+S3method(Nedge, multiPhylo)
+S3method(Nnode, multiPhylo)
+S3method(Ntip, multiPhylo)
S3method(plot, multiPhylo)
S3method(print, multiPhylo)
+S3method(reorder, multiPhylo)
+S3method(root, multiPhylo)
S3method(str, multiPhylo)
S3method(unique, multiPhylo)
+S3method(unroot, multiPhylo)
S3method("+", phylo)
S3method(all.equal, phylo)
@@ -89,13 +104,23 @@ S3method(as.matching, phylo)
S3method(coalescent.intervals, phylo)
S3method(c, phylo)
S3method(cophenetic, phylo)
+S3method(di2multi, phylo)
S3method(identify, phylo)
+S3method(is.binary, phylo)
+S3method(is.rooted, phylo)
+S3method(is.ultrametric, phylo)
S3method(makeLabel, phylo)
+S3method(multi2di, phylo)
+S3method(Nedge, phylo)
+S3method(Nnode, phylo)
+S3method(Ntip, phylo)
S3method(plot, phylo)
S3method(print, phylo)
S3method(reorder, phylo)
+S3method(root, phylo)
S3method(skyline, phylo)
S3method(summary, phylo)
+S3method(unroot, phylo)
S3method(vcv, phylo)
S3method(plot, phymltest)
diff --git a/R/CDF.birth.death.R b/R/CDF.birth.death.R
index 2cced20..dd0cabf 100644
--- a/R/CDF.birth.death.R
+++ b/R/CDF.birth.death.R
@@ -1,8 +1,8 @@
-## CDF.birth.death.R (2015-07-28)
+## CDF.birth.death.R (2016-07-06)
## Functions to Simulate and Fit Time-Dependent Birth-Death Models
-## Copyright 2010-2015 Emmanuel Paradis
+## Copyright 2010-2016 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
@@ -579,7 +579,7 @@ rphylo <-
inc <- 10
x <- t - inc
while (inc > eps) {
- if (Foo(x, t) > P) {
+ if (Foo(t, x) > P) { # fixed by Niko Yasui (2016-07-06)
x <- x + inc
inc <- inc/10
}
diff --git a/R/DNA.R b/R/DNA.R
index 15367d5..0757da6 100644
--- a/R/DNA.R
+++ b/R/DNA.R
@@ -1,4 +1,4 @@
-## DNA.R (2016-02-18)
+## DNA.R (2016-11-26)
## Manipulations and Comparisons of DNA and AA Sequences
@@ -13,6 +13,8 @@ DNAbin2indel <- function(x)
d <- dim(x)
s <- as.integer(d[2])
n <- as.integer(d[1])
+ if (s * n > 2^31 - 1)
+ stop("DNAbin2indel() cannot handle more than 2^31 - 1 bases")
res <- .C(DNAbin2indelblock, x, n, s, integer(n*s), NAOK = TRUE)[[4]]
dim(res) <- d
rownames(res) <- rownames(x)
@@ -106,14 +108,13 @@ as.matrix.DNAbin <- function(x, ...)
dim(x) <- c(1, length(x))
return(x)
}
- s <- unique(unlist(lapply(x, length)))
+ s <- unique(lengths(x, use.names = FALSE))
if (length(s) != 1)
stop("DNA sequences in list not of the same length.")
- nms <- names(x)
n <- length(x)
- y <- matrix(as.raw(0), n, s)
+ y <- matrix(raw(), n, s)
for (i in seq_len(n)) y[i, ] <- x[[i]]
- rownames(y) <- nms
+ rownames(y) <- names(x)
class(y) <- "DNAbin"
y
}
@@ -125,7 +126,7 @@ as.list.DNAbin <- function(x, ...)
else { # matrix
n <- nrow(x)
obj <- vector("list", n)
- for (i in 1:n) obj[[i]] <- x[i, ]
+ for (i in seq_len(n)) obj[[i]] <- x[i, , drop = TRUE]
names(obj) <- rownames(x)
}
class(obj) <- "DNAbin"
@@ -216,8 +217,8 @@ print.DNAbin <- function(x, printlen = 6, digits = 3, ...)
cat("Label:", nms, "\n\n")
} else {
cat(n, "DNA sequences in binary format stored in a list.\n\n")
- tmp <- unlist(lapply(x, length))
- nTot <- sum(tmp)
+ tmp <- lengths(x, use.names = FALSE)
+ nTot <- sum(as.numeric(tmp))
mini <- range(tmp)
maxi <- mini[2]
mini <- mini[1]
@@ -312,19 +313,17 @@ as.character.DNAbin <- function(x, ...)
base.freq <- function(x, freq = FALSE, all = FALSE)
{
- if (!inherits(x, "DNAbin")) stop('base.freq requires an object of class "DNAbin"')
+ if (!inherits(x, "DNAbin"))
+ stop('base.freq requires an object of class "DNAbin"')
- f <- function(x)
- .C(BaseProportion, x, as.integer(length(x)), double(17),
- NAOK = TRUE)[[3]]
+ f <- function(x) .Call(BaseProportion, x)
if (is.list(x)) {
BF <- rowSums(sapply(x, f))
- n <- sum(sapply(x, length))
+ n <- sum(as.double(lengths(x, use.names = FALSE)))
} else {
n <- length(x)
- BF <- .C(BaseProportion, x, as.integer(n), double(17),
- NAOK = TRUE)[[3]]
+ BF <- f(x)
}
names(BF) <- c("a", "c", "g", "t", "r", "m", "w", "s",
@@ -374,16 +373,10 @@ GC.content <- function(x) sum(base.freq(x)[2:3])
seg.sites <- function(x)
{
if (is.list(x)) x <- as.matrix(x)
- if (is.vector(x)) n <- 1
- else { # 'x' is a matrix
- n <- dim(x)
- s <- n[2]
- n <- n[1]
- }
- if (n == 1) return(integer(0))
- ans <- .C(SegSites, x, as.integer(n), as.integer(s),
- integer(s), NAOK = TRUE)
- which(as.logical(ans[[4]]))
+ if (is.vector(x)) return(integer(0))
+ if (nrow(x) == 1) return(integer(0))
+ ans <- .Call(SegSites, x)
+ which(as.logical(ans))
}
dist.dna <- function(x, model = "K80", variance = FALSE, gamma = FALSE,
@@ -411,6 +404,9 @@ dist.dna <- function(x, model = "K80", variance = FALSE, gamma = FALSE,
s <- n[2]
n <- n[1]
+ if (s * n > 2^31 - 1)
+ stop("dist.dna() cannot handle more than 2^31 - 1 bases")
+
if (imod %in% c(4, 6:8)) {
BF <- if (is.null(base.freq)) base.freq(x) else base.freq
} else BF <- 0
@@ -458,6 +454,7 @@ image.DNAbin <- function(x, what, col, bg = "white", xlab = "", ylab = "",
if (missing(what)) c("a", "g", "c", "t", "n", "-") else tolower(what)
if (missing(col))
col <- c("red", "yellow", "green", "blue", "grey", "black")
+ x <- as.matrix(x) # tests if all sequences have the same length
n <- (dx <- dim(x))[1] # number of sequences
s <- dx[2] # number of sites
y <- integer(N <- length(x))
@@ -489,10 +486,10 @@ image.DNAbin <- function(x, what, col, bg = "white", xlab = "", ylab = "",
brks <- c(-0.5, brks)
}
yaxt <- if (show.labels) "n" else "s"
- image.default(1:s, 1:n, t(y), col = co, xlab = xlab,
+ image.default(1:s, 1:n, t(y[n:1, ]), col = co, xlab = xlab,
ylab = ylab, yaxt = yaxt, breaks = brks, ...)
if (show.labels)
- mtext(rownames(x), side = 2, line = 0.1, at = 1:n,
+ mtext(rownames(x), side = 2, line = 0.1, at = n:1,
cex = cex.lab, adj = 1, las = 1)
if (legend) {
psr <- par("usr")
@@ -527,26 +524,24 @@ alview <- function(x, file = "", uppercase = TRUE)
where <- function(x, pattern)
{
pat <- as.DNAbin(strsplit(pattern, NULL)[[1]])
- p <- as.integer(length(pat))
+ p <- length(pat)
- foo <- function(x, pat, p) {
- s <- as.integer(length(x))
+ f <- function(x, pat, p) {
+ s <- length(x)
if (s < p) stop("sequence shorter than the pattern")
- ans <- .C(C_where, x, pat, s, p, integer(s), 0L, NAOK = TRUE)
- n <- ans[[6]]
- if (n) ans[[5]][seq_len(n)] - p + 2L else integer()
+ .Call(C_where, x, pat)
}
- if (is.list(x)) return(lapply(x, foo, pat = pat, p = p))
+ if (is.list(x)) return(lapply(x, f, pat = pat, p = p))
if (is.matrix(x)) {
n <- nrow(x)
res <- vector("list", n)
for (i in seq_len(n))
- res[[i]] <- foo(x[i, , drop = TRUE], pat, p)
+ res[[i]] <- f(x[i, , drop = TRUE], pat, p)
names(res) <- rownames(x)
return(res)
}
- foo(x, pat, p) # if x is a vector
+ f(x, pat, p) # if x is a vector
}
## conversions from BioConductor:
@@ -632,7 +627,7 @@ print.AAbin <- function(x, ...)
cat(n, "amino acid sequence")
if (n > 1) cat("s")
cat(" in a list\n\n")
- tmp <- unlist(lapply(x, length))
+ tmp <- lengths(x, use.names = FALSE)
maxi <- max(tmp)
mini <- min(tmp)
if (mini == maxi)
@@ -668,7 +663,7 @@ as.character.AAbin <- function(x, ...)
{
f <- function(x) strsplit(rawToChar(x), "")[[1]]
if (is.list(x)) return(lapply(x, f))
- if (is.matrix(x)) return(apply(t(x), 1, f))
+ if (is.matrix(x)) return(t(apply(x, 1, f)))
f(x)
}
@@ -712,7 +707,7 @@ dist.aa <- function(x, pairwise.deletion = FALSE, scaled = FALSE)
p <- length(b <- b[!del])
tmp <- sum(a[!del] != b)
k <- k + 1L
- d[k] <- if (scaled) tmp else tmp/p
+ d[k] <- if (scaled) tmp/p else tmp
}
}
}
@@ -774,18 +769,17 @@ image.AAbin <- function(x, what, col, bg = "white", xlab = "", ylab = "",
if (sm == N) {
leg.co <- co <- col
leg.txt <- what
- }
- else {
+ } else {
co <- c(bg, col)
leg.txt <- c(what, "Unknown")
leg.co <- c(col, bg)
brks <- c(-0.5, brks)
}
yaxt <- if (show.labels) "n" else "s"
- image.default(1:s, 1:n, t(y), col = co, xlab = xlab, ylab = ylab,
+ image.default(1:s, 1:n, t(y[n:1, ]), col = co, xlab = xlab, ylab = ylab,
yaxt = yaxt, breaks = brks, ...)
if (show.labels)
- mtext(rownames(x), side = 2, line = 0.1, at = 1:n, cex = cex.lab,
+ mtext(rownames(x), side = 2, line = 0.1, at = n:1, cex = cex.lab,
adj = 1, las = 1)
if (legend) {
psr <- par("usr")
@@ -809,7 +803,7 @@ if (check.gaps) {
else {
rest <- gap.length %% 3
if (any(cond <- rest > 0)) {
- cat("Some gap length are not multiple of 3:", gap.length[cond])
+ cat("Some gap lengths are not multiple of 3:", gap.length[cond])
} else cat("All gap lengths are multiple of 3.")
tab <- tabulate(y, gap.length[length(gap.length)])
@@ -818,6 +812,20 @@ if (check.gaps) {
names(tab) <- gap.length
print(tab)
+ ## find gaps on the borders:
+ col1 <- unique(y[, 1])
+ if (!col1[1]) col1 <- col1[-1]
+ if (length(col1))
+ cat(" => length of gaps on the left border of the alignment:", unique(col1), "\n")
+ else cat(" => no gap on the left border of the alignment\n")
+ i <- which(y != 0, useNames = FALSE)
+ jcol <- i %/% nrow(y) + 1
+ yi <- y[i]
+ j <- yi == s - jcol + 1
+ if (any(j))
+ cat(" => length of gaps on the right border of the alignment:", yi[j], "\n")
+ else cat(" => no gap on the right border of the alignment\n")
+
## find base segments:
A <- B <- numeric()
for (i in seq_len(n)) {
@@ -825,17 +833,17 @@ if (check.gaps) {
if (!length(j)) next
k <- j + y[i, j] # k: start of each base segment in the i-th sequence
if (j[1] != 1) k <- c(1, k) else j <- j[-1]
- if (k[length(k)] > s) k <- k[-length(k)] else j <- c(j, s)
+ if (k[length(k)] > s) k <- k[-length(k)] else j <- c(j, s + 1)
A <- c(A, j)
B <- c(B, k)
}
AB <- unique(cbind(A, B))
not.multiple.of.3 <- (AB[, 1] - AB[, 2]) %% 3 != 0
left.border <- AB[, 2] == 1
- right.border <- AB[, 1] == s
+ right.border <- AB[, 1] == s + 1
Nnot.mult3 <- sum(not.multiple.of.3)
- cat("\nNumber of unique contiguous base segments:", nrow(AB), "\n")
- if (!Nnot.mult3) cat("All segment lengths multiple of 3")
+ cat("\nNumber of unique contiguous base segments defined by gaps:", nrow(AB), "\n")
+ if (!Nnot.mult3) cat("All segment lengths multiple of 3.\n")
else {
Nleft <- sum(not.multiple.of.3 & left.border)
Nright <- sum(not.multiple.of.3 & right.border)
@@ -845,7 +853,7 @@ if (check.gaps) {
if (Nright + Nleft < Nnot.mult3) {
cat(" => positions of these segments inside the alignment: ")
sel <- not.multiple.of.3 & !left.border & !right.border
- cat(paste(AB[sel, 2], AB[sel, 1] - 1, sep = "--"), "\n")
+ cat(paste(AB[sel, 2], AB[sel, 1] - 1, sep = ".."), "\n")
}
}
}
@@ -886,3 +894,99 @@ if (check.gaps) {
plot(1:s, G, "h", xlab = "Sequence position", ylab = "Number of observed bases")
}
}
+
+all.equal.DNAbin <- function(target, current, plot = FALSE, ...)
+{
+ if (identical(target, current)) return(TRUE)
+ name.target <- deparse(substitute(target))
+ name.current <- deparse(substitute(current))
+ st1 <- "convert list as matrix for further comparison."
+# st2 <- ""
+ st3 <- "Subset your data for further comparison."
+ isali1 <- is.matrix(target)
+ isali2 <- is.matrix(current)
+ if (isali1 && !isali2)
+ return(c("1st object is a matrix, 2nd object is a list:", st1))
+ if (!isali1 && isali2)
+ return(c("1st object is a list, 2nd object is a matrix:", st1))
+ if (!isali1 && !isali2)
+ return(c("Both objects are lists:",
+ "convert them as matrices for further comparison."))
+# n1 <- if (isali1) nrow(target) else length(target)
+# n2 <- if (isali2) nrow(current) else length(current)
+ if (ncol(target) != ncol(current))
+ return("Numbers of columns different: comparison stopped here.")
+
+ foo <- function(n) ifelse(n == 1, "sequence", "sequences")
+ doComparison <- function(target, current)
+ which(target != current, arr.ind = TRUE, useNames = FALSE)
+
+ n1 <- nrow(target)
+ n2 <- nrow(current)
+ labs1 <- labels(target)
+ labs2 <- labels(current)
+ if (identical(labs1, labs2)) {
+ res <- "Labels in both objects identical."
+ res <- list(messages = res,
+ different.sites = doComparison(target, current))
+ } else {
+ in12 <- labs1 %in% labs2
+ in21 <- labs2 %in% labs1
+ if (n1 != n2) {
+ res <- c("Number of sequences different:",
+ paste(n1, foo(n1), "in 1st object;", n2,
+ foo(n2), "in 2nd object."),
+ st3)
+ plot <- FALSE
+ } else { # n1 == n2
+ if (any(!in12)) {
+ res <- c("X: 1st object (target), Y: 2nd object (current).",
+ paste("labels in X not in Y:",
+ paste(labs1[!in12], collapse = ", ")),
+ paste("labels in X not in Y:",
+ paste(labs2[!in21], collapse = ", ")),
+ st3)
+ plot <- FALSE
+ } else {
+ res <- c("Labels in both objects identical but not in the same order.",
+ "Comparing sequences after reordering rows of the second matrix.")
+ current <- current[labs1, ]
+ if (identical(target, current)) {
+ res <- c(res, "Sequences are identical.")
+ plot <- FALSE
+ } else {
+ res <- list(messages = res,
+ different.sites = doComparison(target, current))
+ }
+ }
+ }
+ }
+ if (plot) {
+ cols <- unique(res$different.sites[, 2])
+ diff.cols <- diff(cols)
+ j <- which(diff.cols != 1)
+ end <- c(cols[j], cols[length(cols)])
+ start <- c(cols[1], cols[j + 1])
+ v <- cumsum(end - start + 1) + 0.5
+ f <- function(lab) {
+ axis(2, at = seq_len(n1), labels = FALSE)
+ axis(1, at = seq_along(cols), labels = cols)
+ mtext(lab, line = 1, adj = 0, font = 2)
+ }
+ layout(matrix(1:2, 2))
+ par(xpd = TRUE)
+ image(target[, cols], show.labels = FALSE, axes = FALSE, ...)
+ f(name.target)
+ xx <- c(0.5, v)
+ segments(xx, 0.5, xx, n1, lty = 2, col = "white", lwd = 2)
+ segments(xx, 0.5, xx, -1e5, lty = 2, lwd = 2)
+ image(current[, cols], show.labels = FALSE, axes = FALSE, ...)
+ f(name.current)
+ segments(xx, 0.5, xx, n2, lty = 2, col = "white", lwd = 2)
+ segments(xx, 1e5, xx, n2, lty = 2, lwd = 2)
+ #segments(0.5, -5, length(cols) + 0.5, -5, lwd = 5, col = "grey")
+ #rect(0.5, -4, length(cols) + 0.5, -3, col = "grey")
+ #segments(0.5, 0.5, 10, -3)
+ }
+ res
+}
diff --git a/R/MPR.R b/R/MPR.R
index 40cc41b..2d733cc 100644
--- a/R/MPR.R
+++ b/R/MPR.R
@@ -11,7 +11,7 @@ MPR <- function(x, phy, outgroup)
{
if (is.rooted(phy))
stop("the tree must be unrooted")
- if (!is.binary.tree(phy))
+ if (!is.binary.phylo(phy))
stop("the tree must be fully dichotomous")
if (length(outgroup) > 1L)
stop("outgroup must be a single tip")
diff --git a/R/SDM.R b/R/SDM.R
index d2a6686..e58f28d 100644
--- a/R/SDM.R
+++ b/R/SDM.R
@@ -20,7 +20,7 @@ SDM <- function(...)
ROWNAMES <- lapply(st[ONEtoK], rownames)
## the number of rows of each matrix:
- NROWS <- sapply(ROWNAMES, length)
+ NROWS <- lengths(ROWNAMES)
tot <- sum(NROWS)
labels <- unique(unlist(ROWNAMES))
diff --git a/R/SlowinskiGuyer.R b/R/SlowinskiGuyer.R
index 353c467..d9aad9b 100644
--- a/R/SlowinskiGuyer.R
+++ b/R/SlowinskiGuyer.R
@@ -1,8 +1,8 @@
-## SlowinskiGuyer.R (2011-05-05)
+## SlowinskiGuyer.R (2016-10-23)
## Tests of Diversification Shifts with Sister-Clades
-## Copyright 2011 Emmanuel Paradis
+## Copyright 2011-2016 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
@@ -43,22 +43,23 @@ richness.yule.test <- function(x, t)
tb <- c(t, t)
.PrNt.Yule <- function(N, age, birth) {
- tmp <- exp(-birth * age)
- tmp * (1 - tmp)^(N - 1)
+ tmp <- -birth * age
+ tmp + (N - 1) * log(1 - exp(tmp)) # on a log-scale
}
- minusloglik0 <- function(l)
- -sum(log(.PrNt.Yule(n, tb, l)))
+ ## the functions to minimize:
+ minusloglik0 <- function(l) -sum(.PrNt.Yule(n, tb, l))
+ minusloglika <- function(l) -sum(.PrNt.Yule(n1, t, l[1])) - sum(.PrNt.Yule(n2, t, l[2]))
- minusloglika <- function(l)
- -sum(log(.PrNt.Yule(n1, t, l[1]))) - sum(log(.PrNt.Yule(n2, t, l[2])))
+ ## initial values (moment estimators):
+ ipa <- c(mean(log(n1)/t), mean(log(n2)/t))
+ ip0 <- mean(ipa)
- out0 <- nlminb(0.01, minusloglik0, lower = 0, upper = 1)
- outa <- nlminb(rep(out0$par, 2), minusloglika,
- lower = c(0, 0), upper = c(1, 1))
+ out0 <- nlminb(ip0, minusloglik0, lower = 0, upper = 1)
+ outa <- nlminb(ipa, minusloglika, lower = c(0, 0), upper = c(1, 1))
chi <- 2 * (out0$objective - outa$objective)
pval <- pchisq(chi, 1, lower.tail = FALSE)
- data.frame("chisq" = chi, "df" = 1, "P.val" = pval, row.names = "")
+ data.frame(chisq = chi, df = 1, P.val = pval, row.names = "")
}
diversity.contrast.test <-
diff --git a/R/ace.R b/R/ace.R
index 54a2be6..6da4b39 100644
--- a/R/ace.R
+++ b/R/ace.R
@@ -1,8 +1,8 @@
-## ace.R (2015-05-01)
+## ace.R (2016-06-08)
## Ancestral Character Estimation
-## Copyright 2005-2015 Emmanuel Paradis and Ben Bolker
+## Copyright 2005-2016 Emmanuel Paradis and Ben Bolker
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
@@ -319,7 +319,7 @@ AIC.ace <- function(object, ..., k = 2)
anova.ace <- function(object, ...)
{
X <- c(list(object), list(...))
- df <- sapply(lapply(X, "[[", "rates"), length)
+ df <- lengths(lapply(X, "[[", "rates"))
ll <- sapply(X, "[[", "loglik")
## check if models are in correct order
dev <- c(NA, 2*diff(ll))
diff --git a/R/as.phylo.R b/R/as.phylo.R
index f02bb5d..63b6fee 100644
--- a/R/as.phylo.R
+++ b/R/as.phylo.R
@@ -87,7 +87,7 @@ as.phylo.phylog <- function(x, ...)
as.hclust.phylo <- function(x, ...)
{
if (!is.ultrametric(x)) stop("the tree is not ultrametric")
- if (!is.binary.tree(x)) stop("the tree is not binary")
+ if (!is.binary.phylo(x)) stop("the tree is not binary")
if (!is.rooted(x)) stop("the tree is not rooted")
n <- length(x$tip.label)
x$node.label <- NULL # by Jinlong Zhang (2010-12-15)
diff --git a/R/checkValidPhylo.R b/R/checkValidPhylo.R
index d86d3a9..8a22977 100644
--- a/R/checkValidPhylo.R
+++ b/R/checkValidPhylo.R
@@ -1,8 +1,8 @@
-## checkValidPhylo.R (2015-11-17)
+## checkValidPhylo.R (2016-07-26)
## Check the Structure of a "phylo" Object
-## Copyright 2015 Emmanuel Paradis
+## Copyright 2015-2016 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
@@ -76,8 +76,9 @@ checkValidPhylo <- function(phy)
cat(" FATAL: each tip must appear once in 'edge'\n")
if (any(tab[n + 1:m] < 2))
cat(" FATAL: all nodes should appear at least twice in 'edge'\n")
- if (any(tab[n + 2:m] < 3))
- cat(" MODERATE: only the root can be of degree 2, all other nodes should be of degree 3 or higher\n")
+ if (m > 1)
+ if (any(tab[n + 2:m] < 3))
+ cat(" MODERATE: only the root can be of degree 2, all other nodes should be of degree 3 or higher\n")
if (any(phy$edge[, 1] <= n & phy$edge[, 1] > 0))
cat(" FATAL: tips should not appear in the 1st column of 'edge'\n")
if (any(table(phy$edge[, 2]) > 1))
diff --git a/R/chronoMPL.R b/R/chronoMPL.R
index 211a8d5..c859ee8 100644
--- a/R/chronoMPL.R
+++ b/R/chronoMPL.R
@@ -9,7 +9,7 @@
chronoMPL <- function(phy, se = TRUE, test = TRUE)
{
- if (!is.binary.tree(phy)) stop("the tree is not dichotomous.")
+ if (!is.binary.phylo(phy)) stop("the tree is not dichotomous.")
n <- length(phy$tip.label)
m <- phy$Nnode
N <- dim(phy$edge)[1]
diff --git a/R/clustal.R b/R/clustal.R
index 02a9a0a..67a2b41 100644
--- a/R/clustal.R
+++ b/R/clustal.R
@@ -1,4 +1,4 @@
-## clustal.R (2016-01-26)
+## clustal.R (2016-11-29)
## Multiple Sequence Alignment with External Applications
@@ -7,6 +7,15 @@
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
+.errorAlignment <- function(exec, prog)
+{
+ dirs <- strsplit(Sys.getenv("PATH"), .Platform$path.sep)[[1]]
+ paste0("\n cannot find executable ", sQuote(exec), " on your computer.\n",
+ " It is recommended that you place the executable of ", prog, "\n",
+ " in a directory on the PATH of your computer which is:\n",
+ paste(sort(dirs), collapse = "\n"))
+}
+
clustalomega <- function(x, exec = NULL, MoreArgs = "", quiet = TRUE)
{
os <- Sys.info()[1]
@@ -16,7 +25,8 @@ clustalomega <- function(x, exec = NULL, MoreArgs = "", quiet = TRUE)
}
if (missing(x)) {
- system(paste(exec, "-h"))
+ out <- system(paste(exec, "-h"))
+ if (out == 127) stop(.errorAlignment(exec, "Clustal-Omega"))
return(invisible(NULL))
}
@@ -31,7 +41,8 @@ clustalomega <- function(x, exec = NULL, MoreArgs = "", quiet = TRUE)
opts <- paste("-i", inf, "-o", outf, "--force")
if (!quiet) opts <- paste(opts, "-v")
opts <- paste(opts, MoreArgs)
- system(paste(exec, opts), ignore.stdout = quiet)
+ out <- system(paste(exec, opts), ignore.stdout = quiet)
+ if (out == 127) stop(.errorAlignment(exec, "Clustal-Omega"))
res <- read.dna(outf, "fasta")
rownames(res) <- labels.bak
res
@@ -48,7 +59,8 @@ clustal <- function(x, pw.gapopen = 10, pw.gapext = 0.1,
}
if (missing(x)) {
- system(paste(exec, "-help"))
+ out <- system(paste(exec, "-help"))
+ if (out == 127) stop(.errorAlignment(exec, "Clustal"))
return(invisible(NULL))
}
@@ -64,7 +76,8 @@ clustal <- function(x, pw.gapopen = 10, pw.gapext = 0.1,
suffix <- c(inf, pw.gapopen, pw.gapext, gapopen, gapext)
opts <- paste(prefix, suffix, sep = "=", collapse = " ")
opts <- paste(opts, MoreArgs)
- system(paste(exec, opts), ignore.stdout = quiet)
+ out <- system(paste(exec, opts), ignore.stdout = quiet)
+ if (out == 127) stop(.errorAlignment(exec, "Clustal"))
res <- read.dna(outf, "clustal")
if (original.ordering) res <- res[labels(x), ]
rownames(res) <- labels.bak
@@ -74,7 +87,8 @@ clustal <- function(x, pw.gapopen = 10, pw.gapext = 0.1,
muscle <- function(x, exec = "muscle", MoreArgs = "", quiet = TRUE, original.ordering = TRUE)
{
if (missing(x)) {
- system(exec)
+ out <- system(exec)
+ if (out == 127) stop(.errorAlignment(exec, "MUSCLE"))
return(invisible(NULL))
}
@@ -89,7 +103,8 @@ muscle <- function(x, exec = "muscle", MoreArgs = "", quiet = TRUE, original.ord
opts <- paste("-in", inf, "-out", outf)
if (quiet) opts <- paste(opts, "-quiet")
opts <- paste(opts, MoreArgs)
- system(paste(exec, opts))
+ out <- system(paste(exec, opts))
+ if (out == 127) stop(.errorAlignment(exec, "MUSCLE"))
res <- read.dna(outf, "fasta")
if (original.ordering) res <- res[labels(x), ]
rownames(res) <- labels.bak
@@ -99,7 +114,8 @@ muscle <- function(x, exec = "muscle", MoreArgs = "", quiet = TRUE, original.ord
tcoffee <- function(x, exec = "t_coffee", MoreArgs = "", quiet = TRUE, original.ordering = TRUE)
{
if (missing(x)) {
- system(exec)
+ out <- system(exec)
+ if (out == 127) stop(.errorAlignment(exec, "T-Coffee"))
return(invisible(NULL))
}
@@ -114,7 +130,8 @@ tcoffee <- function(x, exec = "t_coffee", MoreArgs = "", quiet = TRUE, original.
write.dna(x, inf, "fasta")
opts <- paste(inf, MoreArgs)
if (quiet) opts <- paste(opts, "-quiet=nothing")
- system(paste(exec, opts))
+ out <- system(paste(exec, opts))
+ if (out == 127) stop(.errorAlignment(exec, "T-Coffee"))
res <- read.dna("input_tcoffee.aln", "clustal")
if (original.ordering) res <- res[labels(x), ]
rownames(res) <- labels.bak
diff --git a/R/coalescent.intervals.R b/R/coalescent.intervals.R
index d73bb76..1d61af3 100644
--- a/R/coalescent.intervals.R
+++ b/R/coalescent.intervals.R
@@ -15,7 +15,7 @@ coalescent.intervals.phylo <- function(x)
if (class(x) != "phylo") stop("object \"x\" is not of class \"phylo\"")
# ensure we have a BINARY tree
- if (!is.binary.tree(x)) stop("object \"x\" is not a binary tree")
+ if (!is.binary.phylo(x)) stop("object \"x\" is not a binary tree")
# ordered branching times
t <- sort(branching.times(x))
lt <- length(t)
diff --git a/R/dist.topo.R b/R/dist.topo.R
index b9c8c6f..4e7f0af 100644
--- a/R/dist.topo.R
+++ b/R/dist.topo.R
@@ -1,4 +1,4 @@
-## dist.topo.R (2016-03-23)
+## dist.topo.R (2016-11-24)
## Topological Distances, Tree Bipartitions,
## Consensus Trees, and Bootstrapping Phylogenies
@@ -42,7 +42,13 @@ dist.topo <- function(x, y = NULL, method = "PH85")
.dist.topo <- function(x, y, method = "PH85")
{
- if (method == "score" && (is.null(x$edge.length) || is.null(y$edge.length)))
+ if (method == "PH85") {
+ bs <- bitsplits(unroot.multiPhylo(c(x, y)))
+ return(sum(bs$freq == 1))
+ }
+
+ ## method == "score":
+ if (is.null(x$edge.length) || is.null(y$edge.length))
stop("trees must have branch lengths for branch score distance.")
nx <- length(x$tip.label)
x <- unroot(x)
@@ -58,48 +64,33 @@ dist.topo <- function(x, y = NULL, method = "PH85")
## End
q1 <- length(bp1)
q2 <- length(bp2)
- if (method == "PH85") {
- p <- 0
- for (i in 1:q1) {
- for (j in 1:q2) {
- if (identical(bp1[[i]], bp2[[j]]) | identical(bp1[[i]], bp2.comp[[j]])) {
- p <- p + 1
- break
- }
- }
- }
- dT <- q1 + q2 - 2 * p # same than:
- ##dT <- if (q1 == q2) 2*(q1 - p) else 2*(min(q1, q2) - p) + abs(q1 - q2)
- }
- if (method == "score") {
- dT <- 0
- found1 <- FALSE
- found2 <- logical(q2)
- found2[1] <- TRUE
- for (i in 2:q1) {
- for (j in 2:q2) {
- if (identical(bp1[[i]], bp2[[j]]) | identical(bp1[[i]], bp2.comp[[j]])) {
- dT <- dT + (x$edge.length[which(x$edge[, 2] == nx + i)] -
- y$edge.length[which(y$edge[, 2] == ny + j)])^2
- found1 <- found2[j] <- TRUE
- break
- }
+
+ dT <- 0
+ found1 <- FALSE
+ found2 <- logical(q2)
+ found2[1] <- TRUE
+ for (i in 2:q1) {
+ for (j in 2:q2) {
+ if (identical(bp1[[i]], bp2[[j]]) | identical(bp1[[i]], bp2.comp[[j]])) {
+ dT <- dT + (x$edge.length[which(x$edge[, 2] == nx + i)] -
+ y$edge.length[which(y$edge[, 2] == ny + j)])^2
+ found1 <- found2[j] <- TRUE
+ break
}
- if (found1) found1 <- FALSE
- else dT <- dT + (x$edge.length[which(x$edge[, 2] == nx + i)])^2
}
- if (!all(found2))
- dT <- dT + sum((y$edge.length[y$edge[, 2] %in% (ny + which(!found2))])^2)
- dT <- sqrt(dT)
+ if (found1) found1 <- FALSE
+ else dT <- dT + (x$edge.length[which(x$edge[, 2] == nx + i)])^2
}
- dT
+ if (!all(found2))
+ dT <- dT + sum((y$edge.length[y$edge[, 2] %in% (ny + which(!found2))])^2)
+ sqrt(dT)
}
-.compressTipLabel <- function(x)
+.compressTipLabel <- function(x, ref = NULL)
{
## 'x' is a list of objects of class "phylo" possibly with no class
if (!is.null(attr(x, "TipLabel"))) return(x)
- ref <- x[[1]]$tip.label
+ if (is.null(ref)) ref <- x[[1]]$tip.label
n <- length(ref)
if (length(unique(ref)) != n)
stop("some tip labels are duplicated in tree no. 1")
@@ -209,40 +200,63 @@ prop.clades <- function(phy, ..., part = NULL, rooted = FALSE)
n
}
-boot.phylo <- function(phy, x, FUN, B = 100, block = 1,
- trees = FALSE, quiet = FALSE, rooted = is.rooted(phy))
+boot.phylo <-
+ function(phy, x, FUN, B = 100, block = 1,
+ trees = FALSE, quiet = FALSE,
+ rooted = is.rooted(phy), jumble = TRUE,
+ mc.cores = 1)
{
if (is.null(dim(x)) || length(dim(x)) != 2)
stop("the data 'x' must have two dimensions (e.g., a matrix or a data frame)")
- boot.tree <- vector("list", B)
-
- if (!quiet) # suggestion by Alastair Potts
- progbar <- txtProgressBar(style = 3)
+ if (anyDuplicated(rownames(x)))
+ stop("some labels are duplicated in the data: you won't be able to analyse tree bipartitions")
+ boot.tree <- vector("list", B)
y <- nc <- ncol(x)
+ nr <- nrow(x)
if (block > 1) {
a <- seq(1, nc - 1, block)
b <- seq(block, nc, block)
y <- mapply(":", a, b, SIMPLIFY = FALSE)
+ getBootstrapIndices <- function() unlist(sample(y, replace = TRUE))
+ } else getBootstrapIndices <- function() sample.int(y, replace = TRUE)
+
+ if (!quiet) {
+ prefix <- "\rRunning bootstraps: "
+ suffix <- paste("/", B)
+ updateProgress <- function(i) cat(prefix, i, suffix)
}
- for (i in 1:B) {
- boot.samp <- unlist(sample(y, replace = TRUE))
- boot.tree[[i]] <- FUN(x[, boot.samp])
- if (!quiet) setTxtProgressBar(progbar, i/B)
+ if (mc.cores == 1) {
+ for (i in 1:B) {
+ boot.samp <- x[, getBootstrapIndices()]
+ if (jumble) boot.samp <- boot.samp[sample.int(nr), ]
+ boot.tree[[i]] <- FUN(boot.samp)
+ if (!quiet && !(i %% 100)) updateProgress(i)
+ }
+ } else {
+ if (!quiet) cat("Running parallel bootstraps...")
+ foo <- function(i) {
+ boot.samp <- x[, getBootstrapIndices()]
+ if (jumble) boot.samp <- boot.samp[sample.int(nr), ]
+ FUN(boot.samp)
+ }
+ boot.tree <- mclapply(1:B, foo, mc.cores = mc.cores)
+ if (!quiet) cat(" done.")
}
- ## for (i in 1:B) storage.mode(boot.tree[[i]]$Nnode) <- "integer"
- ## storage.mode(phy$Nnode) <- "integer"
+ if (!quiet) cat("\nCalculating bootstrap values...")
- if (!quiet) {
- close(progbar)
- cat("Calculating bootstrap values...")
+ ## sort labels after mixed them up
+ if (jumble) {
+ boot.tree <- .compressTipLabel(boot.tree, ref = phy$tip.label)
+ boot.tree <- .uncompressTipLabel(boot.tree)
+ boot.tree <- unclass(boot.tree) # otherwise countBipartitions crashes
}
- if (rooted) {
+ if (rooted) {
pp <- prop.part(boot.tree)
ans <- prop.clades(phy, part = pp, rooted = rooted)
} else {
@@ -252,11 +266,12 @@ boot.phylo <- function(phy, x, FUN, B = 100, block = 1,
ans <- c(B, ans[order(phy$edge[ints, 2])])
}
+ if (!quiet) cat(" done.\n")
+
if (trees) {
class(boot.tree) <- "multiPhylo"
ans <- list(BP = ans, trees = boot.tree)
}
- if (!quiet) cat(" done.\n")
ans
}
@@ -344,8 +359,7 @@ consensus <- function(..., p = 1, check.labels = TRUE)
if (p == 0.5) p <- 0.5000001 # avoid incompatible splits
pp <- pp[attr(pp, "number") >= p * ntree]
## Get the order of the remaining partitions by decreasing size:
- ind <- sort(unlist(lapply(pp, length)), decreasing = TRUE,
- index.return = TRUE)$ix
+ ind <- order(lengths(pp), decreasing = TRUE)
pp <- pp[ind]
n <- length(labels)
m <- length(pp)
diff --git a/R/is.binary.tree.R b/R/is.binary.tree.R
index 1bff336..437dba2 100644
--- a/R/is.binary.tree.R
+++ b/R/is.binary.tree.R
@@ -1,26 +1,29 @@
-## is.binary.tree.R (2002-09-12) [modified by EP 2005-05-31, 2005-08-18,
-## 2006-10-04, 2009-05-10]
+## is.binary.tree.R (2016-11-03)
-## Tests whether a given phylogenetic tree is binary
+## Test for Binary Tree
-## Copyright 2002 Korbinian Strimmer
+## Copyright 2016 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
+is.binary <- function(phy) UseMethod("is.binary")
+
+is.binary.phylo <- function(phy)
+ length(phy$tip.label) - phy$Nnode + is.rooted.phylo(phy) == 2
+
is.binary.tree <- function(phy)
{
- if (!inherits(phy, "phylo")) stop('object "phy" is not of class "phylo"')
- ## modified by EP so that it works without edge lengths too (2005-05-31):
- nb.tip <- length(phy$tip.label)
- nb.node <- phy$Nnode
- ## modified by EP so that it works with both rooted and unrooted
- ## trees (2005-08-18):
- if (is.rooted(phy)) {
- if (nb.tip - 1 == nb.node) return(TRUE)
- else return(FALSE)
- } else {
- if (nb.tip - 2 == nb.node) return(TRUE)
- else return(FALSE)
- }
+ ##warning("is.binary.tree() is deprecated; using is.binary() instead.\n\nis.binary.tree() will be removed soon: see ?is.binary and update your code.")
+ is.binary(phy)
+}
+
+is.binary.multiPhylo <- function(phy)
+{
+ phy <- unclass(phy)
+ n <- length(attr(phy, "TipLabel"))
+ if (n)
+ n - sapply(phy, "[[", "Nnode") + is.rooted.multiPhylo(phy) == 2
+ else
+ sapply(phy, is.binary.phylo)
}
diff --git a/R/is.ultrametric.R b/R/is.ultrametric.R
index cf72142..27a07da 100644
--- a/R/is.ultrametric.R
+++ b/R/is.ultrametric.R
@@ -1,27 +1,24 @@
-## is.ultrametric.R (2012-02-09)
+## is.ultrametric.R (2016-10-04)
## Test if a Tree is Ultrametric
-## Copyright 2003-2012 Emmanuel Paradis
+## Copyright 2003-2016 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
-is.ultrametric <- function(phy, tol = .Machine$double.eps^0.5)
+is.ultrametric <- function(phy, ...) UseMethod("is.ultrametric")
+
+## the main driver code (n = number of tips):
+.is.ultrametric_ape <- function(phy, tol, option, n)
{
- if (!inherits(phy, "phylo"))
- stop('object "phy" is not of class "phylo"')
if (is.null(phy$edge.length))
stop("the tree has no branch lengths")
-
- phy <- reorder(phy)
- n <- length(phy$tip.label)
e1 <- phy$edge[, 1]
e2 <- phy$edge[, 2]
EL <- phy$edge.length
- ## xx: vecteur donnant la distance d'un noeud
- ## ou d'un tip a partir de la racine
+ ## xx: distance from a node or a tip to the root
xx <- numeric(n + phy$Nnode)
## the following must start at the root and follow the
@@ -32,5 +29,31 @@ is.ultrametric <- function(phy, tol = .Machine$double.eps^0.5)
for (i in seq_len(length(e1)))
xx[e2[i]] <- xx[e1[i]] + EL[i]
- isTRUE(all.equal.numeric(var(xx[1:n]), 0, tolerance = tol))
+ xx.tip <- xx[1:n]
+
+ crit <- switch(option, {
+ mn <- min(xx.tip)
+ mx <- max(xx.tip)
+ (mx - mn)/mx
+ }, var(xx.tip))
+
+ isTRUE(all.equal.numeric(crit, 0, tolerance = tol))
+}
+
+is.ultrametric.phylo <- function(phy, tol = .Machine$double.eps^0.5,
+ option = 1, ...)
+{
+ phy <- reorder.phylo(phy)
+ .is.ultrametric_ape(phy, tol, option, length(phy$tip.label))
+}
+
+is.ultrametric.multiPhylo <- function(phy, tol = .Machine$double.eps^0.5,
+ option = 1, ...)
+{
+ phy <- reorder.multiPhylo(phy)
+ labs <- attr(phy, "TipLabel")
+ if (is.null(labs))
+ sapply(phy, is.ultrametric.phylo, tol = tol, option = option)
+ else
+ sapply(phy, .is.ultrametric_ape, tol = tol, option = option, n = length(labs))
}
diff --git a/R/ltt.plot.R b/R/ltt.plot.R
index 98f188d..7792dd0 100644
--- a/R/ltt.plot.R
+++ b/R/ltt.plot.R
@@ -10,7 +10,7 @@
ltt.plot.coords <- function(phy, backward = TRUE, tol = 1e-6)
{
if (is.ultrametric(phy, tol)) {
- if (is.binary.tree(phy)) {
+ if (is.binary.phylo(phy)) {
N <- numeric(phy$Nnode + 1)
N[] <- 1
} else {
@@ -21,9 +21,9 @@ ltt.plot.coords <- function(phy, backward = TRUE, tol = 1e-6)
names(bt) <- NULL
o <- order(bt, decreasing = TRUE)
time <- c(-bt[o], 0)
- if (!is.binary.tree(phy)) N <- c(1, N[o])
+ if (!is.binary.phylo(phy)) N <- c(1, N[o])
} else {
- if (!is.binary.tree(phy)) phy <- multi2di(phy)
+ if (!is.binary.phylo(phy)) phy <- multi2di(phy)
n <- Ntip(phy)
m <- phy$Nnode
ROOT <- n + 1L
diff --git a/R/makeLabel.R b/R/makeLabel.R
index fb688e9..f34dd64 100644
--- a/R/makeLabel.R
+++ b/R/makeLabel.R
@@ -1,8 +1,8 @@
-## makeLabel.R (2015-11-18)
+## makeLabel.R (2016-06-08)
## Label Management
-## Copyright 2010-2015 Emmanuel Paradis
+## Copyright 2010-2016 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
@@ -150,7 +150,7 @@ abbreviateGenus <- function(x, genus = TRUE, species = FALSE, sep = NULL)
if (genus) x <- sub(paste0("[[:lower:]]{1,}", sep), paste0(".", sep), x)
if (!species) return(x)
x <- strsplit(x, sep)
- k <- which(sapply(x, length) > 1)
+ k <- which(lengths(x, use.names = FALSE) > 1)
for (i in k)
x[[i]][2] <- paste0(substr(x[[i]][2], 1, 1), ".")
sapply(x, paste, collapse = sep)
diff --git a/R/multi2di.R b/R/multi2di.R
index fc598ca..364b0ea 100644
--- a/R/multi2di.R
+++ b/R/multi2di.R
@@ -1,19 +1,21 @@
-## multi2di.R (2014-10-17)
+## multi2di.R (2016-10-16)
-## Collapse and Resolve Multichotomies
+## Collapse or Resolve Multichotomies
-## Copyright 2005-2014 Emmanuel Paradis
+## Copyright 2005-2016 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
-multi2di <- function(phy, random = TRUE)
+multi2di <- function(phy, ...) UseMethod("multi2di")
+
+.multi2di_ape <- function(phy, random, n)
{
+ ## n: number of tips of phy
degree <- tabulate(phy$edge[, 1])
target <- which(degree > 2)
if (!length(target)) return(phy)
nb.edge <- dim(phy$edge)[1]
- n <- length(phy$tip.label)
nextnode <- n + phy$Nnode + 1L
new.edge <- edge2delete <- NULL
wbl <- FALSE
@@ -89,12 +91,31 @@ multi2di <- function(phy, random = TRUE)
phy
}
-di2multi <- function(phy, tol = 1e-8)
+multi2di.phylo <- function (phy, random = TRUE, ...)
+ .multi2di_ape(phy, random, length(phy$tip.label))
+
+multi2di.multiPhylo <- function(phy, random = TRUE, ...)
+{
+ labs <- attr(phy, "TipLabel")
+ oc <- oldClass(phy)
+ class(phy) <- NULL
+ if (is.null(labs)) phy <- lapply(phy, multi2di.phylo, random = random)
+ else {
+ phy <- lapply(phy, .multi2di_ape, random = random, n = length(labs))
+ attr(phy, "TipLabel") <- labs
+ }
+ class(phy) <- oc
+ phy
+}
+
+di2multi <- function(phy, ...) UseMethod("di2multi")
+
+.di2multi_ape <- function(phy, tol = 1e-8, ntips)
{
if (is.null(phy$edge.length)) stop("the tree has no branch length")
## We select only the internal branches which are
## significantly small:
- ind <- which(phy$edge.length < tol & phy$edge[, 2] > length(phy$tip.label))
+ ind <- which(phy$edge.length < tol & phy$edge[, 2] > ntips)
n <- length(ind)
if (!n) return(phy)
## recursive function to `propagate' node #'s in case
@@ -120,6 +141,23 @@ di2multi <- function(phy, tol = 1e-8)
for (i in which(sel))
phy$edge[i] <- phy$edge[i] - sum(node2del < phy$edge[i])
if (!is.null(phy$node.label))
- phy$node.label <- phy$node.label[-(node2del - length(phy$tip.label))]
+ phy$node.label <- phy$node.label[-(node2del - ntips)]
+ phy
+}
+
+di2multi.phylo <- function (phy, tol = 1e-08, ...)
+ .di2multi_ape(phy, tol, length(phy$tip.label))
+
+di2multi.multiPhylo <- function(phy, tol = 1e-08, ...)
+{
+ labs <- attr(phy, "TipLabel")
+ oc <- oldClass(phy)
+ class(phy) <- NULL
+ if (is.null(labs)) phy <- lapply(phy, di2multi.phylo, tol = tol)
+ else {
+ phy <- lapply(phy, .di2multi_ape, tol = tol, ntips = length(labs))
+ attr(phy, "TipLabel") <- labs
+ }
+ class(phy) <- oc
phy
}
diff --git a/R/node.dating.R b/R/node.dating.R
new file mode 100644
index 0000000..bafe3cf
--- /dev/null
+++ b/R/node.dating.R
@@ -0,0 +1,277 @@
+## node.dating.R (2016-06-21)
+## This file is part of the R-package `ape'.
+## See the file ../COPYING for licensing issues.
+
+# Copyright (c) 2016, Bradley R. Jones, BC Centre for Excellence in HIV/AIDS
+# All rights reserved.
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+# * Redistributions of source code must retain the above copyright
+# notice, this list of conditions and the following disclaimer.
+# * Redistributions in binary form must reproduce the above copyright
+# notice, this list of conditions and the following disclaimer in the
+# documentation and/or other materials provided with the distribution.
+# * Neither the name of the BC Centre for Excellence in HIV/AIDS nor the
+# names of its contributors may be used to endorse or promote products
+# derived from this software without specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+# DISCLAIMED. IN NO EVENT SHALL The BC Centre for Excellence in HIV/AIDS BE LIABLE FOR ANY
+# DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+# (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+# LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+# SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+# Estimate the mutation rate and node dates based on tip dates.
+#
+# Felsenstein, Joseph. "Evolutionary trees from DNA sequences: A maximum
+# likelihood approach." Journal of Molecular Evolution 17 (1981):368-376.
+#
+# Rambaut, Andrew. "Estimating the rate of molecular evolution: incorporating
+# non-contemporaneous sequences into maximum likelihood phylogenies."
+# Bioinformatics 16.4 (2000): 395-399.
+
+# Estimate the mutation rate of a phylogenetic tree from the tip dates using
+# linear regression. This model assumes that the tree follows a molecular
+# clock. It will produce a warning if the regression cannot reject the null
+# hypothesis.
+#
+# t: rooted tree with edge lengths equal to genetic distance
+#
+# tip.dates: vector of dates for the tips, in the same order as t$tip.label.
+# Tip dates can be censored with NA values
+#
+# p: p-value cutoff for failed regression (default=0.05)
+#
+# returns the mutation rate as a double
+estimate.mu <- function(t, node.dates, p = 0.05) {
+ # fit linear model
+ g <- glm(node.depth.edgelength(t)[1:length(node.dates)] ~ node.dates, na.action=na.omit)
+ null.g <- glm(node.depth.edgelength(t)[1:length(node.dates)] ~ 1, na.action=na.omit)
+
+ # test fit
+ if ((1 - pchisq(AIC(null.g) - AIC(g) + 2, df=1)) > p) {
+ warning(paste("Cannot reject null hypothesis (p=", (1 - pchisq(AIC(null.g) - AIC(g) + 2, df=1)), ")", sep=""))
+ }
+
+ coef(g)[[2]]
+}
+
+# Estimate the dates of the internal nodes of a phylogenetic tree.
+#
+# t: rooted tree with edge lengths equal to genetic distance
+#
+# node.dates: either a vector of dates for the tips, in the same order as
+# t$tip.label; or a vector of dates to initalize each node
+#
+# mu: mutation rate
+#
+# min.date: the minimum date that a node can have (needed for optimize()). The
+# default is -.Machine$double.xmax
+#
+# show.steps: set to print the log likelihood every show.steps. Set to 0 to
+# supress output
+#
+# opt.tol: tolerance for optimization precision. By default, the optimize()
+# function uses a tolerance of .Machine$double.eps^0.25 (see ?optimize)
+#
+# lik.tol: tolerance for likelihood comparison. estimate.dates will stop when
+# the log likelihood between successive trees is less than like.tol. If
+# 0 will stop after nsteps steps.
+#
+# nsteps: the maximum number of steps to run. If 0 will run until the log
+# likelihood between successive runs is less than lik.tol. The default
+# is 1000.
+#
+# is.binary: if the phylogentic tree is binary, setting is.binary to TRUE, will
+# run a optimization method
+#
+# If lik.tol and nsteps are both 0 then estimate.dates will only run the inital
+# step.
+#
+# returns a vector of all of the dates of the tips and nodes
+estimate.dates <- function(t, node.dates, mu = estimate.mu(t, node.dates), min.date = -.Machine$double.xmax, show.steps = 0, opt.tol = 1e-8, nsteps = 1000, lik.tol = if (nsteps <= 0) opt.tol else 0, is.binary = is.binary.phylo(t)) {
+ if (mu < 0)
+ stop(paste("mu (", mu, ") less than 0", sep=""))
+
+ # init vars
+ n.tips <- length(t$tip.label)
+ nodes <- unique(reorder(t)$edge[,1])
+ dates <- if (length(node.dates) == n.tips) {
+ c(node.dates, rep(NA, t$Nnode))
+ } else {
+ node.dates
+ }
+
+ # Don't count initial step if all values are seeded
+ iter.step <- if (any(is.na(dates))) {
+ 0
+ } else {
+ 1
+ }
+
+ children <- lapply(1:t$Nnode,
+ function(x) {
+ t$edge[,1] == x + n.tips
+ })
+ parent <- lapply(1:t$Nnode,
+ function(x) {
+ t$edge[,2] == x + n.tips
+ })
+
+ # to process children before parents
+ nodes <- c(1)
+ i <- 1
+
+ while (i <= length(nodes)) {
+ nodes <- c(nodes, t$edge[t$edge[,1] == nodes[i] + n.tips, 2] - n.tips)
+
+ i <- i + 1
+ }
+
+ nodes <- nodes[nodes > 0]
+ nodes <- rev(nodes)
+
+ # calculate likelihood
+ scale.lik <- sum(-lgamma(t$edge.length+1)+(t$edge.length+1)*log(mu))
+
+ calc.Like <- function(ch.node, ch.edge, x) {
+ tim <- ch.node - x
+
+ ch.edge*log(tim)-mu*tim
+ }
+
+ opt.fun <- function(x, ch, p, ch.edge.length, p.edge.length, use.parent=T) {
+ sum(if (!use.parent || length(dates[p]) == 0 || is.na(dates[p])) {
+ calc.Like(dates[ch], ch.edge.length, x)
+ } else {
+ calc.Like(c(dates[ch], x), c(ch.edge.length, p.edge.length), c(rep(x, length(dates[ch])), dates[p]))
+ })
+ }
+
+ solve.bin <- function(bounds, ch.times, ch.edge.length) {
+ a <- 2 * mu
+ b <- ch.edge.length[1] + ch.edge.length[2] - 2 * mu * (ch.times[1] + ch.times[2])
+ c.0 <- 2*mu*ch.times[1] * ch.times[2] - ch.times[1] * ch.edge.length[2] - ch.times[2] * ch.edge.length[1]
+
+ b ^ 2 - 4 * a * c.0 < 0
+
+ if (b ^ 2 - 4 * a * c.0 < 0) {
+ return(bounds[1 + (sum(calc.Like(ch.times, ch.edge.length, bounds[2] - opt.tol)) > sum(calc.Like(ch.times, ch.edge.length, bounds[1] + opt.tol)))])
+ }
+ else {
+ x.1 <- (-b + sqrt(b ^ 2 - 4 * a * c.0)) / (2 * a)
+ x.2 <- (-b - sqrt(b ^ 2 - 4 * a * c.0)) / (2 * a)
+
+ x <- c(bounds[1] + opt.tol, bounds[2] - opt.tol)
+ if (bounds[1] < x.1 && x.1 < bounds[2])
+ x <- c(x, x.1)
+ if (bounds[1] < x.2 && x.2 < bounds[2])
+ x <- c(x, x.2)
+
+ return(x[which.max(unlist(lapply(x, function(y) sum(calc.Like(ch.times, ch.edge.length, y)))))])
+ }
+ }
+
+ solve.cube <- function(bounds, ch.times, ch.edge.length, par.time, par.edge.length) {
+ a <- mu
+ b <- sum(ch.edge.length) + par.edge.length - mu * (sum(ch.times) + par.time)
+ c.0 <- mu * (ch.times[1] * ch.times[2] + ch.times[1] * par.time + ch.times[2] * par.time) - (ch.times[1] + ch.times[2]) * par.edge.length - (ch.times[1] + par.time) * ch.edge.length[2] - (ch.times[2] + par.time) * ch.edge.length[1]
+ d <- ch.edge.length[1] * ch.times[2] * par.time + ch.edge.length[2] * ch.times[1] * par.time + par.edge.length * ch.times[1] * ch.times[2] - mu * prod(ch.times) * par.time
+
+ delta.0 <- complex(real=b^2 - 3 * a * c.0)
+ delta.1 <- complex(real=2 * b^3 - 9 * a * b * c.0 + 27 * a^2 * d)
+ C <- ((delta.1 + sqrt(delta.1^2 - 4 * delta.0^3)) / 2)^(1/3)
+
+ x.1 <- Re(-1 / (3 * a) * (b + complex(real=1) * C + delta.0 / (complex(real=1) * C)))
+ x.2 <- Re(-1 / (3 * a) * (b + complex(real=-1/2, imaginary=sqrt(3)/2) * C + delta.0 / (complex(real=-1/2, imaginary=sqrt(3)/2) * C)))
+ x.3 <- Re(-1 / (3 * a) * (b + complex(real=-1/2, imaginary=-sqrt(3)/2) * C + delta.0 / (complex(real=-1/2, imaginary=-sqrt(3)/2) * C)))
+
+ x <- c(bounds[1] + opt.tol, bounds[2] - opt.tol)
+ if (x.1 && bounds[1] < x.1 && x.1 < bounds[2])
+ x <- c(x, x.1)
+ if (bounds[1] < x.2 && x.2 < bounds[2])
+ x <- c(x, x.2)
+ if (bounds[1] < x.3 && x.3 < bounds[2])
+ x <- c(x, x.3)
+
+ return(x[which.max(unlist(lapply(x, function(y) sum(calc.Like(c(ch.times, y), c(ch.edge.length, par.edge.length), c(y, y, par.time))))))])
+ }
+
+ estimate <- function(node) {
+ ch <- t$edge[children[[node]], 2]
+ ch.edge.length <- t$edge.length[children[[node]]]
+
+ p <- t$edge[parent[[node]], 1]
+ p.edge.length <- t$edge.length[parent[[node]]]
+
+ m <- if (length(p) == 0 || is.na(dates[p])) {
+ min.date
+ } else {
+ dates[p]
+ }
+
+ if (is.binary) {
+ if (length(dates[p]) == 0 || is.na(dates[p]))
+ solve.bin(c(m, min(dates[ch])), dates[ch], ch.edge.length)
+ else
+ solve.cube(c(m, min(dates[ch])), dates[ch], ch.edge.length, dates[p], p.edge.length)
+ }
+ else {
+ res <- optimize(opt.fun, c(m, min(dates[ch])), ch, p, ch.edge.length, p.edge.length, maximum=T)
+
+ res$maximum
+ }
+ }
+
+ # iterate to estimate dates
+ lik <- NA
+
+ repeat
+ {
+ for (n in nodes) {
+ dates[n + n.tips] <- estimate(n)
+ }
+
+ new.lik <- sum(unlist(lapply(nodes, function(node) {
+ ch <- t$edge[children[[node]], 2]
+ ch.edge.length <- t$edge.length[children[[node]]]
+
+ p <- t$edge[parent[[node]], 1]
+ p.edge.length <- t$edge.length[parent[[node]]]
+
+ opt.fun(dates[node+n.tips], ch, p, ch.edge.length, p.edge.length, use.parent=F)
+ }))) + scale.lik
+
+ if (show.steps > 0 && ((iter.step %% show.steps) == 0)) {
+ cat(paste("Step: ", iter.step, ", Likelihood: ", new.lik, "\n", sep=""))
+ }
+
+ if ((lik.tol > 0 && (!is.na(lik) && (is.infinite(lik) || is.infinite(new.lik) || new.lik - lik < lik.tol))) || (nsteps > 0 && iter.step >= nsteps) || (lik.tol <= 0 && nsteps <= 0)) {
+ if (is.infinite(lik) || is.infinite(new.lik)) {
+ warning("Likelihood infinite")
+ }
+ else if (!is.na(lik) && new.lik < lik - lik.tol) {
+ warning("Likelihood less than previous estimate")
+ }
+
+ break
+ } else {
+ lik <- new.lik
+ }
+
+ iter.step <- iter.step + 1
+ }
+
+ if (show.steps > 0) {
+ cat(paste("Step: ", iter.step, ", Likelihood: ", new.lik, "\n", sep=""))
+ }
+
+ dates
+}
+
diff --git a/R/pic.R b/R/pic.R
index e66042b..21e1c20 100644
--- a/R/pic.R
+++ b/R/pic.R
@@ -66,7 +66,7 @@ pic.ortho <- function(x, phy, var.contrasts = FALSE, intra = FALSE)
xx <- unlist(lapply(x, mean)) # 'x' in Felsenstein's paper
xx <- c(xx, numeric(m))
delta.v <- numeric(n + m)
- s <- 1/unlist(lapply(x, length))
+ s <- 1/lengths(x)
s <- c(s, numeric(m))
contrast <- var.cont <- numeric(m)
diff --git a/R/plot.phylo.R b/R/plot.phylo.R
index 8669c10..3c351e0 100644
--- a/R/plot.phylo.R
+++ b/R/plot.phylo.R
@@ -1,4 +1,4 @@
-## plot.phylo.R (2016-02-09)
+## plot.phylo.R (2016-07-01)
## Plot Phylogenies
@@ -51,7 +51,17 @@ plot.phylo <-
"unrooted", "radial"))
direction <- match.arg(direction, c("rightwards", "leftwards",
"upwards", "downwards"))
- if (is.null(x$edge.length)) use.edge.length <- FALSE
+ if (is.null(x$edge.length)) {
+ use.edge.length <- FALSE
+ } else {
+ if (use.edge.length && type != "radial") {
+ tmp <- sum(is.na(x$edge.length))
+ if (tmp) {
+ warning(paste(tmp, "branch length(s) NA(s): branch lengths ignored in the plot"))
+ use.edge.length <- FALSE
+ }
+ }
+ }
if (is.numeric(align.tip.label)) {
align.tip.label.lty <- align.tip.label
@@ -399,6 +409,13 @@ if (plot) {
if (type == "unrooted") "horizontal" else "axial"
} else match.arg(lab4ut, c("horizontal", "axial"))
+ xx.tips <- xx[1:Ntip]
+ yy.tips <- yy[1:Ntip]
+ if (label.offset) {
+ xx.tips <- xx.tips + label.offset * cos(angle)
+ yy.tips <- yy.tips + label.offset * sin(angle)
+ }
+
if (lab4ut == "horizontal") {
y.adj <- x.adj <- numeric(Ntip)
sel <- abs(angle) > 0.75 * pi
@@ -409,12 +426,10 @@ if (plot) {
y.adj[sel] <- strheight(x$tip.label)[sel] / 2
sel <- angle < -pi / 4 & angle > -0.75 * pi
y.adj[sel] <- -strheight(x$tip.label)[sel] * 0.75
- text(xx[1:Ntip] + x.adj * cex, yy[1:Ntip] + y.adj * cex,
+ text(xx.tips + x.adj * cex, yy.tips + y.adj * cex,
x$tip.label, adj = c(adj, 0), font = font,
srt = srt, cex = cex, col = tip.color)
} else { # if lab4ut == "axial"
- xx.tips <- xx[1:Ntip]
- yy.tips <- yy[1:Ntip]
if (align.tip.label) {
POL <- rect2polar(xx.tips, yy.tips)
POL$r[] <- max(POL$r)
@@ -423,10 +438,6 @@ if (plot) {
yy.tips <- REC$y
segments(xx[1:Ntip], yy[1:Ntip], xx.tips, yy.tips, lty = align.tip.label.lty)
}
- if (label.offset) {
- xx.tips <- xx.tips + label.offset * cos(angle)
- yy.tips <- yy.tips + label.offset * sin(angle)
- }
if (type == "unrooted") {
adj <- abs(angle) > pi/2
angle <- angle * 180/pi # switch to degrees
diff --git a/R/plot.phyloExtra.R b/R/plot.phyloExtra.R
index f56225b..938de98 100644
--- a/R/plot.phyloExtra.R
+++ b/R/plot.phyloExtra.R
@@ -1,4 +1,4 @@
-## plot.phyloExtra.R (2016-02-18)
+## plot.phyloExtra.R (2016-08-09)
## Extra Functions for Plotting and Annotating
@@ -28,3 +28,33 @@ drawSupportOnEdges <- function(value, ...)
i <- match(nodes, lastPP$edge[, 2])
edgelabels(value, i, ...)
}
+
+plotTreeTime <- function(phy, tip.dates, show.tip.label = FALSE, y.lim = NULL, ...)
+{
+ n <- Ntip(phy)
+ if (length(tip.dates) != n)
+ stop("number of dates does not match number of tips of the tree")
+ if (is.null(y.lim)) y.lim <- c(-n/4, n)
+ plot(phy, show.tip.label = show.tip.label, y.lim = y.lim, ...)
+ psr <- par("usr")
+
+ if (anyNA(tip.dates)) {
+ s <- which(!is.na(tip.dates))
+ tip.dates <- tip.dates[s]
+ } else s <- 1:n
+ range.dates <- range(as.numeric(tip.dates))
+
+ diff.range <- range.dates[2] - range.dates[1]
+ footrans <- function(x)
+ psr[2] * (as.numeric(x) - range.dates[1]) / diff.range
+
+ x1 <- footrans(tip.dates)
+ y1 <- psr[3]
+ lastPP <- get("last_plot.phylo", envir = .PlotPhyloEnv)
+ x2 <- lastPP$xx[s]
+ y2 <- lastPP$yy[s]
+ x1.scaled <- x1 / max(x1)
+ segments(x1, y1, x2, y2, col = rgb(x1.scaled, 0, 1 - x1.scaled, alpha = .5))
+ at <- pretty(tip.dates)
+ axis(1, at = footrans(at), labels = at)
+}
diff --git a/R/read.GenBank.R b/R/read.GenBank.R
index f224144..d2a44c2 100644
--- a/R/read.GenBank.R
+++ b/R/read.GenBank.R
@@ -1,53 +1,53 @@
-## read.GenBank.R (2014-07-03)
+## read.GenBank.R (2016-07-27)
## Read DNA Sequences from GenBank via Internet
-## Copyright 2002-2014 Emmanuel Paradis
+## Copyright 2002-2016 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
-read.GenBank <-
- function(access.nb, seq.names = access.nb, species.names = TRUE,
- gene.names = FALSE, as.character = FALSE)
+read.GenBank <- function(access.nb, seq.names = access.nb, species.names = TRUE,
+ gene.names = FALSE, as.character = FALSE)
{
N <- length(access.nb)
- ## If there are more than 400 sequences, we need to break down the
- ## requests, otherwise there is a segmentation fault.
- nrequest <- N %/% 400 + as.logical(N %% 400)
- X <- character(0)
- for (i in 1:nrequest) {
- a <- (i - 1) * 400 + 1
- b <- 400 * i
- if (i == nrequest) b <- N
- URL <- paste("http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nucleotide&id=",
- paste(access.nb[a:b], collapse = ","),
- "&rettype=gb&retmode=text", sep = "")
- X <- c(X, scan(file = URL, what = "", sep = "\n", quiet = TRUE))
+ ## if more than 400 sequences, we break down the requests
+ a <- 1L
+ b <- if (N > 400) 400L else N
+ fl <- tempfile()
+ repeat {
+ URL <- paste0("https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nucleotide&id=",
+ paste(access.nb[a:b], collapse = ","), "&rettype=fasta&retmode=text")
+ X <- scan(file = URL, what = "", sep = "\n", quiet = TRUE)
+ cat(X, sep = "\n", file = fl, append = TRUE)
+ if (b == N) break
+ a <- b + 1L
+ b <- b + 400L
+ if (b > N) b <- N
}
- FI <- grep("^ORIGIN ", X) + 1 # fix by Sofia Sal Bregua (2014-04-02) + Ingo Michalak (2014-07-03)
- LA <- which(X == "//") - 1
- obj <- vector("list", N)
- for (i in 1:N) {
- ## remove all spaces and digits
- tmp <- gsub("[[:digit:] ]", "", X[FI[i]:LA[i]])
- obj[[i]] <- unlist(strsplit(tmp, NULL))
- }
- names(obj) <- seq.names
- if (!as.character) obj <- as.DNAbin(obj)
+ res <- read.FASTA(fl)
+ attr(res, "description") <- names(res)
+ names(res) <- access.nb
+
+ if (as.character) res <- as.character(res)
+
if (species.names) {
- tmp <- character(N)
- sp <- grep("ORGANISM", X)
- for (i in 1:N)
- tmp[i] <- unlist(strsplit(X[sp[i]], " +ORGANISM +"))[2]
- attr(obj, "species") <- gsub(" ", "_", tmp)
- }
- if (gene.names) {
- tmp <- character(N)
- sp <- grep(" +gene +<", X)
- for (i in 1:N)
- tmp[i] <- unlist(strsplit(X[sp[i + 1L]], " +/gene=\""))[2]
- attr(obj, "gene") <- gsub("\"$", "", tmp)
+ a <- 1L
+ b <- if (N > 400) 400L else N
+ sp <- character(0)
+ repeat {
+ URL <- paste("https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nucleotide&id=",
+ paste(access.nb[a:b], collapse = ","), "&rettype=gb&retmode=text", sep = "")
+ X <- scan(file = URL, what = "", sep = "\n", quiet = TRUE, n = -1)
+ sp <- c(sp, gsub(" +ORGANISM +", "", grep("ORGANISM", X, value = TRUE)))
+ if (b == N) break
+ a <- b + 1L
+ b <- b + 400L
+ if (b > N) b <- N
+ }
+ attr(res, "species") <- gsub(" ", "_", sp)
}
- obj
+ if (gene.names)
+ warning("you used 'gene.names = TRUE': this option is obsolete; please update your code.")
+ res
}
diff --git a/R/read.dna.R b/R/read.dna.R
index c72609b..79ce558 100644
--- a/R/read.dna.R
+++ b/R/read.dna.R
@@ -1,8 +1,8 @@
-## read.dna.R (2014-05-29)
+## read.dna.R (2016-06-21)
## Read DNA Sequences in a File
-## Copyright 2003-2014 Emmanuel Paradis
+## Copyright 2003-2016 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
@@ -14,10 +14,14 @@ read.FASTA <- function(file)
file <- tempfile()
download.file(url, file)
}
- sz <- file.info(file)$size
+ sz <- file.size(file)
x <- readBin(file, "raw", sz)
- icr <- which(x == as.raw(0x0d)) # CR
- if (length(icr)) x <- x[-icr]
+ ## if the file is larger than 2 Gb we assume that it is UNIX-encoded
+ ## and skip the search-replacement of carriage returns
+ if (sz < 1e9) {
+ icr <- which(x == as.raw(0x0d)) # CR
+ if (length(icr)) x <- x[-icr]
+ }
res <- .Call(rawStreamToDNAbin, x)
names(res) <- sub("^ +", "", names(res)) # to permit phylosim
class(res) <- "DNAbin"
@@ -114,7 +118,7 @@ read.dna <- function(file, format = "interleaved", skip = 0,
rownames(obj) <- taxa
if (!as.character) obj <- as.DNAbin(obj)
} else {
- LENGTHS <- unique(unlist(lapply(obj, length)))
+ LENGTHS <- unique(lengths(obj, use.names = FALSE))
allSameLength <- length(LENGTHS) == 1
if (is.logical(as.matrix)) {
if (as.matrix && !allSameLength)
diff --git a/R/read.gff.R b/R/read.gff.R
new file mode 100644
index 0000000..6490719
--- /dev/null
+++ b/R/read.gff.R
@@ -0,0 +1,22 @@
+## read.gff.R (2016-10-06)
+
+## Read GFF Files
+
+## Copyright 2016 Emmanuel Paradis
+
+## This file is part of the R-package `ape'.
+## See the file ../COPYING for licensing issues.
+
+read.gff <- function(file, na.strings = c(".", "?"))
+{
+ w <- list("", "", "", 0L, 0L, 0, "", "", "")
+ x <- scan(file, w, sep = "\t", quote = "", quiet = TRUE,
+ na.strings = na.strings, comment.char = "#")
+ for (i in c(1, 2, 3, 7, 8)) x[[i]] <- factor(x[[i]])
+ names(x) <- c("seqid", "source", "type", "start", "end",
+ "score", "strand", "phase", "attributes")
+ n <- length(x[[1]])
+ attr(x, "row.names") <- as.character(seq_len(n))
+ class(x) <- "data.frame"
+ x
+}
diff --git a/R/reconstruct.R b/R/reconstruct.R
index 50055ea..dad9e9d 100644
--- a/R/reconstruct.R
+++ b/R/reconstruct.R
@@ -1,8 +1,8 @@
-## reconstruct.R (2014-10-24)
+## reconstruct.R (2016-07-15)
## Ancestral Character Estimation
-## Copyright 2014 Manuela Royer-Carenzi, Didier Gilles
+## Copyright 2014-2016 Manuela Royer-Carenzi, Didier Gilles
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
@@ -12,6 +12,8 @@ racine <- function(arbre) {
Ntip(arbre) + 1
}
+
+
# renvoie une liste dont la premiere composante est l'arbre renumerote
# de telle sorte que l'index d'un enfant est superieur a celui de son pere,
# la seconde compopsante est la fonction de l'index initial vers le second,
@@ -55,10 +57,11 @@ renumeroteArbre <- function(arbre) {
l
}
-#calcule la matrice C selon le modele BM
+
+#calcule la matrice C selon le modele BM ou ABM
#
-calculeC <- function(arbre) {
- m <- Ntip(arbre) + Nnode(arbre)
+calculeC_ABM <- function(arbre) {
+ m <- max(arbre[["edge"]])
C <- matrix(0,nrow=m,ncol=m)
for(i in 1:(m)) {
l <- which(arbre$edge[, 2] == i)
@@ -72,6 +75,28 @@ calculeC <- function(arbre) {
t(C)
}
+#calcule la matrice C selon le modele OU ou OU*
+#
+calculeC_OU <- function(arbre, a) {
+ m <- max(arbre[["edge"]])
+ C <- matrix(0,nrow=m,ncol=m)
+ for(i in 1:(m)) {
+ l <- which(arbre$edge[, 2] == i)
+ if(length(l)>0){
+ for(j in 1:(m)) {
+ C[j,i] <- C[j, arbre$edge[l[1], 1]]*exp(-a*arbre$edge.length[l[1]])
+ }
+ }
+ C[i,i]<-1;
+ }
+ t(C)
+}
+
+#calcule la matrice C selon le modele type qui vaut ABM ou OU
+calculeC <- function(type, arbre, a) {
+ switch(type, ABM = calculeC_ABM(arbre), OU = calculeC_OU(arbre, a))
+}
+
@@ -151,49 +176,301 @@ getREMLHessian <- function(value, arbre, sigma2) {
-reconstruct <- function (x, phyInit, method = "ML", CI = TRUE) {
+glsBM <- function (phy, x, CI=TRUE)
+{
+ obj <- list()
+ nb.tip <- length(phy$tip.label)
+ nb.node <- phy$Nnode
+ nbTotN <- nb.tip+nb.node
+ sigmaMF <- 1
+ TsTemps <- dist.nodes(phy)
+ TempsRacine <- TsTemps[(nb.tip+1),]
+ IndicesMRCA <- mrca(phy, full=T)
+ M <- matrix(NA, ncol=nbTotN, nrow=nbTotN)
+ for (i in 1:nbTotN)
+ {
+ for (j in 1:nbTotN)
+ {
+ M[i,j] <- sigmaMF^2 * TempsRacine[IndicesMRCA[i,j]]
+ }
+ }
+# M = SigmaZ
+ varAL <- M[-(1:nb.tip), 1:nb.tip]
+ varAA <- M[-(1:nb.tip), -(1:nb.tip)]
+ varLL <- M[(1:nb.tip), 1:nb.tip]
+ invVarLL <- solve(varLL)
+ UL <- rep(1, length=nb.tip)
+ UA <- rep(1, length=nb.node)
+ TL <- TempsRacine[1:nb.tip]
+ TA <- TempsRacine[(nb.tip+1):(nb.tip+nb.node)]
+#
+ IVL_Z <- invVarLL %*% x
+ IVL_T <- invVarLL %*% TL
+ IVL_U <- invVarLL %*% UL
+ U_IVL_U <- t(UL) %*% IVL_U
+ U_IVL_Z <- t(UL) %*% IVL_Z
+ DeltaU <- UA - varAL %*% IVL_U
+#
+ Racine_chap <- solve(U_IVL_U) %*% U_IVL_Z
+ Racine_chap <- as.numeric(Racine_chap)
+ Anc_chap <- Racine_chap * DeltaU + varAL %*% IVL_Z
+ Anc_chap <- as.vector(Anc_chap)
+ obj$ace <- Anc_chap
+ names(obj$ace) <- (nb.tip + 1):(nb.tip + nb.node)
+#
+ if (CI) {
+ Vec <- x - Racine_chap
+ Num <- t(Vec) %*% invVarLL %*% Vec
+ Num <- as.numeric(Num)
+ sigma2_chap <- Num / (nb.tip-1)
+ obj$sigma2 <- sigma2_chap
+ se <- sqrt((varAA - varAL %*% invVarLL %*% t(varAL))[cbind(1:nb.node,
+ 1:nb.node)])
+ se <- se * sqrt(sigma2_chap)
+ tmp <- se * qt(0.025, df=nb.tip-1)
+ obj$CI95 <- cbind(lower=obj$ace + tmp, upper=obj$ace - tmp)
+ }
+ obj
+}
+
+glsABM <- function (phy, x, CI=TRUE)
+{
+ obj <- list()
+ nb.tip <- length(phy$tip.label)
+ nb.node <- phy$Nnode
+ nbTotN <- nb.tip+nb.node
+ sigmaMF <- 1
+ TsTemps <- dist.nodes(phy)
+ TempsRacine <- TsTemps[(nb.tip+1),]
+ IndicesMRCA <- mrca(phy, full=T)
+ M <- matrix(NA, ncol=nbTotN, nrow=nbTotN)
+ for (i in 1:nbTotN)
+ {
+ for (j in 1:nbTotN)
+ {
+ M[i,j] <- sigmaMF^2 * TempsRacine[IndicesMRCA[i,j]]
+ }
+ }
+# M = SigmaZ
+ varAL <- M[-(1:nb.tip), 1:nb.tip]
+ varAA <- M[-(1:nb.tip), -(1:nb.tip)]
+ varLL <- M[(1:nb.tip), 1:nb.tip]
+ invVarLL <- solve(varLL)
+ UL <- rep(1, length=nb.tip)
+ UA <- rep(1, length=nb.node)
+ TL <- TempsRacine[1:nb.tip]
+ TA <- TempsRacine[(nb.tip+1):(nb.tip+nb.node)]
+#
+ IVL_Z <- invVarLL %*% x
+ IVL_T <- invVarLL %*% TL
+ IVL_U <- invVarLL %*% UL
+ U_IVL_U <- t(UL) %*% IVL_U
+ T_IVL_T <- t(TL) %*% IVL_T
+ U_IVL_T <- t(UL) %*% IVL_T
+ U_IVL_Z <- t(UL) %*% IVL_Z
+ T_IVL_Z <- t(TL) %*% IVL_Z
+ DeltaT <- TA - varAL %*% IVL_T
+ DeltaU <- UA - varAL %*% IVL_U
+#
+ Den <- U_IVL_U * T_IVL_T - U_IVL_T^2
+ Den <- as.numeric(Den)
+ Mu_chap <- (U_IVL_U * T_IVL_Z - U_IVL_T * U_IVL_Z) / Den
+ Mu_chap <- as.numeric(Mu_chap)
+ Racine_chap <- (T_IVL_T * U_IVL_Z - U_IVL_T * T_IVL_Z) / Den
+ Racine_chap <- as.numeric(Racine_chap)
+ Anc_chap <- Mu_chap * DeltaT + Racine_chap * DeltaU + varAL %*% IVL_Z
+ Anc_chap <- as.vector(Anc_chap)
+ obj$ace <- Anc_chap
+ names(obj$ace) <- (nb.tip + 1):(nb.tip + nb.node)
+ obj$mu <- Mu_chap
+#
+ if (CI) {
+ Vec <- x - Racine_chap - Mu_chap * TL
+ Num <- t(Vec) %*% invVarLL %*% Vec
+ Num <- as.numeric(Num)
+ sigma2_chap <- Num / (nb.tip-2)
+ obj$sigma2 <- sigma2_chap
+ se <- sqrt((varAA - varAL %*% invVarLL %*% t(varAL))[cbind(1:nb.node,
+ 1:nb.node)])
+ se <- se * sqrt(sigma2_chap)
+ tmp <- se * qt(0.025, df=nb.tip-2)
+ obj$CI95 <- cbind(lower=obj$ace + tmp, upper=obj$ace - tmp)
+ }
+ obj
+}
+
+# theta = z0
+glsOUv1 <- function (phy, x, alpha, CI=TRUE)
+{
+ obj <- list()
+ nb.tip <- length(phy$tip.label)
+ nb.node <- phy$Nnode
+ nbTotN <- nb.tip+nb.node
+ sigmaMF <- 1
+ alphaM <- alpha
+ nbTotN <- nb.tip+nb.node
+ TsTemps <- dist.nodes(phy)
+ TempsRacine <- TsTemps[(nb.tip+1),]
+ IndicesMRCA <- mrca(phy, full=T)
+ M <- matrix(NA, ncol=nbTotN, nrow=nbTotN)
+ for (i in 1:nbTotN)
+ {
+ for (j in 1:nbTotN)
+ {
+ Tempsm <- TempsRacine[IndicesMRCA[i,j]]
+ Tempsi <- TempsRacine[i]
+ Tempsj <- TempsRacine[j]
+ M[i,j] <- sigmaMF^2 * exp(-alphaM * (Tempsi+Tempsj-2*Tempsm)) * (1-exp(-2*alphaM * Tempsm)) / (2 * alphaM)
+ }
+ }
+# M = SigmaZ
+ varAL <- M[-(1:nb.tip), 1:nb.tip]
+ varAA <- M[-(1:nb.tip), -(1:nb.tip)]
+ varLL <- M[(1:nb.tip), 1:nb.tip]
+ invVarLL <- solve(varLL)
+ UL <- rep(1, length=nb.tip)
+ UA <- rep(1, length=nb.node)
+ TL <- TempsRacine[1:nb.tip]
+ TA <- TempsRacine[(nb.tip+1):(nb.tip+nb.node)]
+#
+ IVL_Z <- invVarLL %*% x
+ IVL_T <- invVarLL %*% TL
+ IVL_U <- invVarLL %*% UL
+ U_IVL_U <- t(UL) %*% IVL_U
+ U_IVL_Z <- t(UL) %*% IVL_Z
+ DeltaU <- UA - varAL %*% IVL_U
+#
+ Racine_chap <- solve(U_IVL_U) %*% U_IVL_Z
+ Racine_chap <- as.numeric(Racine_chap)
+ Anc_chap <- Racine_chap * DeltaU + varAL %*% IVL_Z
+ Anc_chap <- as.vector(Anc_chap)
+ obj$ace <- Anc_chap
+ names(obj$ace) <- (nb.tip + 1):(nb.tip + nb.node)
+#
+# vraisemblance
+#
+ mL <- Racine_chap
+ Num <- t(x-mL) %*% invVarLL %*% (x-mL)
+ Num <- as.numeric(Num)
+ sigma2_chap <- Num / (nb.tip-1)
+ obj$sigma <- sqrt(sigma2_chap)
+ VL <- sigma2_chap * varLL
+ invVL <- invVarLL / sigma2_chap
+ LVrais <- - t(x-mL) %*% invVL %*% (x-mL) /2 - nb.tip * log(2*pi)/2 - log(det(VL))/2
+ obj$LLik <- as.numeric(LVrais)
+#
+ if (CI) {
+ se <- sqrt((varAA - varAL %*% invVarLL %*% t(varAL))[cbind(1:nb.node,
+ 1:nb.node)])
+ se <- se * sqrt(sigma2_chap)
+ tmp <- se * qt(0.025, df=nb.tip-1)
+ obj$CI95 <- cbind(lower=obj$ace + tmp, upper=obj$ace - tmp)
+ }
+ obj
+}
+
+
+# theta pas egal a z0
+glsOUv2 <- function (phy, x, alpha, CI=TRUE)
+{
+ obj <- list()
+ nb.tip <- length(phy$tip.label)
+ nb.node <- phy$Nnode
+ nbTotN <- nb.tip+nb.node
+ sigmaMF <- 1
+ nbTotN <- nb.tip+nb.node
+ TsTemps <- dist.nodes(phy)
+ TempsRacine <- TsTemps[(nb.tip+1),]
+ IndicesMRCA <- mrca(phy, full=T)
+ M <- matrix(NA, ncol=nbTotN, nrow=nbTotN)
+ for (i in 1:nbTotN)
+ {
+ for (j in 1:nbTotN)
+ {
+ Tempsm <- TempsRacine[IndicesMRCA[i,j]]
+ Tempsi <- TempsRacine[i]
+ Tempsj <- TempsRacine[j]
+ M[i,j] <- sigmaMF^2 * exp(-alpha * (Tempsi+Tempsj-2*Tempsm)) * (1-exp(-2*alpha * Tempsm)) / (2 * alpha)
+ }
+ }
+# M = SigmaZ
+ varAL <- M[-(1:nb.tip), 1:nb.tip]
+ varAA <- M[-(1:nb.tip), -(1:nb.tip)]
+ varLL <- M[(1:nb.tip), 1:nb.tip]
+ invVarLL <- solve(varLL)
+ vecW <- exp(-alpha * TempsRacine)
+ UL <- vecW[1:nb.tip]
+ UA <- vecW[(nb.tip+1):(nb.tip+nb.node)]
+ TL <- 1-UL
+ TA <- 1-UA
+#
+#
+ IVL_Z <- invVarLL %*% x
+ IVL_T <- invVarLL %*% TL
+ IVL_U <- invVarLL %*% UL
+ U_IVL_U <- t(UL) %*% IVL_U
+ T_IVL_T <- t(TL) %*% IVL_T
+ U_IVL_T <- t(UL) %*% IVL_T
+ U_IVL_Z <- t(UL) %*% IVL_Z
+ T_IVL_Z <- t(TL) %*% IVL_Z
+ DeltaT <- TA - varAL %*% IVL_T
+ DeltaU <- UA - varAL %*% IVL_U
+#
+ Den <- U_IVL_U * T_IVL_T - U_IVL_T^2
+ Den <- as.numeric(Den)
+ Theta_chap <- (U_IVL_U * T_IVL_Z - U_IVL_T * U_IVL_Z) / Den
+ Theta_chap <- as.numeric(Theta_chap)
+ Racine_chap <- (T_IVL_T * U_IVL_Z - U_IVL_T * T_IVL_Z) / Den
+ Racine_chap <- as.numeric(Racine_chap)
+ Anc_chap <- Theta_chap * DeltaT + Racine_chap * DeltaU + varAL %*% IVL_Z
+ Anc_chap <- as.vector(Anc_chap)
+ obj$ace <- Anc_chap
+ names(obj$ace) <- (nb.tip + 1):(nb.tip + nb.node)
+ obj$theta <- Theta_chap
+#
+# vraisemblance
+#
+ mL <- (Racine_chap * UL + Theta_chap * TL)
+ Num <- t(x-mL) %*% invVarLL %*% (x-mL)
+ Num <- as.numeric(Num)
+ sigma2_chap <- Num / (nb.tip-2)
+ obj$sigma <- sqrt(sigma2_chap)
+ VL <- sigma2_chap * varLL
+ invVL <- invVarLL / sigma2_chap
+ LVrais <- - t(x-mL) %*% invVL %*% (x-mL) /2 - nb.tip * log(2*pi)/2 - log(det(VL))/2
+ obj$LLik <- as.numeric(LVrais)
+#
+ if (CI) {
+ se <- sqrt((varAA - varAL %*% invVarLL %*% t(varAL))[cbind(1:nb.node,
+ 1:nb.node)])
+ se <- se * sqrt(sigma2_chap)
+ tmp <- se * qt(0.025, df=nb.tip-2)
+ obj$CI95 <- cbind(lower=obj$ace + tmp, upper=obj$ace - tmp)
+ }
+ obj
+}
+
+
+reconstruct <- function (x, phyInit, method = "ML", alpha = NULL, CI = TRUE) {
+ if(!is.null(alpha)) {
+ if(alpha<=0)
+ stop("alpha has to be positive.")
+ }
if (!inherits(phyInit, "phylo"))
stop("object \"phy\" is not of class \"phylo\"")
if (is.null(phyInit$edge.length))
stop("tree has no branch lengths")
nN <- phyInit$Nnode
nT <- length(x)
-#renumber tree
- transf <- renumeroteArbre(phyInit)
- phy <- transf$arbre
- vY <- x
- for (iF in 1:nT) {
- indT <- transf$cod[iF] - nN
- vY[indT] <- x[iF]
- names(vY)[indT] <- names(x)[iF]
- }
- LongBranche <- phy$edge.length
- Fils <- phy$edge[,2]
- Tau <- phy$edge.length[order(phy$edge[,2])]
- vJ <- rep(1, length=nT)
- sigmaAcc <- Tau
- V2Acc <- diag(sigmaAcc)
- CABM <- calculeC(phy)
- C <- CABM[2:(nN+nT), 2:(nN+nT)]
- sigmaNoeuds <- C %*% V2Acc %*% t(C)
- sigma11 <- sigmaNoeuds[(1:(nN-1)), (1:(nN-1))]
- sigma22 <- sigmaNoeuds[(nN:(nN+nT-1)) , (nN:(nN+nT-1)) ]
- sigma12 <- sigmaNoeuds[(1:(nN-1)), (nN:(nN+nT-1)) ]
- sigma21 <- sigmaNoeuds[(nN:(nN+nT-1)) , (1:(nN-1))]
- GM <- solve( t(vJ) %*% solve(sigma22) %*% vJ ) %*% t(vJ) %*% solve(sigma22) %*% vY
- ZA <- rep(GM, length=nN-1)+ sigma12 %*% solve(sigma22) %*% (vY-rep(GM, length=nT))
- Intern <- c(GM, ZA)
- ValueTmp <- c(Intern, vY)
- Value <- ValueTmp
- for(iF in 1:(nN+nT)) {
- Value[transf$dec[iF]] <- ValueTmp[iF]
- }
switch(method,
ML = {
+ Intern <- glsBM(phy=phyInit, x=x, CI=F)$ace
+ Value <- c(x, Intern)
Hessian <- getMLHessian(Value, phyInit)
se <- sqrt(diag(solve(Hessian)))
se <- se[-1]
tmp <- se*qt(0.025, nN)
+ CI95 <- cbind(lower=Intern+tmp, upper=Intern-tmp)
},
REML={
minusLogLik <- function(sig2) {
@@ -203,38 +480,66 @@ reconstruct <- function (x, phyInit, method = "ML", CI = TRUE) {
logdet <- sum(log(eigen(V, symmetric = TRUE, only.values = TRUE)$values))
(nT * log(2 * pi) + logdet + distval)/2
}
+ Intern <- glsBM(phy=phyInit, x=x, CI=F)$ace
+ Value <- c(x, Intern)
+ GM <- Intern[1]
mu <- rep(GM, nT)
out <- nlm(minusLogLik, 1, hessian = FALSE)
sigma2 <- out$estimate
Hessian <- getREMLHessian(Value, phyInit, sigma2)
se <- sqrt(diag(solve(Hessian)))
tmp <- se*qt(0.025, nN)
+ CI95 <- cbind(lower=Intern+tmp, upper=Intern-tmp)
},
GLS = {
- Hessian <- (sigma11 - sigma12 %*% solve(sigma22) %*% t(sigma12))
- seTmp <- c(0, sqrt(diag(Hessian)))
- se <- seTmp
- for(iF in 1:nN) {
- se[iF] <- seTmp[transf$cod[iF+nT]]
- }
- tmp <- se*qnorm(0.025)
+ result <- glsBM(phy=phyInit, x=x, CI=T)
+ Intern <- result$ace
+ CI95 <- result$CI95
+ },
+ GLS_ABM = {
+ result <- glsABM(phy=phyInit, x=x, CI=T)
+ Intern <- result$ace
+ CI95 <- result$CI95
+ },
+ GLS_OUS = {
+ if(is.null(alpha)) {
+ funOpt1 <- function(alpha)
+ {
+ -glsOUv1(phy=phyInit, x=x, alpha, CI=F)$LLik
+ }
+ calOp <- optim(par=0.25, fn=funOpt1, method="Brent", lower=0.001, upper=1)
+ if (calOp$convergence == 0)
+ {
+ alpha <- calOp$par
+ } else {
+ stop("Estimation error for alpha")
+ }
+ }
+ result <- glsOUv1(phy=phyInit, x=x, alpha=alpha, CI=T)
+ Intern <- result$ace
+ CI95 <- result$CI95
+ },
+ GLS_OU = {
+ if(is.null(alpha)) {
+ funOpt2 <- function(alpha)
+ {
+ -glsOUv2(phy=phyInit, x=x, alpha, CI=F)$LLik
+ }
+ calOp <- optim(par=0.25, fn=funOpt2, method="Brent", lower=0.001, upper=1)
+ if (calOp$convergence == 0)
+ {
+ alpha <- calOp$par
+ } else {
+ stop("Estimation error for alpha")
+ }
+ }
+ result <- glsOUv2(phy=phyInit, x=x, alpha=alpha, CI=T)
+ Intern <- result$ace
+ CI95 <- result$CI95
}
)
- InternOP <- Intern
- for (iF in 1:nN) {
- InternOP[iF] <- Value[iF+nT]
- }
- CI95 <- cbind(lower=InternOP+tmp, upper=InternOP-tmp)
if (CI==TRUE)
- list(ace=InternOP, CI95=CI95)
+ list(ace=Intern, CI95=CI95)
else
- list(ace=InternOP)
+ list(ace=Intern)
}
-
-
-
-
-
-
-
-
diff --git a/R/reorder.phylo.R b/R/reorder.phylo.R
index 9cf74ed..718aadb 100644
--- a/R/reorder.phylo.R
+++ b/R/reorder.phylo.R
@@ -1,18 +1,14 @@
-## reorder.phylo.R (2015-08-24)
+## reorder.phylo.R (2016-10-04)
## Internal Reordering of Trees
-## Copyright 2006-2015 Emmanuel Paradis
+## Copyright 2006-2016 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
-reorder.phylo <- function(x, order = "cladewise", index.only = FALSE, ...)
+.reorder_ape <- function(x, order, index.only, nb.tip, io)
{
- ORDER <- c("cladewise", "postorder", "pruningwise")
- io <- pmatch(order, ORDER)
- if (is.na(io)) stop("ambiguous order")
- order <- ORDER[io]
nb.edge <- dim(x$edge)[1]
if (!is.null(attr(x, "order")))
if (attr(x, "order") == order)
@@ -21,8 +17,6 @@ reorder.phylo <- function(x, order = "cladewise", index.only = FALSE, ...)
if (nb.node == 1)
if (index.only) return(1:nb.edge) else return(x)
- nb.tip <- length(x$tip.label)
-
## I'm adding the next check for badly conformed trees to avoid R
## crashing (2013-05-17):
if (nb.node >= nb.tip)
@@ -50,6 +44,31 @@ reorder.phylo <- function(x, order = "cladewise", index.only = FALSE, ...)
x
}
+reorder.phylo <- function(x, order = "cladewise", index.only = FALSE, ...)
+{
+ ORDER <- c("cladewise", "postorder", "pruningwise")
+ io <- pmatch(order, ORDER)
+ if (is.na(io)) stop("ambiguous order")
+ order <- ORDER[io]
+ .reorder_ape(x, order, index.only, length(x$tip.label), io)
+}
+
+reorder.multiPhylo <- function(x, order = "cladewise", ...)
+{
+ ORDER <- c("cladewise", "postorder", "pruningwise")
+ io <- pmatch(order, ORDER)
+ if (is.na(io)) stop("ambiguous order")
+ order <- ORDER[io]
+ oc <- oldClass(x)
+ class(x) <- NULL
+ labs <- attr(x, "TipLabel")
+ x <-
+ if (is.null(labs)) lapply(x, reorder.phylo, order = order)
+ else lapply(x, .reorder_ape, order = order, nb.tip = length(labs), io = io)
+ class(x) <- oc
+ x
+}
+
rotateConstr <- function(phy, constraint)
{
D <- match(phy$tip.label, constraint)
diff --git a/R/root.R b/R/root.R
index 9c89319..379ffaf 100644
--- a/R/root.R
+++ b/R/root.R
@@ -1,26 +1,36 @@
-## root.R (2016-04-05)
+## root.R (2016-11-16)
-## Root of Phylogenetic Trees
+## Roots Phylogenetic Trees
## Copyright 2004-2016 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
-is.rooted <- function(phy)
+is.rooted <- function(phy) UseMethod("is.rooted")
+
+.is.rooted_ape <- function(phy, ntips)
{
- if (!inherits(phy, "phylo"))
- stop('object "phy" is not of class "phylo"')
- if (!is.null(phy$root.edge)) TRUE
- else
- if (tabulate(phy$edge[, 1])[length(phy$tip.label) + 1] > 2)
- FALSE else TRUE
+ if (!is.null(phy$root.edge)) return(TRUE)
+ if (tabulate(phy$edge[, 1])[ntips + 1] > 2) FALSE else TRUE
}
-unroot <- function(phy)
+is.rooted.phylo <- function (phy)
+ .is.rooted_ape(phy, length(phy$tip.label))
+
+is.rooted.multiPhylo <- function(phy)
{
- if (!inherits(phy, "phylo"))
- stop('object "phy" is not of class "phylo"')
+ phy <- unclass(phy)
+ labs <- attr(phy, "TipLabel")
+ if (is.null(labs)) sapply(phy, is.rooted.phylo)
+ else sapply(phy, .is.rooted_ape, ntips = length(labs))
+}
+
+unroot <- function(phy) UseMethod("unroot")
+
+.unroot_ape <- function(phy, n)
+{
+ ## n: number of tips of phy
N <- dim(phy$edge)[1]
if (N < 3)
stop("cannot unroot a tree with less than three edges.")
@@ -28,13 +38,12 @@ unroot <- function(phy)
## delete FIRST the root.edge (in case this is sufficient to
## unroot the tree, i.e. there is a multichotomy at the root)
if (!is.null(phy$root.edge)) phy$root.edge <- NULL
- if (!is.rooted(phy)) return(phy)
+ if (!.is.rooted_ape(phy, n)) return(phy)
- n <- length(phy$tip.label)
ROOT <- n + 1L
-### EDGEROOT[1]: the edge to delete
-### EDGEROOT[2]: the target where to stick the edge to delete
+### EDGEROOT[1]: the edge to be deleted
+### EDGEROOT[2]: the target where to stick the edge to be deleted
### If the tree is in pruningwise (or postorder) order, then
### the last two edges are those connected to the root; the node
@@ -44,6 +53,9 @@ unroot <- function(phy)
if (!is.null(ophy) && ophy != "cladewise") {
NEWROOT <- phy$edge[N - 2L, 1L]
EDGEROOT <- c(N, N - 1L)
+ ## make sure EDGEROOT is ordered as described above:
+ if (phy$edge[EDGEROOT[1L], 2L] != NEWROOT)
+ EDGEROOT <- EDGEROOT[2:1]
} else {
### ... otherwise, we remove one of the edges coming from
@@ -55,13 +67,13 @@ unroot <- function(phy)
### origin of the tree).
EDGEROOT <- which(phy$edge[, 1L] == ROOT)
- NEWROOT <- ROOT + 1L
+##### NEWROOT <- ROOT + 1L
+ ## make sure EDGEROOT is ordered as described above:
+ if (phy$edge[EDGEROOT[1L], 2L] <= n)
+ EDGEROOT <- EDGEROOT[2:1]
+ NEWROOT <- phy$edge[EDGEROOT[1L], 2L]
}
- ## make sure EDGEROOT is ordered as described above:
- if (phy$edge[EDGEROOT[1L], 2L] != NEWROOT)
- EDGEROOT <- EDGEROOT[2:1]
-
phy$edge <- phy$edge[-EDGEROOT[1L], ]
s <- phy$edge == NEWROOT # renumber the new root
@@ -91,8 +103,27 @@ unroot <- function(phy)
phy
}
-root <- function(phy, outgroup, node = NULL,
- resolve.root = FALSE, interactive = FALSE)
+unroot.phylo <- function(phy)
+ .unroot_ape(phy, length(phy$tip.label))
+
+unroot.multiPhylo <- function(phy)
+{
+ oc <- oldClass(phy)
+ class(phy) <- NULL
+ labs <- attr(phy, "TipLabel")
+ if (is.null(labs)) phy <- lapply(phy, unroot.phylo)
+ else {
+ phy <- lapply(phy, .unroot_ape, n = length(labs))
+ attr(phy, "TipLabel") <- labs
+ }
+ class(phy) <- oc
+ phy
+}
+
+root <- function(phy, ...) UseMethod("root")
+
+root.phylo <- function(phy, outgroup, node = NULL, resolve.root = FALSE,
+ interactive = FALSE, edgelabel = FALSE, ...)
{
if (!inherits(phy, "phylo"))
stop('object not of class "phylo"')
@@ -201,7 +232,7 @@ root <- function(phy, outgroup, node = NULL,
if (!fuseRoot) INV[i] <- TRUE
## bind the other clades...
-# for (j in 1:Nclade) { # fix by EP (2016-01-04)
+
if (fuseRoot) { # do we have to fuse the two basal edges?
k <- which(e1 == ROOT)
k <- if (k[2] > w) k[2] else k[1]
@@ -209,11 +240,17 @@ root <- function(phy, outgroup, node = NULL,
if (wbl)
phy$edge.length[k] <- phy$edge.length[k] + phy$edge.length[i]
}
-# }
if (fuseRoot) phy$Nnode <- oldNnode - 1L
+ ## added after discussion with Jaime Huerta Cepas (2016-07-30):
+ if (edgelabel) {
+ phy$node.label[e1[INV] - n] <- phy$node.label[e2[INV] - n]
+ phy$node.label[newroot - n] <- ""
+ }
+
phy$edge[INV, ] <- phy$edge[INV, 2:1]
+
if (fuseRoot) {
phy$edge <- phy$edge[-i, ]
if (wbl) phy$edge.length <- phy$edge.length[-i]
@@ -268,3 +305,19 @@ root <- function(phy, outgroup, node = NULL,
attr(phy, "order") <- NULL
reorder.phylo(phy)
}
+
+root.multiPhylo <- function(phy, outgroup, ...)
+{
+ oc <- oldClass(phy)
+ class(phy) <- NULL
+ labs <- attr(phy, "TipLabel")
+ if (!is.null(labs))
+ for (i in seq_along(phy)) phy[[i]]$tip.label <- labs
+ phy <- lapply(phy, root.phylo, outgroup = outgroup, ...)
+ if (!is.null(labs)) {
+ for (i in seq_along(phy)) phy[[i]]$tip.label <- NULL
+ attr(phy, "TipLabel") <- labs
+ }
+ class(phy) <- oc
+ phy
+}
diff --git a/R/rtt.R b/R/rtt.R
index 14095ce..e6f186c 100644
--- a/R/rtt.R
+++ b/R/rtt.R
@@ -7,13 +7,13 @@
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
-rtt <- function (t, tip.dates, ncpu = 1, objective = "correlation",
+rtt <- function (t, tip.dates, ncpu = 1, objective = "correlation",
opt.tol = .Machine$double.eps^0.25) {
- if (objective == "correlation")
+ if (objective == "correlation")
objective <- function(x, y) cor.test(y, x)$estimate
- else if (objective == "rsquared")
+ else if (objective == "rsquared")
objective <- function(x, y) summary(lm(y ~ x))$r.squared
- else if (objective == "rms")
+ else if (objective == "rms")
objective <- function(x, y) -summary(lm(y ~ x))$sigma^2
else stop("objective must be one of \"correlation\", \"rsquared\", or \"rms\"")
@@ -25,8 +25,8 @@ rtt <- function (t, tip.dates, ncpu = 1, objective = "correlation",
objective(tip.dates, edge.dist)
}
- obj.edge <- if (ncpu > 1)
- unlist(parallel::mclapply(1:nrow(ut$edge), function (e) {
+ obj.edge <- if (ncpu > 1)
+ unlist(mclapply(1:nrow(ut$edge), function (e) {
opt.fun <- function (x) f(x, ut$edge[e,1], ut$edge[e,2])
optimize(opt.fun, c(0, 1), maximum = TRUE, tol = opt.tol)$objective
}, mc.cores=ncpu))
@@ -44,10 +44,10 @@ rtt <- function (t, tip.dates, ncpu = 1, objective = "correlation",
opt.fun <- function (x) f(x, best.edge.parent, best.edge.child)
best.pos <- optimize(opt.fun, c(0, 1), maximum = TRUE, tol = opt.tol)$maximum
- new.root <- list(edge = matrix(c(2L, 1L), 1, 2), tip.label = "new.root",
+ new.root <- list(edge = matrix(c(2L, 1L), 1, 2), tip.label = "new.root",
edge.length = 1, Nnode = 1L, root.edge = 1)
class(new.root) <- "phylo"
- ut <- bind.tree(ut, new.root, where = best.edge.child, position = best.pos *
+ ut <- bind.tree(ut, new.root, where = best.edge.child, position = best.pos *
best.edge.length)
ut <- collapse.singles(ut)
ut <- root(ut, "new.root")
diff --git a/R/summary.phylo.R b/R/summary.phylo.R
index f390b24..e8d4817 100644
--- a/R/summary.phylo.R
+++ b/R/summary.phylo.R
@@ -1,34 +1,44 @@
-## summary.phylo.R (2011-08-04)
+## summary.phylo.R (2016-11-24)
## Print Summary of a Phylogeny and "multiPhylo" operators
-## Copyright 2003-2011 Emmanuel Paradis, and 2006 Ben Bolker
+## Copyright 2003-2016 Emmanuel Paradis, 2006 Ben Bolker, and Klaus Schliep 2016
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
-Ntip <- function(phy)
+Ntip <- function(phy) UseMethod("Ntip")
+
+Ntip.phylo <- function(phy) length(phy$tip.label)
+
+Ntip.multiPhylo <- function(phy)
{
- if (!inherits(phy, "phylo"))
- stop('object "phy" is not of class "phylo"')
- length(phy$tip.label)
+ labs <- attr(phy, "TipLabel")
+ if (is.null(labs)) sapply(unclass(phy), Ntip.phylo)
+ else setNames(rep(length(labs), length(phy)), names(phy))
}
-Nnode <- function(phy, internal.only = TRUE)
+Nnode <- function(phy, ...) UseMethod("Nnode")
+
+Nnode.phylo <- function(phy, internal.only = TRUE, ...)
{
- if (!inherits(phy, "phylo"))
- stop('object "phy" is not of class "phylo"')
if (internal.only) return(phy$Nnode)
phy$Nnode + length(phy$tip.label)
}
-Nedge <- function(phy)
+Nnode.multiPhylo <- function(phy, internal.only = TRUE, ...)
{
- if (!inherits(phy, "phylo"))
- stop('object "phy" is not of class "phylo"')
- dim(phy$edge)[1]
+ res <- sapply(unclass(phy), "[[", "Nnode")
+ if (internal.only) return(res)
+ res + Ntip.multiPhylo(phy)
}
+Nedge <- function(phy) UseMethod("Nedge")
+
+Nedge.phylo <- function(phy) dim(phy$edge)[1]
+
+Nedge.multiPhylo <- function(phy) sapply(unclass(phy), Nedge.phylo)
+
summary.phylo <- function(object, ...)
{
cat("\nPhylogenetic tree:", deparse(substitute(object)), "\n\n")
@@ -134,27 +144,79 @@ str.multiPhylo <- function(object, ...)
str(object, ...)
}
-c.phylo <- function(..., recursive = FALSE)
- structure(list(...), class = "multiPhylo")
-## only the first object in '...' is checked for its class,
-## but that should be OK for the moment
+.c_phylo_single <- function(phy) structure(list(phy), class = "multiPhylo")
-c.multiPhylo <- function(..., recursive = FALSE)
+c.phylo <- function(..., recursive = TRUE)
{
obj <- list(...)
+ classes <- lapply(obj, class)
+ isphylo <- sapply(classes, function(x) "phylo" %in% x)
+ if (all(isphylo)) {
+ class(obj) <- "multiPhylo"
+ return(obj)
+ }
+ if (!recursive) return(obj)
+ ismulti <- sapply(classes, function(x) "multiPhylo" %in% x)
+ if (all(isphylo | ismulti)) {
+ for (i in which(isphylo)) obj[[i]] <- .c_phylo_single(obj[[i]])
+ ## added by Klaus:
+ for (i in which(ismulti)) obj[[i]] <- .uncompressTipLabel(obj[[i]])
+ obj <- .makeMultiPhyloFromObj(obj)
+ } else {
+ warning('some objects not of class "phylo" or "multiPhylo": argument recursive=TRUE ignored')
+ }
+ obj
+}
+
+# this is an option to avoid growing the list, better check it also
+# not really as important as long the list of trees is short (by Klaus)
+.makeMultiPhyloFromObj <- function(obj)
+{
n <- length(obj)
- x <- obj[[1L]]
- N <- length(x)
- i <- 2L
- while (i <= n) {
- a <- N + 1L
- N <- N + length(obj[[i]])
- ## x is of class "multiPhylo", so this uses the operator below:
- x[a:N] <- obj[[i]]
- i <- i + 1L
+ N <- lengths(obj, FALSE)
+ cs <- c(0, cumsum(N))
+ x <- vector("list", cs[length(cs)])
+ for (i in 1:n) {
+ a <- cs[i] + 1L
+ b <- cs[i + 1L]
+ x[a:b] <- obj[[i]]
}
+ class(x) <- "multiPhylo"
x
}
+## the original code:
+##.makeMultiPhyloFromObj <- function(obj)
+##{
+## n <- length(obj)
+## x <- obj[[1L]]
+## N <- length(x)
+## i <- 2L
+## while (i <= n) {
+## a <- N + 1L
+## N <- N + length(obj[[i]])
+## ## x is of class "multiPhylo", so this uses the operator below:
+## x[a:N] <- obj[[i]]
+## i <- i + 1L
+## }
+## x
+##}
+
+c.multiPhylo <- function(..., recursive = TRUE)
+{
+ obj <- list(...)
+ if (!recursive) return(obj)
+ classes <- lapply(obj, class)
+ isphylo <- sapply(classes, function(x) "phylo" %in% x)
+ ismulti <- sapply(classes, function(x) "multiPhylo" %in% x)
+ if (!all(isphylo | ismulti)) {
+ warning('some objects not of class "phylo" or "multiPhylo": argument recursive=TRUE ignored')
+ return(obj)
+ }
+ for (i in which(isphylo)) obj[[i]] <- .c_phylo_single(obj[[i]])
+ ## added by Klaus
+ for (i in which(ismulti)) obj[[i]] <- .uncompressTipLabel(obj[[i]])
+ .makeMultiPhyloFromObj(obj)
+}
.uncompressTipLabel <- function(x)
{
diff --git a/R/write.dna.R b/R/write.dna.R
index 3f53cb9..632c286 100644
--- a/R/write.dna.R
+++ b/R/write.dna.R
@@ -26,7 +26,7 @@ write.dna <- function(x, file, format = "interleaved", append = FALSE,
rm(xx)
} else {
N <- length(x)
- S <- unique(unlist(lapply(x, length)))
+ S <- unique(lengths(x, use.names = FALSE))
if (length(S) > 1) aligned <- FALSE
}
if (is.null(names(x))) names(x) <- as.character(1:N)
diff --git a/R/yule.R b/R/yule.R
index abd719c..edabb8c 100644
--- a/R/yule.R
+++ b/R/yule.R
@@ -12,7 +12,7 @@
yule <- function(phy, use.root.edge = FALSE)
{
- if (!is.binary.tree(phy))
+ if (!is.binary.phylo(phy))
stop("tree must be dichotomous to fit the Yule model.")
X <- sum(phy$edge.length)
diff --git a/build/vignette.rds b/build/vignette.rds
index 597d4e5..bd50957 100644
Binary files a/build/vignette.rds and b/build/vignette.rds differ
diff --git a/data/hivtree.newick.rda b/data/hivtree.newick.rda
index b731460..60c5da9 100644
Binary files a/data/hivtree.newick.rda and b/data/hivtree.newick.rda differ
diff --git a/data/landplants.newick.rda b/data/landplants.newick.rda
deleted file mode 100644
index 284b0e1..0000000
Binary files a/data/landplants.newick.rda and /dev/null differ
diff --git a/data/opsin.newick.rda b/data/opsin.newick.rda
deleted file mode 100644
index b6d33b9..0000000
Binary files a/data/opsin.newick.rda and /dev/null differ
diff --git a/inst/doc/MoranI.pdf b/inst/doc/MoranI.pdf
index 846cec6..aa82ee1 100644
Binary files a/inst/doc/MoranI.pdf and b/inst/doc/MoranI.pdf differ
diff --git a/man/AAbin.Rd b/man/AAbin.Rd
index 7f46f25..a03d95b 100644
--- a/man/AAbin.Rd
+++ b/man/AAbin.Rd
@@ -98,4 +98,5 @@ data(woodmouse)
AA <- trans(woodmouse, 2)
seg.sites(woodmouse)
AAsubst(AA)
-}
\ No newline at end of file
+}
+\keyword{manip}
diff --git a/man/all.equal.DNAbin.Rd b/man/all.equal.DNAbin.Rd
new file mode 100644
index 0000000..1bcbdf0
--- /dev/null
+++ b/man/all.equal.DNAbin.Rd
@@ -0,0 +1,65 @@
+\name{all.equal.DNAbin}
+\alias{all.equal.DNAbin}
+\title{Compare DNA Sets}
+\description{
+ Comparison of DNA sequence sets, particularly when aligned.
+}
+\usage{
+\method{all.equal}{DNAbin}(target, current, plot = FALSE, ...)
+}
+\arguments{
+ \item{target, current}{the two sets of sequences to be compared.}
+ \item{plot}{a logical value specifying whether to plot the sites that
+ are different (only if the labels of both alignments are the same).}
+ \item{\dots}{further arguments passed to \code{\link{image.DNAbin}}.}
+}
+\details{
+ If the two sets of DNA sequences are exactly identical, this function
+ returns \code{TRUE}. Otherwise, a detailed comparison is made only if
+ the labels (i.e., rownames) of \code{target} and \code{current} are the
+ same (possibly in different orders). In all other cases, a brief
+ description of the differences is returned (sometimes with
+ recommendations to make further comparisons).
+
+ This function can be used for testing in programs using
+ \code{\link[base]{isTRUE}} (see examples below).
+}
+\value{
+ \code{TRUE} if the two sets are identical; a list with two elements
+ (message and different.sites) if a detailed comparison is done; or a
+ vector of mode character.
+}
+\author{Emmanuel Paradis}
+\seealso{
+ \code{\link{image.DNAbin}}, \code{\link{clustal}},
+ \code{\link{checkAlignment}},
+ the generic function: \code{\link[base]{all.equal}}
+}
+\examples{
+data(woodmouse)
+woodm2 <- woodmouse
+woodm2[1, c(1:5, 10:12, 30:40)] <- as.DNAbin("g")
+res <- all.equal(woodmouse, woodm2, plot = TRUE)
+str(res)
+
+## if used for testing in R programs:
+isTRUE(all.equal(woodmouse, woodmouse)) # TRUE
+isTRUE(all.equal(woodmouse, woodm2)) # FALSE
+
+all.equal(woodmouse, woodmouse[15:1, ])
+all.equal(woodmouse, woodmouse[-1, ])
+all.equal(woodmouse, woodmouse[, -1])
+
+\dontrun{
+## To run the followings you need internet and Clustal and MUSCLE
+## correctly installed.
+## Data from Johnson et al. (2006, Science)
+refs <- paste("DQ082", 505:545, sep = "")
+DNA <- read.GenBank(refs)
+DNA.clustal <- clustal(DNA)
+DNA.muscle <- muscle(DNA)
+isTRUE(all.equal(DNA.clustal, DNA.muscle)) # FALSE
+all.equal(DNA.clustal, DNA.muscle, TRUE)
+}
+}
+\keyword{manip}
diff --git a/man/alview.Rd b/man/alview.Rd
index deaca18..97809d1 100644
--- a/man/alview.Rd
+++ b/man/alview.Rd
@@ -22,7 +22,7 @@ The first line of the output shows the position of the last column of the printe
\author{Emmanuel Paradis}
\seealso{
\code{\link{DNAbin}}, \code{\link{image.DNAbin}}, \code{\link{alex}},
- \code{\link{clustal}}, \code{\link{checkAlignment}}
+ \code{\link{clustal}}, \code{\link{checkAlignment}}, \code{\link{all.equal.DNAbin}}
}
\examples{
data(woodmouse)
diff --git a/man/ape-internal.Rd b/man/ape-internal.Rd
index d320282..29ce7bc 100644
--- a/man/ape-internal.Rd
+++ b/man/ape-internal.Rd
@@ -72,6 +72,12 @@
\alias{racine}
\alias{renumeroteArbre}
\alias{calculeC}
+\alias{calculeC_ABM}
+\alias{calculeC_OU}
+\alias{glsABM}
+\alias{glsBM}
+\alias{glsOUv1}
+\alias{glsOUv2}
\alias{getSumSquare}
\alias{getMLHessian}
\alias{getREMLHessian}
diff --git a/man/ape-package.Rd b/man/ape-package.Rd
index 46c7b61..996b3f8 100644
--- a/man/ape-package.Rd
+++ b/man/ape-package.Rd
@@ -34,5 +34,9 @@ Analyses of Phylogenetics and Evolution
Paradis, E., Claude, J. and Strimmer, K. (2004) APE: analyses of
phylogenetics and evolution in R language. \emph{Bioinformatics},
\bold{20}, 289--290.
+
+ Popescu, A.-A., Huber, K. T. and Paradis, E. (2012) ape 3.0: new tools
+ for distance based phylogenetics and evolutionary analysis in
+ R. \emph{Bioinformatics}, \bold{28}, 1536-1537.
}
\keyword{package}
diff --git a/man/boot.phylo.Rd b/man/boot.phylo.Rd
index d685e31..602ae8a 100644
--- a/man/boot.phylo.Rd
+++ b/man/boot.phylo.Rd
@@ -22,7 +22,9 @@
}
\usage{
boot.phylo(phy, x, FUN, B = 100, block = 1,
- trees = FALSE, quiet = FALSE, rooted = is.rooted(phy))
+ trees = FALSE, quiet = FALSE,
+ rooted = is.rooted(phy), jumble = TRUE,
+ mc.cores = 1)
prop.part(..., check.labels = TRUE)
prop.clades(phy, ..., part = NULL, rooted = FALSE)
\method{print}{prop.part}(x, ...)
@@ -43,6 +45,11 @@ prop.clades(phy, ..., part = NULL, rooted = FALSE)
\item{quiet}{a logical: a progress bar is displayed by default.}
\item{rooted}{a logical specifying whether the trees should be treated
as rooted or not.}
+ \item{jumble}{a logical value. By default, the rows of \code{x} are
+ randomized to avoid artificially too large bootstrap values
+ associated with very short branches.}
+ \item{mc.cores}{the number of cores (CPUs) to be used (passed to
+ \pkg{parallel}).}
\item{\dots}{either (i) a single object of class \code{"phylo"}, (ii) a
series of such objects separated by commas, or (iii) a list
containing such objects. In the case of \code{plot} further
diff --git a/man/c.phylo.Rd b/man/c.phylo.Rd
index 98497f3..283fffa 100644
--- a/man/c.phylo.Rd
+++ b/man/c.phylo.Rd
@@ -8,21 +8,26 @@
These functions help to build lists of trees of class \code{"multiPhylo"}.
}
\usage{
-\method{c}{phylo}(..., recursive = FALSE)
-\method{c}{multiPhylo}(..., recursive = FALSE)
-.compressTipLabel(x)
+\method{c}{phylo}(..., recursive = TRUE)
+\method{c}{multiPhylo}(..., recursive = TRUE)
+.compressTipLabel(x, ref = NULL)
.uncompressTipLabel(x)
}
\arguments{
- \item{\dots}{one or several objects of class \code{"phylo"} or
+ \item{\dots}{one or several objects of class \code{"phylo"} and/or
\code{"multiPhylo"}.}
\item{recursive}{for compatibily with the generic (do not change).}
\item{x}{an object of class \code{"phylo"} or \code{"multiPhylo"}.}
+ \item{ref}{an optional vector of mode character to constrain the order
+ of the tips. By default, the order from the first tree is used.}
}
\details{
- These \code{c} methods do not check all the arguments, so it is the
- user's responsibility to make sure that only objects of the same class
- (either \code{"phylo"} or \code{"multiPhylo"}) are used.
+ These \code{c} methods check all the arguments, and return by default
+ a list of single trees unless some objects are not trees or lists of
+ trees, in which case \code{recursive} is switched to FALSE and a
+ warning message is given. If \code{recursive = FALSE}, the objects are
+ simply concatenated into a list. Before \pkg{ape} 4.0, \code{recursive}
+ was always set to FALSE.
\code{.compressTipLabel} transforms an object of class
\code{"multiPhylo"} by checking that all trees have the same tip
@@ -37,9 +42,7 @@
An object of class \code{"multiPhylo"}.
}
\author{Emmanuel Paradis}
-\seealso{
- \code{\link{summary.phylo}}, \code{\link{multiphylo}}
-}
+\seealso{\code{\link{summary.phylo}}, \code{\link{multiphylo}}}
\examples{
x <- c(rtree(4), rtree(2))
x
diff --git a/man/checkAlignment.Rd b/man/checkAlignment.Rd
index a64fe15..bfc6a17 100644
--- a/man/checkAlignment.Rd
+++ b/man/checkAlignment.Rd
@@ -34,7 +34,7 @@ checkAlignment(x, check.gaps = TRUE, plot = TRUE, what = 1:4)
\value{NULL}
\author{Emmanuel Paradis}
\seealso{
-\code{\link{alview}}, \code{\link{image.DNAbin}}
+\code{\link{alview}}, \code{\link{image.DNAbin}}, \code{\link{all.equal.DNAbin}}
}
\examples{
data(woodmouse)
diff --git a/man/clustal.Rd b/man/clustal.Rd
index a0d77df..f39a3ef 100644
--- a/man/clustal.Rd
+++ b/man/clustal.Rd
@@ -6,7 +6,10 @@
\title{Multiple Sequence Alignment with External Applications}
\description{
These functions call their respective program from \R to align a set of
- nucleotide sequences of class \code{"DNAbin"}.
+ nucleotide sequences of class \code{"DNAbin"}. The application(s) must
+ be installed seperately and it is highly recommended to do this so
+ that the executables are in a directory located on the PATH of the
+ system.
}
\usage{
clustal(x, pw.gapopen = 10, pw.gapext = 0.1,
@@ -33,15 +36,24 @@ tcoffee(x, exec = "t_coffee", MoreArgs = "", quiet = TRUE,
aligned sequences in the same order than in \code{x}.}
}
\details{
+ It is highly recommended to install the executables properly so that
+ they are in a directory located on the PATH (i.e., accessible from any
+ other directory). Alternatively, the full path to the executable
+ may be given (e.g., \code{exec = "~/muscle/muscle"}), or a (symbolic)
+ link may be copied in the working directory. For Debian and its
+ derivatives (e.g., Ubuntu), it is recommended to use the binaries
+ distributed by Debian.
+
\code{clustal} tries to guess the name of the executable program
depending on the operating system. Specifically, the followings are
used: ``clustalw'' under Linux, ``clustalw2'' under MacOS, or
``clustalw2.exe'' under Windows.
- For \code{clustalomega}, ``clustalo'' is the default on all systems
- (with no specific path) since it seems there is no Windows installer.
+ For \code{clustalomega}, ``clustalo[.exe]'' is the default on all
+ systems (with no specific path) since it seems there is no Windows
+ installer.
- The calculations are done in a temporary directory which is deleted
+ The computations are done in a temporary directory which is deleted
when \R is quit. So it is possible to find the files created by the
last call in the directory printed by \code{tempdir()}.
@@ -49,6 +61,9 @@ tcoffee(x, exec = "t_coffee", MoreArgs = "", quiet = TRUE,
function prints the options of the program which may be passed to
\code{MoreArgs}.
}
+\note{
+ \code{clustalomega} is actually for AA sequences only.
+}
\value{
an object of class \code{"DNAbin"} with the aligned sequences.
}
@@ -67,7 +82,7 @@ tcoffee(x, exec = "t_coffee", MoreArgs = "", quiet = TRUE,
Notredame, C., Higgins, D. and Heringa, J. (2000) T-Coffee: A novel
method for multiple sequence alignments. \emph{Journal of Molecular
Biology}, \bold{302}, 205--217.
- \url{http://www.tcoffee.org/Documentation/t_coffee/t_coffee_technical.htm}
+ \url{http://www.tcoffee.org/}
Sievers, F., Wilm, A., Dineen, D., Gibson, T. J., Karplus, K., Li, W.,
Lopez, R., McWilliam, H., Remmert, M., S\"oding, J., Thompson,
@@ -78,10 +93,8 @@ tcoffee(x, exec = "t_coffee", MoreArgs = "", quiet = TRUE,
\author{Emmanuel Paradis}
\seealso{
\code{\link{image.DNAbin}}, \code{\link{del.gaps}},
- \code{\link{alex}}, \code{\link{alview}}, \code{\link{checkAlignment}}
-
- The package \pkg{ips} which has similar functions for MAFFT and
- Prank.
+ \code{\link{all.equal.DNAbin}}, \code{\link{alex}},
+ \code{\link{alview}}, \code{\link{checkAlignment}}
}
\examples{
\dontrun{
diff --git a/man/dist.topo.Rd b/man/dist.topo.Rd
index 6c7814c..8c91553 100644
--- a/man/dist.topo.Rd
+++ b/man/dist.topo.Rd
@@ -20,14 +20,16 @@ dist.topo(x, y = NULL, method = "PH85")
a single numeric value.
}
\details{
- Two methods are available: the one by Penny and Hendy (1985), and the
- branch length score by Kuhner and Felsenstein (1994). The trees are
- always considered as unrooted.
+ Two methods are available: the one by Penny and Hendy (1985,
+ originally from Robinson and Foulds 1981), and the branch length score
+ by Kuhner and Felsenstein (1994). The trees are always considered as
+ unrooted.
The topological distance is defined as twice the number of internal
- branches defining different bipartitions of the tips (Penny and Hendy
- 1985). Rzhetsky and Nei (1992) proposed a modification of the original
- formula to take multifurcations into account.
+ branches defining different bipartitions of the tips (Robinson and
+ Foulds 1981; Penny and Hendy 1985). Rzhetsky and Nei (1992) proposed a
+ modification of the original formula to take multifurcations into
+ account.
The branch length score may be seen as similar to the previous
distance but taking branch lengths into account. Kuhner and
@@ -54,6 +56,9 @@ dist.topo(x, y = NULL, method = "PH85")
Penny, D. and Hendy, M. D. (1985) The use of tree comparison
metrics. \emph{Systemetic Zoology}, \bold{34}, 75--82.
+ Robinson, D. F. and Foulds, L. R. (1981) Comparison of phylogenetic
+ trees. \emph{Mathematical Biosciences}, \bold{53}, 131--147.
+
Rzhetsky, A. and Nei, M. (1992) A simple method for estimating and
testing minimum-evolution trees. \emph{Molecular Biology and
Evolution}, \bold{9}, 945--967.
@@ -66,7 +71,7 @@ dist.topo(x, y = NULL, method = "PH85")
\examples{
ta <- rtree(30)
tb <- rtree(30)
-dist.topo(ta, ta) # = 0
-dist.topo(ta, tb) # This is unlikely to be 0!
+dist.topo(ta, ta) # 0
+dist.topo(ta, tb) # unlikely to be 0
}
\keyword{manip}
diff --git a/man/image.DNAbin.Rd b/man/image.DNAbin.Rd
index fedb652..f5fdfec 100644
--- a/man/image.DNAbin.Rd
+++ b/man/image.DNAbin.Rd
@@ -44,8 +44,9 @@
\author{Emmanuel Paradis}
\seealso{
\code{\link{DNAbin}}, \code{\link{del.gaps}}, \code{\link{alex}},
- \code{\link{alview}}, \code{\link{clustal}},
- \code{\link[graphics]{grid}}, \code{\link{image.AAbin}}
+ \code{\link{alview}}, \code{\link{all.equal.DNAbin}},
+ \code{\link{clustal}}, \code{\link[graphics]{grid}},
+ \code{\link{image.AAbin}}
}
\examples{
data(woodmouse)
diff --git a/man/is.binary.tree.Rd b/man/is.binary.tree.Rd
index 1463da4..d868313 100644
--- a/man/is.binary.tree.Rd
+++ b/man/is.binary.tree.Rd
@@ -1,35 +1,43 @@
-\name{is.binary.tree}
+\name{is.binary}
+\alias{is.binary}
+\alias{is.binary.phylo}
+\alias{is.binary.multiPhylo}
\alias{is.binary.tree}
\title{Test for Binary Tree}
\description{
This function tests whether a phylogenetic tree is binary.
}
\usage{
-is.binary.tree(phy)
+is.binary(phy)
+\method{is.binary}{phylo}(phy)
+\method{is.binary}{multiPhylo}(phy)
+\method{is.binary}{tree}(phy)
}
\arguments{
- \item{phy}{an object of class \code{"phylo"}.}
+ \item{phy}{an object of class \code{"phylo"} or \code{"multiPhylo"}.}
}
\details{
- The test differs slightly whether the tree is rooted or not. An
- urooted tree is considered binary if all its nodes are of degree three
- (i.e., three edges connect to each node). A rooted tree is considered
- binary if all nodes (including the root node) have exactly two
- descendant nodes, so that they are of degree three expect the root
- which is of degree 2.
+ The test differs whether the tree is rooted or not. An urooted tree is
+ considered binary if all its nodes are of degree three (i.e., three
+ edges connect to each node). A rooted tree is considered binary if all
+ nodes (including the root node) have exactly two descendant nodes, so
+ that they are of degree three expect the root which is of degree 2.
+
+ \code{is.binary.tree} is deprecated and will be removed soon:
+ currently it calls \code{is.binary}.
}
\value{
- \code{is.binary.tree} returns \code{TRUE} if \code{tree} is a fully
- binary phylogenetic tree, otherwise it returns \code{FALSE}.
+ a logical vector.
}
\seealso{
\code{\link{is.rooted}}, \code{\link{is.ultrametric}}, \code{\link{multi2di}}
}
-\author{Korbinian Strimmer}
+\author{Emmanuel Paradis}
\examples{
-data("hivtree.newick") # example tree in NH format
-tree.hiv <- read.tree(text = hivtree.newick)
-is.binary.tree(tree.hiv)
-plot(tree.hiv)
+is.binary(rtree(10))
+is.binary(rtree(10, rooted = FALSE))
+is.binary(stree(10))
+x <- setNames(rmtree(10, 10), LETTERS[1:10])
+is.binary(x)
}
\keyword{logic}
diff --git a/man/is.ultrametric.Rd b/man/is.ultrametric.Rd
index b3fdab4..663c0ee 100644
--- a/man/is.ultrametric.Rd
+++ b/man/is.ultrametric.Rd
@@ -1,29 +1,40 @@
\name{is.ultrametric}
\alias{is.ultrametric}
+\alias{is.ultrametric.phylo}
+\alias{is.ultrametric.multiPhylo}
\title{Test if a Tree is Ultrametric}
+\description{
+ This function tests whether a tree is ultrametric using the distances
+ from each tip to the root.
+}
\usage{
-is.ultrametric(phy, tol = .Machine$double.eps^0.5)
+is.ultrametric(phy, ...)
+\method{is.ultrametric}{phylo}(phy, tol = .Machine$double.eps^0.5, option = 1, ...)
+\method{is.ultrametric}{multiPhylo}(phy, tol = .Machine$double.eps^0.5, option = 1, ...)
}
\arguments{
- \item{phy}{an object of class \code{"phylo"}.}
+ \item{phy}{an object of class \code{"phylo"} or \code{"multiPhylo"}.}
\item{tol}{a numeric >= 0, variation below this value are considered
- non-significant (see details).}
+ non-significant.}
+ \item{option}{an integer (1 or 2; see details).}
+ \item{\dots}{arguments passed among methods.}
}
-\description{
- This function computes the distances from each tip to the root: if the
- variance of these distances is null, the tree is considered as
- ultrametric.
+\details{
+ The test is based on the distances from each tip to the root and a
+ criterion: if \code{option = 1}, the criterion is the scaled range
+ ((max - min/max)), if \code{option = 2}, the variance is used (this
+ was the method used until ape 3.5). The default criterion is invariant
+ to linear changes of the branch lengths.
}
\value{
- a logical: \code{TRUE} if the tree is ultrametric, \code{FALSE}
- otherwise.
-}
-\details{
- The default value for \code{tol} is based on the numerical
- characteristics of the machine \R is running on.
+ a logical vector.
}
\author{Emmanuel Paradis}
\seealso{
- \code{\link{is.binary.tree}}, \code{\link[base]{.Machine}}
+ \code{\link{is.binary}}, \code{\link[base]{.Machine}}
+}
+\examples{
+is.ultrametric(rtree(10))
+is.ultrametric(rcoal(10))
}
\keyword{utilities}
diff --git a/man/landplants.Rd b/man/landplants.Rd
deleted file mode 100644
index 2b870e3..0000000
--- a/man/landplants.Rd
+++ /dev/null
@@ -1,36 +0,0 @@
-\name{landplants}
-\alias{landplants}
-\alias{landplants.newick}
-\title{Gene Tree of 36 Landplant rbcL Sequences}
-\description{
- This data set describes a gene tree estimated from 36 landplant
- \emph{rbc}L sequences.
-}
-\usage{
-data(landplants.newick)
-}
-\format{
- \code{landplants.newick} is a string with the tree in Newick format.
-}
-\source{
- This tree is described in Sanderson (1997) and is also a
- data example in the software package r8s.
-}
-\references{
- Sanderson, M. J. (1997) A nonparametric approach to estimating
- divergence times in the absence of rate constancy. \emph{Molecular
- Biology and Evolution}, \bold{14}, 1218--1231.
-}
-\examples{
-# example tree in NH format (a string)
-data("landplants.newick")
-landplants.newick
-
-# get corresponding phylo object
-tree.landplants <- read.tree(text = landplants.newick)
-
-# plot tree
-plot(tree.landplants, label.offset = 0.001)
-}
-\keyword{datasets}
-
diff --git a/man/mixedFontLabel.Rd b/man/mixedFontLabel.Rd
index 1455ff2..8b73672 100644
--- a/man/mixedFontLabel.Rd
+++ b/man/mixedFontLabel.Rd
@@ -29,7 +29,7 @@ mixedFontLabel(..., sep = " ", italic = NULL, bold = NULL,
The idea is to have different bits of text in different vectors that
are put together to make a vector of \R expressions. This vector is
interpreted by graphical functions to format the text. A simple use
- may be \code{mixedFontLabel(genus, species), italic = 1:2}, but it is
+ may be \code{mixedFontLabel(genus, species, italic = 1:2)}, but it is
more interesting when mixing fonts (see examples).
To have an element in bolditalics, its number must given in both
diff --git a/man/multi2di.Rd b/man/multi2di.Rd
index 2dd5f03..0fd0b97 100644
--- a/man/multi2di.Rd
+++ b/man/multi2di.Rd
@@ -1,22 +1,31 @@
\name{multi2di}
\alias{multi2di}
+\alias{multi2di.phylo}
+\alias{multi2di.multiPhylo}
\alias{di2multi}
+\alias{di2multi.phylo}
+\alias{di2multi.multiPhylo}
\title{Collapse and Resolve Multichotomies}
\description{
These two functions collapse or resolve multichotomies in phylogenetic
trees.
}
\usage{
-multi2di(phy, random = TRUE)
-di2multi(phy, tol = 1e-08)
+multi2di(phy, ...)
+\method{multi2di}{phylo}(phy, random = TRUE, ...)
+\method{multi2di}{multiPhylo}(phy, random = TRUE, ...)
+di2multi(phy, ...)
+\method{di2multi}{phylo}(phy, tol = 1e-08, ...)
+\method{di2multi}{multiPhylo}(phy, tol = 1e-08, ...)
}
\arguments{
- \item{phy}{an object of class \code{"phylo"}.}
+ \item{phy}{an object of class \code{"phylo"} or \code{"multiPhylo"}.}
\item{random}{a logical value specifying whether to resolve the
multichotomies randomly (the default) or in the order they appear in
the tree (if \code{random = FALSE}).}
\item{tol}{a numeric value giving the tolerance to consider a branch
length significantly greater than zero.}
+ \item{\dots}{arguments passed among methods.}
}
\details{
\code{multi2di} transforms all multichotomies into a series of
@@ -25,17 +34,15 @@ di2multi(phy, tol = 1e-08)
\code{di2multi} deletes all branches smaller than \code{tol} and
collapses the corresponding dichotomies into a multichotomy.
}
-\seealso{
-\code{\link{is.binary.tree}}
-}
+\seealso{\code{\link{is.binary}}}
\author{Emmanuel Paradis}
\value{
- Both functions return an object of class \code{"phylo"}.
+ an object of the same class than the input.
}
\examples{
data(bird.families)
-is.binary.tree(bird.families)
-is.binary.tree(multi2di(bird.families))
+is.binary(bird.families)
+is.binary(multi2di(bird.families))
all.equal(di2multi(multi2di(bird.families)), bird.families)
### To see the results of randomly resolving a trichotomy:
tr <- read.tree(text = "(a:1,b:1,c:1);")
diff --git a/man/node.dating.Rd b/man/node.dating.Rd
new file mode 100644
index 0000000..0667b26
--- /dev/null
+++ b/man/node.dating.Rd
@@ -0,0 +1,107 @@
+\name{node.dating}
+\alias{estimate.mu}
+\alias{estimate.dates}
+\title{node.dating}
+\description{
+ Estimate the dates of a rooted phylogenetic tree from the tip dates.
+}
+\usage{
+estimate.mu(t, node.dates, p = 0.05)
+estimate.dates(t, node.dates, mu = estimate.mu(t, node.dates),
+ min.date = -.Machine$double.xmax, show.steps = 0,
+ opt.tol = 1e-8, nsteps = 1000,
+ lik.tol = if (nsteps <= 0) opt.tol else 0,
+ is.binary = is.binary.phylo(t))
+}
+\arguments{
+ \item{t}{an object of class '"phylo"'}
+ \item{node.dates}{a numeric vector of dates for the tips, in the same order as
+ 't$tip.label' or a vector of dates for all of the nodes.}
+ \item{p}{p-value cutoff for failed regression.}
+ \item{mu}{mutation rate.}
+ \item{min.date}{the minimum bound on the dates of nodes.}
+ \item{show.steps}{print the log likelihood every show.steps. If 0 will supress
+ output.}
+ \item{opt.tol}{tolerance for optimization precision.}
+ \item{lik.tol}{tolerance for likelihood comparison.}
+ \item{nsteps}{the maximum number of steps to run.}
+ \item{is.binary}{if TRUE, will run a faster optimization method that only works if
+ the tree is binary; otherwise will use optimize() as the optimization method.}
+}
+\value{
+ The estimated mutation rate as a numeric vector of length one for estimate.mu.
+
+ The estimated dates of all of the nodes of the tree as a numeric vector with
+ length equal to the number of nodes in the tree.
+}
+\details{
+ This code duplicates the functionality of the program Tip.Dates (see references).
+ The dates of the internal nodes of 't' are estimated using a maximum likelihood
+ approach.
+
+ 't' must be rooted and have branch lengths in units of expected substitutions per
+ site.
+
+ 'node.dates' can be either a numeric vector of dates for the tips or a numeric
+ vector for all of the nodes of 't'. 'estimate.mu' will use all of the values
+ given in 'node.dates' to estimate the mutation rate. Dates can be censored with
+ NA. 'node.dates' must contain all of the tip dates when it is a parameter of
+ 'estimate.dates'. If only tip dates are given, then 'estimate.dates' will run an
+ initial step to estimate the dates of the internal nodes. If 'node.dates'
+ contains dates for some of the nodes, 'estimate.dates' will use those dates as
+ priors in the inital step. If all of the dates for nodes are given, then
+ 'estimate.dates' will not run the inital step.
+
+ If 'is.binary' is set to FALSE, 'estimate.dates' uses the '"optimize"' function as
+ the optimization method. By default, R's '"optimize"' function uses a precision
+ of '".Machine$double.eps^0.25"', which is about 0.0001 on a 64-bit system. This
+ should be set to a smaller value if the branch lengths of 't' are very short. If
+ 'is.binary' is set to TRUE, estimate dates uses calculus to deterimine the maximum
+ likelihood at each step, which is faster. The bounds of permissible values are
+ reduced by 'opt.tol'.
+
+ 'estimate.dates' has several criteria to decide how many steps it will run. If
+ 'lik.tol' and 'nsteps' are both 0, then 'estimate.dates' will only run the initial
+ step. If 'lik.tol' is greater than 0 and 'nsteps' is 0, then 'estimate.dates'
+ will run until the difference between successive steps is less than 'lik.tol'. If
+ 'lik.tol' is 0 and 'nsteps' is greater than 0, then 'estimate.dates' will run the
+ inital step and then 'nsteps' steps. If 'lik.tol' and 'nsteps' are both greater
+ than 0, then 'estimate.dates' will run the inital step and then either 'nsteps'
+ steps or until the difference between successive steps is less than 'lik.tol'.
+}
+\note{
+ This model assumes that the tree follows a molecular clock. It only performs a
+ rudimentary statistical test of the molecular clock hypothesis.
+}
+\author{Bradley R. Jones <email: brj1 at sfu.ca>}
+\references{
+ Felsenstein, J. (1981) Evolutionary trees from DNA sequences: a maximum likelihood
+ approach. \emph{Journal of Molecular Evolution}, \bold{17}, 368--376.
+
+ Rambaut, A. (2000) Estimating the rate of molecular evolution:
+ incorporating non-contemporaneous sequences into maximum likelihood
+ phylogenies. \emph{Bioinformatics}, \bold{16}, 395--399.
+
+ Jones, Bradley R., and Poon, Art F. Y. (2016) \emph{node.dating:
+ dating ancestors in phylogenetic trees in R}. Preprint.
+}
+\seealso{
+ \code{\link[stats]{optimize}}, \code{\link{rtt}},
+ \code{\link{plotTreeTime}}
+}
+\examples{
+t <- rtree(100)
+tip.date <- rnorm(t$tip.label, mean = node.depth.edgelength(t)[1:Ntip(t)])^2
+t <- rtt(t, tip.date)
+mu <- estimate.mu(t, tip.date)
+
+## Run for 100 steps
+node.date <- estimate.dates(t, tip.date, mu, nsteps = 100)
+
+## Run until the difference between successive log likelihoods is
+## less than $10^{-4}$ starting with the 100th step's results
+node.date <- estimate.dates(t, node.date, mu, nsteps = 0, lik.tol = 1e-4)
+
+## To rescale the tree over time
+t$edge.length <- node.date[t$edge[, 2]] - node.date[t$edge[, 1]]
+}
diff --git a/man/opsin.Rd b/man/opsin.Rd
deleted file mode 100644
index 8b56c7e..0000000
--- a/man/opsin.Rd
+++ /dev/null
@@ -1,35 +0,0 @@
-\name{opsin}
-\alias{opsin}
-\alias{opsin.newick}
-\title{Gene Tree of 32 opsin Sequences}
-\description{
- This data set describes a gene tree estimated from 32 opsin
- sequences.
-}
-\usage{
-data(opsin.newick)
-}
-\format{
- \code{opsin.newick} is a string with the tree in Newick format.
-}
-\source{
- This tree is described in Misawa and Tajima (2000) as an example
- for application of the Klastorin (1982) classification method.
-}
-\references{
- Misawa, K. (2000) A simple method for classifying genes and a bootstrap
- test for classifications. \emph{Molecular
- Biology and Evolution}, \bold{17}, 1879--1884.
-}
-\examples{
-# example tree in NH format (a string)
-data("opsin.newick")
-opsin.newick
-
-# get corresponding phylo object
-tree.opsin <- read.tree(text = opsin.newick)
-
-# plot tree
-plot(tree.opsin, label.offset = 0.01)
-}
-\keyword{datasets}
diff --git a/man/plot.phyloExtra.Rd b/man/plot.phyloExtra.Rd
index b839972..78bcfa3 100644
--- a/man/plot.phyloExtra.Rd
+++ b/man/plot.phyloExtra.Rd
@@ -28,7 +28,7 @@ drawSupportOnEdges(value, ...)
\author{Emmanuel Paradis}
\seealso{
\code{\link{plot.phylo}}, \code{\link{edgelabels}},
- \code{\link{boot.phylo}}
+ \code{\link{boot.phylo}}, \code{\link{plotTreeTime}}
}
\examples{
tr <- rtree(10)
diff --git a/man/plotTreeTime.Rd b/man/plotTreeTime.Rd
new file mode 100644
index 0000000..494c0bc
--- /dev/null
+++ b/man/plotTreeTime.Rd
@@ -0,0 +1,44 @@
+\name{plotTreeTime}
+\alias{plotTreeTime}
+\title{Plot Tree With Time Axis}
+\description{
+ This function plots a non-ultrametric tree where the tips are not
+ contemporary together with their dates on the x-axis.
+}
+\usage{
+plotTreeTime(phy, tip.dates, show.tip.label = FALSE, y.lim = NULL, ...)
+}
+\arguments{
+ \item{phy}{an object of class \code{"phylo"}.}
+ \item{tip.dates}{a vector of the same length than the number of tips
+ in \code{phy} (see details).}
+ \item{show.tip.label}{a logical value; see \code{\link{plot.phylo}}.}
+ \item{y.lim}{by default, one fifth of the plot is left below the tree;
+ use this option to change this behaviour.}
+ \item{\dots}{other arguments to be passed to \code{plot.phylo}.}
+}
+\details{
+ The vector \code{tip.dates} may be numeric or of class
+ \dQuote{\link[base]{Date}}. In either case, the time axis is set
+ accordingly. The length of this vector must be equal to the number of
+ tips of the tree: the dates are matched to the tips numbers. Missing
+ values are allowed.
+}
+\value{NULL}
+\author{Emmanuel Paradis}
+\seealso{
+ \code{\link{plot.phylo}}, \code{\link{estimate.dates}}
+}
+\examples{
+dates <- as.Date(.leap.seconds)
+tr <- rtree(length(dates))
+plotTreeTime(tr, dates)
+
+## handling NA's:
+dates[11:26] <- NA
+plotTreeTime(tr, dates)
+
+## dates can be on an arbitrary scale, e.g., [-1, 1]:
+plotTreeTime(tr, runif(Ntip(tr), -1, 1))
+}
+\keyword{hplot}
diff --git a/man/read.GenBank.Rd b/man/read.GenBank.Rd
index 42ca2e5..6790fdd 100644
--- a/man/read.GenBank.Rd
+++ b/man/read.GenBank.Rd
@@ -5,39 +5,41 @@
read.GenBank(access.nb, seq.names = access.nb, species.names = TRUE,
gene.names = FALSE, as.character = FALSE)
}
+\description{
+ This function connects to the GenBank database, and reads nucleotide
+ sequences using accession numbers given as arguments.
+}
\arguments{
\item{access.nb}{a vector of mode character giving the accession numbers.}
\item{seq.names}{the names to give to each sequence; by default the
accession numbers are used.}
\item{species.names}{a logical indicating whether to attribute the
species names to the returned object.}
- \item{gene.names}{a logical indicating whether to attribute the
- gene names to the returned object. It is \code{FALSE} by default
- because this will work correctly only when reading sequences with a
- single gene.}
+ \item{gene.names}{obsolete (will be removed soon).}
\item{as.character}{a logical controlling whether to return the
sequences as an object of class \code{"DNAbin"} (the default).}
}
-\description{
- This function connects to the GenBank database, and reads nucleotide
- sequences using accession numbers given as arguments.
-}
-\value{
- A list of DNA sequences made of vectors of class \code{"DNAbin"}, or
- of single characters (if \code{as.character = "TRUE"}).
-}
\details{
The function uses the site \url{http://www.ncbi.nlm.nih.gov/} from
- where the sequences are downloaded.
+ where the sequences are retrieved.
If \code{species.names = TRUE}, the returned list has an attribute
\code{"species"} containing the names of the species taken from the
field ``ORGANISM'' in GenBank.
- If \code{gene.names = TRUE}, the returned list has an attribute
- \code{"gene"} containing the names of the gene. This will not work
- correctly if reading a sequence with multiple genes (e.g., a
- mitochondrial genome).
+ Since ape 3.6, this function retrieves the sequences in FASTA format:
+ this is more efficient and more flexible (scaffolds and contigs can be
+ read). The option \code{gene.names} is obsolete and will be removed;
+ this information is also present in the description.
+
+ Setting \code{species.names = FALSE} is quite faster (could be useful
+ if you read a series of scaffolds or contigs, or if you already have
+ the species names).
+}
+\value{
+ A list of DNA sequences made of vectors of class \code{"DNAbin"}, or
+ of single characters (if \code{as.character = TRUE}) with two
+ attributes (species and description).
}
\seealso{
\code{\link{read.dna}}, \code{\link{write.dna}},
@@ -45,26 +47,26 @@ read.GenBank(access.nb, seq.names = access.nb, species.names = TRUE,
}
\author{Emmanuel Paradis}
\examples{
-### This won't work if your computer is not connected
-### to the Internet
-###
-### Get the 8 sequences of tanagers (Ramphocelus)
-### as used in Paradis (1997)
+## This won't work if your computer is not connected
+## to the Internet
+
+## Get the 8 sequences of tanagers (Ramphocelus)
+## as used in Paradis (1997)
ref <- c("U15717", "U15718", "U15719", "U15720",
"U15721", "U15722", "U15723", "U15724")
-### Copy/paste or type the following commands if you
-### want to try them.
+## Copy/paste or type the following commands if you
+## want to try them.
\dontrun{
Rampho <- read.GenBank(ref)
-### get the species names:
+## get the species names:
attr(Rampho, "species")
-### build a matrix with the species names and the accession numbers:
+## build a matrix with the species names and the accession numbers:
cbind(attr(Rampho, "species"), names(Rampho))
-### print the first sequence
-### (can be done with `Rampho$U15717' as well)
+## print the first sequence
+## (can be done with `Rampho$U15717' as well)
Rampho[[1]]
-### print the first sequence in a cleaner way
-cat(Rampho[[1]], "\n", sep = "")
+## the description from each FASTA sequence:
+attr(Rampho, "description")
}
}
\keyword{IO}
diff --git a/man/read.dna.Rd b/man/read.dna.Rd
index caba0d2..00fee97 100644
--- a/man/read.dna.Rd
+++ b/man/read.dna.Rd
@@ -6,7 +6,7 @@
These functions read DNA sequences in a file, and returns a matrix or a
list of DNA sequences with the names of the taxa read in the file as
rownames or names, respectively. By default, the sequences are stored
- in binary format, otherwise (if \code{as.character = "TRUE"}) in lower
+ in binary format, otherwise (if \code{as.character = TRUE}) in lower
case.
}
\usage{
@@ -133,7 +133,7 @@ cat(">No305",
"ATTCGAAAAACACACCCACTACTAAAAATTATCAATCACT",
file = "exdna.txt", sep = "\n")
ex.dna4 <- read.dna("exdna.txt", format = "fasta")
-### The first three are the same!
+### They are the same!
identical(ex.dna, ex.dna2)
identical(ex.dna, ex.dna3)
identical(ex.dna, ex.dna4)
diff --git a/man/read.gff.Rd b/man/read.gff.Rd
new file mode 100644
index 0000000..11bd395
--- /dev/null
+++ b/man/read.gff.Rd
@@ -0,0 +1,46 @@
+\name{read.gff}
+\alias{read.gff}
+\title{Read GFF Files}
+\description{
+ This function reads a file in general feature format version 3 (GFF3)
+ and returns a data frame.
+}
+\usage{
+read.gff(file, na.strings = c(".", "?"))
+}
+\arguments{
+ \item{file}{a file name specified by a character string.}
+ \item{na.strings}{the strings in the GFF file that will be converted
+ as NA's (missing values).}
+}
+\details{
+ The returned data frame has its (column) names correctly set (see
+ References) and the categorical variables (seqid, source, type,
+ strand, and phase) set as factors.
+
+ This function should be more efficient than using \code{read.delim}.
+
+ GFF2 files can also be read but the field names conform to GFF3.
+
+ The file can be gz-compressed (see examples), but not zipped.
+}
+\value{NULL}
+\author{Emmanuel Paradis}
+\references{
+ \url{https://en.wikipedia.org/wiki/General_feature_format}
+}
+\examples{
+\dontrun{
+## requires to be connected on Internet
+d <- "ftp://ftp.ensembl.org/pub/release-86/gff3/homo_sapiens/"
+f <- "Homo_sapiens.GRCh38.86.chromosome.MT.gff3.gz"
+download.file(paste0(d, f), "mt_gff3.gz")
+gff.mito <- read.gff("mt_gff3.gz")
+## the lengths of the sequence features:
+gff.mito$end - (gff.mito$start - 1)
+table(gff.mito$type)
+## where the exons start:
+gff.mito$start[gff.mito$type == "exon"]
+}
+}
+\keyword{IO}
diff --git a/man/reconstruct.Rd b/man/reconstruct.Rd
index fd56490..8168d63 100644
--- a/man/reconstruct.Rd
+++ b/man/reconstruct.Rd
@@ -3,33 +3,56 @@
\title{Continuous Ancestral Character Estimation}
\description{
This function estimates ancestral character states, and the associated
- uncertainty, for continuous characters. It mainly works as the
- \code{\link{ace}} function, from which it differs in the fact that
- optimisations are not performed by numerical algorithms but through
- matrix computations.
+ uncertainty, for continuous characters. It mainly works as the ace
+ function, from which it differs, first, in the fact that computations
+ are not performed by numerical optimisation but through matrix
+ calculus. Second, besides classical Brownian-based reconstruction
+ methods, it reconstructs ancestral states under Arithmetic Brownian
+ Motion (ABM, i.e. Brownian with linear trend) and Ornstein-Uhlenbeck
+ process (OU, i.e. Brownian with an attractive optimum).
}
\usage{
-reconstruct(x, phyInit, method = "ML", CI = TRUE)
+reconstruct(x, phyInit, method = "ML", alpha = NULL, CI = TRUE)
}
\arguments{
- \item{x}{a vector or a factor.}
+ \item{x}{a numerical vector.}
\item{phyInit}{an object of class \code{"phylo"}.}
\item{method}{a character specifying the method used for
- estimation. Three choices are possible: \code{"ML"}, \code{"REML"}
- or \code{"GLS"}.}
+ estimation. Six choices are possible: \code{"ML"}, \code{"REML"},
+ \code{"GLS"}, \code{"GLS_ABM"}, \code{"GLS_OU"} or
+ \code{"GLS_OUS"}.}
+ \item{alpha}{a numerical value which accounts for the attractive
+ strength parameter of \code{"GLS_OU"} or \code{"GLS_OUS"} (used only
+ in these cases). If alpha = NULL (the default), then it is estimated
+ by maximum likelihood using \code{optim} which may lead to
+ convergence issue.}
\item{CI}{a logical specifying whether to return the 95\% confidence
intervals of the ancestral state estimates.}
}
\details{
- The default model is Brownian motion where characters evolve randomly
- following a random walk. This model can be fitted by maximum
- likelihood (Felsenstein 1973, Schluter et al. 1997 - the default),
- residual maximum likelihood, or generalized least squares
- (\code{method = "GLS"}, Martins and Hansen 1997, Cunningham et
- al. 1998).
+ For \code{"ML"}, \code{"REML"} and \code{"GLS"}, the default model is
+ Brownian motion. This model can be fitted by maximum likelihood
+ (\code{method = "ML"}, Felsenstein 1973, Schluter et al. 1997) - the
+ default, residual maximum likelihood (\code{method = "REML"}), or
+ generalized least squares (\code{method = "GLS"}, Martins and Hansen
+ 1997, Garland T and Ives AR 2000).
+
+ \code{"GLS_ABM"} is based on Brownian motion with trend model. Both
+ \code{"GLS_OU"} and \code{"GLS_OUS"} are based on Ornstein-Uhlenbeck
+ model.
+
+ \code{"GLS_OU"} and \code{"GLS_OUS"} differs in the fact that
+ \code{"GLS_OUS"} assume that the process starts from the optimum,
+ while the root state has to be estimated for \code{"GLS_OU"}, which
+ may rise some issues (see Royer-Carenzi and Didier, 2016). Users may
+ provide the attractive strength parameter \code{alpha}, for these two
+ models.
+
+ \code{"GLS_ABM"}, \code{"GLS_OU"} and \code{"GLS_OUS"} are all fitted
+ by generalized least squares (Royer-Carenzi and Didier, 2016).
}
\value{
- a list with the following elements:
+ an object of class \code{"ace"} with the following elements:
\item{ace}{the estimates of the ancestral character values.}
\item{CI95}{the estimated 95\% confidence intervals.}
@@ -38,19 +61,24 @@ reconstruct(x, phyInit, method = "ML", CI = TRUE)
\item{loglik}{if \code{method = "ML"}, the maximum log-likelihood.}
}
\references{
- Cunningham, C. W., Omland, K. E. and Oakley, T. H. (1998)
- Reconstructing ancestral character states: a critical
- reappraisal. \emph{Trends in Ecology & Evolution}, \bold{13},
- 361--366.
-
Felsenstein, J. (1973) Maximum likelihood estimation of evolutionary
trees from continuous characters. \emph{American Journal of Human
Genetics}, \bold{25}, 471--492.
+ Garland T. and Ives A.R. (2000) Using the past to predict the present:
+ confidence intervals for regression equations in phylogenetic
+ comparative methods. \emph{American Naturalist}, \bold{155},
+ 346--364.
+
Martins, E. P. and Hansen, T. F. (1997) Phylogenies and the
comparative method: a general approach to incorporating phylogenetic
information into the analysis of interspecific data. \emph{American
- Naturalist}, \bold{149}, 646--667.
+ Naturalist}, \bold{149}, 646--667.
+
+ Royer-Carenzi, M. and Didier, G. (2016) A comparison of ancestral
+ state reconstruction methods for quantitative
+ characters. \emph{Journal of Theoretical Biology}, \bold{404},
+ 126--142.
Schluter, D., Price, T., Mooers, A. O. and Ludwig, D. (1997)
Likelihood of ancestor states in adaptive radiation. \emph{Evolution},
@@ -58,10 +86,6 @@ reconstruct(x, phyInit, method = "ML", CI = TRUE)
Yang, Z. (2006) \emph{Computational Molecular Evolution}. Oxford:
Oxford University Press.
-
- Royer-Carenzi, M. and Didier, G. (2014) Comparison of ancestral state
- reconstruction methods for continuous characters under directional
- evolution. Submitted.
}
\author{Manuela Royer-Carenzi, Gilles Didier}
\seealso{
@@ -70,12 +94,17 @@ reconstruct(x, phyInit, method = "ML", CI = TRUE)
Reconstruction of ancestral sequences can be done with the package
\pkg{phangorn} (see function \code{?ancestral.pml}).
}
+\note{
+\code{GLS_ABM} should not be used on ultrametric tree.
+
+\code{GLS_OU} may lead to aberrant reconstructions.
+}
\examples{
### Some random data...
data(bird.orders)
-x <- rnorm(23)
-### Compare the three methods for continuous characters:
+x <- rnorm(23, m=100)
+### Reconstruct ancestral quantitative characters:
reconstruct(x, bird.orders)
-reconstruct(x, bird.orders, method = "REML")
+reconstruct(x, bird.orders, method = "GLS_OUS", alpha=NULL)
}
\keyword{models}
diff --git a/man/reorder.phylo.Rd b/man/reorder.phylo.Rd
index c81beee..d34654f 100644
--- a/man/reorder.phylo.Rd
+++ b/man/reorder.phylo.Rd
@@ -1,5 +1,6 @@
\name{reorder.phylo}
\alias{reorder.phylo}
+\alias{reorder.multiPhylo}
\title{Internal Reordering of Trees}
\description{
This function changes the internal structure of a phylogeny stored as
@@ -8,9 +9,10 @@
}
\usage{
\method{reorder}{phylo}(x, order = "cladewise", index.only = FALSE, ...)
+\method{reorder}{multiPhylo}(x, order = "cladewise", ...)
}
\arguments{
- \item{x}{an object of class \code{"phylo"}.}
+ \item{x}{an object of class \code{"phylo"} or \code{"multiPhylo"}.}
\item{order}{a character string: either \code{"cladewise"} (the
default), \code{"postorder"}, \code{"pruningwise"}, or any
unambiguous abbreviation of these.}
@@ -39,7 +41,9 @@
}
\value{
an object of class \code{"phylo"} (with the attribute \code{"order"}
- set accordingly), or a numeric vector if \code{index.only = TRUE}.
+ set accordingly), or a numeric vector if \code{index.only = TRUE}; if
+ \code{x} is of class \code{"multiPhylo"}, then an object of the same
+ class.
}
\references{
Valiente, G. (2002) \emph{Algorithms on Trees and Graphs.} New York:
diff --git a/man/root.Rd b/man/root.Rd
index b20fe46..046cea3 100644
--- a/man/root.Rd
+++ b/man/root.Rd
@@ -1,32 +1,52 @@
\name{root}
\alias{root}
+\alias{root.phylo}
+\alias{root.multiPhylo}
\alias{unroot}
+\alias{unroot.phylo}
+\alias{unroot.multiPhylo}
\alias{is.rooted}
+\alias{is.rooted.phylo}
+\alias{is.rooted.multiPhylo}
\title{Roots Phylogenetic Trees}
+\description{
+ \code{root} reroots a phylogenetic tree with respect to the specified
+ outgroup or at the node specified in \code{node}.
+
+ \code{unroot} unroots a phylogenetic tree, or returns it unchanged if
+ it is already unrooted.
+
+ \code{is.rooted} tests whether a tree is rooted.
+}
\usage{
-root(phy, outgroup, node = NULL, resolve.root = FALSE, interactive = FALSE)
+root(phy, ...)
+\method{root}{phylo}(phy, outgroup, node = NULL, resolve.root = FALSE,
+ interactive = FALSE, edgelabel = FALSE, ...)
+\method{root}{multiPhylo}(phy, outgroup, ...)
+
unroot(phy)
+\method{unroot}{phylo}(phy)
+\method{unroot}{multiPhylo}(phy)
+
is.rooted(phy)
+\method{is.rooted}{phylo}(phy)
+\method{is.rooted}{multiPhylo}(phy)
}
\arguments{
- \item{phy}{an object of class \code{"phylo"}.}
+ \item{phy}{an object of class \code{"phylo"} or \code{"multiPhylo"}.}
\item{outgroup}{a vector of mode numeric or character specifying the
new outgroup.}
- \item{node}{alternatively, a node number where to root the tree (this
- should give the MRCA of the ingroup).}
+ \item{node}{alternatively, a node number where to root the tree.}
\item{resolve.root}{a logical specifying whether to resolve the new
root as a bifurcating node.}
\item{interactive}{if \code{TRUE} the user is asked to select the node
by clicking on the tree which must be plotted.}
-}
-\description{
- \code{root} reroots a phylogenetic tree with respect to the specified
- outgroup or at the node specified in \code{node}.
-
- \code{unroot} unroots a phylogenetic tree, or returns it unchanged if
- it is already unrooted.
-
- \code{is.rooted} tests whether a tree is rooted.
+ \item{edgelabel}{a logical value specifying whether to treat node
+ labels as edge labels and thus eventually switching them so that
+ they are associated with the correct edges when using
+ \code{\link{drawSupportOnEdges}} (see Czech et al. 2016).}
+ \item{\dots}{arguments passed among methods (e.g., when rooting lists
+ of trees).}
}
\details{
The argument \code{outgroup} can be either character or numeric. In
@@ -55,7 +75,7 @@ is.rooted(phy)
The use of \code{resolve.root = TRUE} together with \code{node = }
gives an error if the specified node is the current root of the
tree. This is because there is an ambiguity when resolving a node in
- an unrooted tree if no explicit outgroup. If the node is not the
+ an unrooted tree with no explicit outgroup. If the node is not the
current root, the ambiguity is solved arbitrarily by considering the
clade on the right of \code{node} (when the tree is plotted by
default) as the ingroup. See a detailed explanation there:
@@ -63,8 +83,13 @@ is.rooted(phy)
\url{http://www.mail-archive.com/r-sig-phylo@r-project.org/msg03805.html}.
}
\value{
- an object of class \code{"phylo"} for \code{root} and \code{unroot}; a
- single logical value for \code{is.rooted}.
+ an object of class \code{"phylo"} or \code{"multiPhylo"} for
+ \code{root} and \code{unroot}; a logical vector for \code{is.rooted}.
+}
+\references{
+ Czech, L., Huerta-Cepas, J. and Stamatakis, A. (2016) A critical
+ review on the use of support values in tree viewers and bioinformatics
+ toolkits. \url{http://dx.doi.org/10.1101/035360}
}
\author{Emmanuel Paradis}
\seealso{
@@ -77,8 +102,8 @@ plot(root(bird.orders, 1))
plot(root(bird.orders, 1:5))
tr <- root(bird.orders, 1)
-is.rooted(bird.orders) # yes!
-is.rooted(tr) # no!
+is.rooted(bird.orders) # yes
+is.rooted(tr) # no
### This is because the tree has been unrooted first before rerooting.
### You can delete the outgroup...
is.rooted(drop.tip(tr, "Struthioniformes"))
@@ -88,5 +113,8 @@ is.rooted(root(bird.orders, 1, r = TRUE))
### To keep the basal trichotomy but forcing the tree as rooted:
tr$root.edge <- 0
is.rooted(tr)
+
+x <- setNames(rmtree(10, 10), LETTERS[1:10])
+is.rooted(x)
}
\keyword{manip}
diff --git a/man/summary.phylo.Rd b/man/summary.phylo.Rd
index 984b65a..b00608d 100644
--- a/man/summary.phylo.Rd
+++ b/man/summary.phylo.Rd
@@ -1,17 +1,33 @@
\name{summary.phylo}
\alias{summary.phylo}
\alias{Ntip}
+\alias{Ntip.phylo}
+\alias{Ntip.multiPhylo}
\alias{Nnode}
+\alias{Nnode.phylo}
+\alias{Nnode.multiPhylo}
\alias{Nedge}
+\alias{Nedge.phylo}
+\alias{Nedge.multiPhylo}
\title{Print Summary of a Phylogeny}
\usage{
\method{summary}{phylo}(object, \dots)
+
Ntip(phy)
-Nnode(phy, internal.only = TRUE)
+\method{Ntip}{phylo}(phy)
+\method{Ntip}{multiPhylo}(phy)
+
+Nnode(phy, ...)
+\method{Nnode}{phylo}(phy, internal.only = TRUE, ...)
+\method{Nnode}{multiPhylo}(phy, internal.only = TRUE, ...)
+
Nedge(phy)
+\method{Nedge}{phylo}(phy)
+\method{Nedge}{multiPhylo}(phy)
}
\arguments{
- \item{object, phy}{an object of class \code{"phylo"}.}
+ \item{object, phy}{an object of class \code{"phylo"} or
+ \code{"multiPhylo"}.}
\item{\dots}{further arguments passed to or from other methods.}
\item{internal.only}{a logical indicating whether to return the number
of internal nodes only (the default), or of internal and terminal
diff --git a/src/ape.c b/src/ape.c
index 9e6c7a3..dc0a726 100644
--- a/src/ape.c
+++ b/src/ape.c
@@ -1,4 +1,4 @@
-/* ape.c 2016-02-18 */
+/* ape.c 2016-06-21 */
/* Copyright 2011-2016 Emmanuel Paradis, and 2007 R Development Core Team */
@@ -90,7 +90,6 @@ SEXP bitsplits_multiPhylo(SEXP x, SEXP n, SEXP nr);
static R_CMethodDef C_entries[] = {
{"C_additive", (DL_FUNC) &C_additive, 4},
- {"BaseProportion", (DL_FUNC) &BaseProportion, 3},
{"C_bionj", (DL_FUNC) &C_bionj, 5},
{"C_bionjs", (DL_FUNC) &C_bionjs, 6},
{"delta_plot", (DL_FUNC) &delta_plot, 5},
@@ -113,12 +112,10 @@ static R_CMethodDef C_entries[] = {
{"node_height_clado", (DL_FUNC) &node_height_clado, 7},
{"C_pic", (DL_FUNC) &C_pic, 10},
{"C_rTraitCont", (DL_FUNC) &C_rTraitCont, 9},
- {"SegSites", (DL_FUNC) &SegSites, 4},
{"C_treePop", (DL_FUNC) &C_treePop, 7},
{"C_triangMtd", (DL_FUNC) &C_triangMtd, 5},
{"C_triangMtds", (DL_FUNC) &C_triangMtds, 5},
{"C_ultrametric", (DL_FUNC) &C_ultrametric, 4},
- {"C_where", (DL_FUNC) &C_where, 6},
{"bitsplits_phylo", (DL_FUNC) &bitsplits_phylo, 6},
{"CountBipartitionsFromTrees", (DL_FUNC) &CountBipartitionsFromTrees, 8},
{"DNAbin2indelblock", (DL_FUNC) &DNAbin2indelblock, 4},
@@ -133,6 +130,9 @@ static R_CallMethodDef Call_entries[] = {
{"seq_root2tip", (DL_FUNC) &seq_root2tip, 3},
{"treeBuildWithTokens", (DL_FUNC) &treeBuildWithTokens, 1},
{"bitsplits_multiPhylo", (DL_FUNC) &bitsplits_multiPhylo, 3},
+ {"BaseProportion", (DL_FUNC) &BaseProportion, 1},
+ {"SegSites", (DL_FUNC) &SegSites, 1},
+ {"C_where", (DL_FUNC) &C_where, 2},
{NULL, NULL, 0}
};
diff --git a/src/dist_dna.c b/src/dist_dna.c
index 41a9bbb..a90ff2c 100644
--- a/src/dist_dna.c
+++ b/src/dist_dna.c
@@ -1,4 +1,4 @@
-/* dist_dna.c 2016-02-18 */
+/* dist_dna.c 2016-11-28 */
/* Copyright 2005-2016 Emmanuel Paradis */
@@ -1039,45 +1039,65 @@ void distDNA_ParaLin_pairdel(unsigned char *x, int *n, int *s, double *d,
}
/* a hash table is much faster than switch (2012-01-10) */
-void BaseProportion(unsigned char *x, int *n, double *BF)
+SEXP BaseProportion(SEXP x)
{
- int i;
- double count[256];
-
- memset(count, 0, 256*sizeof(double));
-
- for (i = 0; i < *n; i++) count[x[i]]++;
-
- BF[0] = count[136];
- BF[1] = count[40];
- BF[2] = count[72];
- BF[3] = count[24];
- BF[4] = count[192];
- BF[5] = count[160];
- BF[6] = count[144];
- BF[7] = count[96];
- BF[8] = count[80];
- BF[9] = count[48];
- BF[10] = count[224];
- BF[11] = count[176];
- BF[12] = count[208];
- BF[13] = count[112];
- BF[14] = count[240];
- BF[15] = count[4];
- BF[16] = count[2];
+ long i;
+ unsigned char *p;
+ double n, count[256], *BF;
+ SEXP res;
+
+ PROTECT(x = coerceVector(x, RAWSXP));
+ memset(count, 0, 256*sizeof(double));
+ n = XLENGTH(x);
+ p = RAW(x);
+
+ for (i = 0; i < n; i++) count[p[i]]++;
+
+ PROTECT(res = allocVector(REALSXP, 17));
+ BF = REAL(res);
+ BF[0] = count[136];
+ BF[1] = count[40];
+ BF[2] = count[72];
+ BF[3] = count[24];
+ BF[4] = count[192];
+ BF[5] = count[160];
+ BF[6] = count[144];
+ BF[7] = count[96];
+ BF[8] = count[80];
+ BF[9] = count[48];
+ BF[10] = count[224];
+ BF[11] = count[176];
+ BF[12] = count[208];
+ BF[13] = count[112];
+ BF[14] = count[240];
+ BF[15] = count[4];
+ BF[16] = count[2];
+ UNPROTECT(2);
+ return res;
}
#define SEGCOL seg[j] = 1; done = 1; break
-void SegSites(unsigned char *x, int *n, int *s, int *seg)
+SEXP SegSites(SEXP DNASEQ)
{
- int i, j, done, end;
- unsigned char base;
+ int n, s, j, done, *seg;
+ long i, end;
+ unsigned char base, *x;
+ SEXP ans;
- for (j = 0; j < *s; j++) {
+ PROTECT(DNASEQ = coerceVector(DNASEQ, RAWSXP));
+ x = RAW(DNASEQ);
+ n = nrows(DNASEQ);
+ s = ncols(DNASEQ);
+
+ PROTECT(ans = allocVector(INTSXP, s));
+ seg = INTEGER(ans);
+ memset(seg, 0, s*sizeof(int));
+
+ for (j = 0; j < s; j++) {
- i = *n * j; /* start */
- end = i + *n - 1;
+ i = (long) n * j; /* start */
+ end = (long) i + n - 1;
base = x[i];
done = 0;
@@ -1124,6 +1144,9 @@ void SegSites(unsigned char *x, int *n, int *s, int *seg)
i++;
}
}
+
+ UNPROTECT(2);
+ return ans;
}
void GlobalDeletionDNA(unsigned char *x, int *n, int *s, int *keep)
@@ -1181,25 +1204,44 @@ void dist_dna(unsigned char *x, int *n, int *s, int *model, double *d,
}
}
-void C_where(unsigned char *x, unsigned char *pat, int *s, int *p,
- int *ans, int *n)
+SEXP C_where(SEXP DNASEQ, SEXP PAT)
{
- int i, j, k, ln;
-
- ln = 0; /* local n */
-
- for (i = 0; i <= *s - *p; i++) {
+ int p, j, nans;
+ double s, *buf, *a;
+ long i, k;
+ unsigned char *x, *pat;
+ SEXP ans;
+
+ PROTECT(DNASEQ = coerceVector(DNASEQ, RAWSXP));
+ PROTECT(PAT = coerceVector(PAT, RAWSXP));
+ x = RAW(DNASEQ);
+ pat = RAW(PAT);
+ s = XLENGTH(DNASEQ);
+ p = LENGTH(PAT);
+ nans = 0;
+
+ buf = (double *)R_alloc(s, sizeof(double));
+
+ for (i = 0; i <= s - p; i++) {
k = i; j = 0;
while (1) {
if (x[k] != pat[j]) break;
j++; k++;
- if (j == *p) {
- ans[ln++] = k - 1;
+ if (j == p) {
+ buf[nans++] = i + 1;
break;
}
}
}
- *n = ln;
+
+ PROTECT(ans = allocVector(REALSXP, nans));
+ if (nans) {
+ a = REAL(ans);
+ for (i = 0; i < nans; i++) a[i] = buf[i];
+ }
+
+ UNPROTECT(3);
+ return ans;
}
unsigned char codon2aa_Code1(unsigned char x, unsigned char y, unsigned char z)
diff --git a/src/read_dna.c b/src/read_dna.c
index 1edc4f8..c1ef4f0 100644
--- a/src/read_dna.c
+++ b/src/read_dna.c
@@ -1,6 +1,6 @@
-/* read_dna.c 2014-03-05 */
+/* read_dna.c 2016-06-21 */
-/* Copyright 2013-2014 Emmanuel Paradis */
+/* Copyright 2013-2016 Emmanuel Paradis */
/* This file is part of the R-package `ape'. */
/* See the file ../COPYING for licensing issues. */
@@ -81,12 +81,14 @@ static const unsigned char lineFeed = 0x0a;
SEXP rawStreamToDNAbin(SEXP x)
{
- int N, i, j, k, n, startOfSeq;
+ int k, startOfSeq;
+ long i, j, n;
+ double N;
unsigned char *xr, *rseq, *buffer, tmp;
SEXP obj, nms, seq;
PROTECT(x = coerceVector(x, RAWSXP));
- N = LENGTH(x);
+ N = XLENGTH(x);
xr = RAW(x);
/* do a 1st pass to find the number of sequences
@@ -94,19 +96,19 @@ SEXP rawStreamToDNAbin(SEXP x)
this code should be robust to '>' present inside
a label or in the header text before the sequences */
- n = j = 0; /* use j as a flag */
+ n = 0;
+ k = 0; /* use k as a flag */
if (xr[0] == hook) {
- j = 1;
+ k = 1;
startOfSeq = 0;
}
- i = 1;
for (i = 1; i < N; i++) {
- if (j && xr[i] == lineFeed) {
+ if (k && xr[i] == lineFeed) {
n++;
- j = 0;
+ k = 0;
} else if (xr[i] == hook) {
if (!n) startOfSeq = i;
- j = 1;
+ k = 1;
}
}
@@ -114,9 +116,9 @@ SEXP rawStreamToDNAbin(SEXP x)
PROTECT(nms = allocVector(STRSXP, n));
/* Refine the way the size of the buffer is set? */
- buffer = (unsigned char *)R_alloc(N, sizeof(unsigned char *));
+ buffer = (unsigned char *)R_alloc(N, sizeof(unsigned char));
- i = startOfSeq;
+ i = (long) startOfSeq;
j = 0; /* gives the index of the sequence */
while (i < N) {
/* 1st read the label... */
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/r-cran-ape.git
More information about the debian-med-commit
mailing list