[med-svn] [r-cran-ape] 01/02: Imported Upstream version 4.1
Dylan Aïssi
bob.dybian-guest at moszumanska.debian.org
Wed Feb 22 21:41:49 UTC 2017
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 4ef917e38fbca0dad8540c3706b17a2a2ecba350
Author: Dylan Aïssi <bob.dybian at gmail.com>
Date: Wed Feb 22 22:41:29 2017 +0100
Imported Upstream version 4.1
---
DESCRIPTION | 12 +--
MD5 | 40 ++++-----
NAMESPACE | 2 +-
R/DNA.R | 47 ++++++-----
R/apetools.R | 83 ++++++++++++++++++
R/chronos.R | 7 +-
R/drop.tip.R | 83 ++++++------------
R/multi2di.R | 6 +-
R/nj.R | 6 +-
R/node.dating.R | 232 ++++++++++++++++++++++++++++++---------------------
R/read.GenBank.R | 11 ++-
R/read.dna.R | 49 +++++++++--
R/read.nexus.R | 8 +-
inst/doc/MoranI.pdf | Bin 238589 -> 238354 bytes
man/apetools.Rd | 47 +++++++++++
man/as.bitsplits.Rd | 2 +-
man/c.phylo.Rd | 2 +-
man/is.compatible.Rd | 1 +
man/node.dating.Rd | 15 ++--
man/read.dna.Rd | 111 ++++++++++++++++++------
man/read.nexus.Rd | 10 ++-
src/read_dna.c | 11 ++-
22 files changed, 524 insertions(+), 261 deletions(-)
diff --git a/DESCRIPTION b/DESCRIPTION
index c8452d6..c8f74d5 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
Package: ape
-Version: 4.0
-Date: 2016-11-29
+Version: 4.1
+Date: 2017-02-14
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")),
@@ -11,6 +11,7 @@ Authors at R: c(person("Emmanuel", "Paradis", role = c("aut", "cre", "cph"), email
person("Gilles", "Didier", role = c("aut", "cph")),
person("Benoit", "Durand", role = c("aut", "cph")),
person("Julien", "Dutheil", role = c("aut", "cph")),
+ person("RJ", "Ewing", role = c("aut", "cph")),
person("Olivier", "Gascuel", role = c("aut", "cph")),
person("Christoph", "Heibl", role = c("aut", "cph")),
person("Anthony", "Ives", role = c("aut", "cph")),
@@ -32,11 +33,11 @@ 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 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 [...]
+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-11-30 09:33:29 UTC; paradis
+Packaged: 2017-02-14 16:20:43 UTC; paradis
Author: Emmanuel Paradis [aut, cre, cph],
Simon Blomberg [aut, cph],
Ben Bolker [aut, cph],
@@ -46,6 +47,7 @@ Author: Emmanuel Paradis [aut, cre, cph],
Gilles Didier [aut, cph],
Benoit Durand [aut, cph],
Julien Dutheil [aut, cph],
+ RJ Ewing [aut, cph],
Olivier Gascuel [aut, cph],
Christoph Heibl [aut, cph],
Anthony Ives [aut, cph],
@@ -64,4 +66,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-12-01 17:40:32
+Date/Publication: 2017-02-14 19:15:21
diff --git a/MD5 b/MD5
index 6a8e3b2..e958f98 100644
--- a/MD5
+++ b/MD5
@@ -1,12 +1,12 @@
eb723b61539feef013de476e68b5c50a *COPYING
-e28f0a5656abc88ffa75eb720a68e0ff *DESCRIPTION
-6936565af4821e47088568525c3bda01 *NAMESPACE
+04fa612328ddf8b15a99cc4f6eb2cf5a *DESCRIPTION
+c08ae7d7f93544183572e7c378f67c74 *NAMESPACE
854a025cb7e5da3e4fe230c0be950d08 *NEWS
0c7bc9101516fd26fb3ddbedbe25b6a3 *R/CADM.global.R
e042efca1d8a00d4178204f5f481b64d *R/CADM.post.R
a10ebb0e9ecf861eea483eb41da1bd23 *R/CDF.birth.death.R
fdb9cfa0cbda82bda982b290693e44e3 *R/Cheverud.R
-32741db1a4135fdd540d71a103d78d6f *R/DNA.R
+ba0115131695f6fb2b7fb58f4d63300e *R/DNA.R
2187a1289b767066d1efe1ebbe7c3b0c *R/MPR.R
bb95af56d882b6162aa517a77140175e *R/MoranI.R
218b8cf3c13a700c757b7b2303dc897e *R/PGLS.R
@@ -16,6 +16,7 @@ bb95af56d882b6162aa517a77140175e *R/MoranI.R
4ce79cf3f3ff49bef989454d86d0c891 *R/additive.R
9fca35b1005cce7b2a260b553dd971a3 *R/alex.R
9fe874382f024a98f62a0ccfcd6d09ac *R/all.equal.phylo.R
+aa844e5d40f27d3dbbe45cd932b39ad6 *R/apetools.R
06b049b681ab2e594124a2b010c5cddd *R/as.bitsplits.R
c94018d5e792c72e20ce84085b2df9e7 *R/as.matching.R
86d55d04b20c1ecad73ea817866d632d *R/as.phylo.R
@@ -30,7 +31,7 @@ b4557270916972166e13e76391bb0a69 *R/checkValidPhylo.R
e43b5dec7eae6d4bf9371e50117bf6ed *R/cherry.R
5fd65efbdae6e0934e8bb42191142d4b *R/chronoMPL.R
74e1019810b06458e808a447bb099a91 *R/chronopl.R
-8ae1af4a30c40d15676d70e180f3d593 *R/chronos.R
+f68d5d19e6a080ff70e03f41b7e43400 *R/chronos.R
2a9983307f686facf6d9ec843d214f38 *R/clustal.R
dedf66f0595977a11a7dbc2b0d2f56bb *R/coalescent.intervals.R
c3d44d960238a711764aaeebe9570559 *R/collapse.singles.R
@@ -49,7 +50,7 @@ dfd5bb35f1cb1fd9154d023e0e4cfc2b *R/dist.gene.R
c59b82a465d826b821896e0da82f1017 *R/dist.topo.R
b28ced504fedeb7f991f7eba10ad06df *R/diversi.gof.R
8b2ec4004022afdc7e2cb42f2657b628 *R/diversi.time.R
-ea5fcbd0f207ab1d874c0d032392967d *R/drop.tip.R
+d83a976aa9a582ab31ea12f65d7c6374 *R/drop.tip.R
1c4bdd6cf163aa4d605ed67249883ad4 *R/evonet.R
fceafc86fae624fd9037403ad301d35a *R/ewLasso.R
aa09abeb90ef891384128f978ffce843 *R/extract.popsize.R
@@ -71,11 +72,11 @@ d2c16632492bfafd2ee18f2fe3d3d64a *R/matexpo.R
61021f7af1175c46a743c7fee4cdc87e *R/me.R
e1327a592199bac7ac8f670c3a068d46 *R/mrca.R
a078728fb5907565f85b54b30e5bf83f *R/mst.R
-93f3fe244219d565bb2e111ff249b540 *R/multi2di.R
+dea8394974d847a254f810883e1f246e *R/multi2di.R
0850fdd19c01d37ac632fc308632a463 *R/mvr.R
-282308c27ac1e78e330711d635d2e572 *R/nj.R
+2f302fa091736fd5b1cb067c6ecd1215 *R/nj.R
e3f22d0f260c43be87a17b0ab091e2bb *R/njs.R
-5e43665c088ced872dd2a03f7adc2add *R/node.dating.R
+3efe5029d25c09e238b3d45aa9b3a140 *R/node.dating.R
e4a0712fa45754cd02f5d3da6ef23204 *R/nodelabels.R
ae2aeb0e8aef7f8d4b19939ca61b3482 *R/nodepath.R
d9fd8e402e6fce6939a05332823a9390 *R/parafit.R
@@ -90,11 +91,11 @@ e579ec65c808c29e1ecaae1501784b37 *R/plot.popsize.R
1e2485437566ca9af99d93b4580cbbc2 *R/print.lmorigin.R
d0e8bd41d5acc217fdee3578adcf635b *R/print.parafit.R
8c401518738b9cda403fa9f0eb382757 *R/rTrait.R
-9c261f833a2f3325949b00f49a3d9d61 *R/read.GenBank.R
+d05e33bba295c6208e302275668b3a27 *R/read.GenBank.R
b13dfb8f455b1c9e74a364085f72dbce *R/read.caic.R
-740c43ccf599f0ea029fddf4361f5738 *R/read.dna.R
+ec0d133f24b736de2d3a66edad978f83 *R/read.dna.R
b25674d229dfa4d2bf76a50e1874c690 *R/read.gff.R
-12cc225d7dbe722e9e61cd3911de0b5a *R/read.nexus.R
+d8045fd7da3d9aabaeb11216098995d0 *R/read.nexus.R
13ce7f5c7d1bcb7101469d12651e99c8 *R/read.nexus.data.R
88e633c80886f5a6c9f2b34918452d15 *R/read.tree.R
df2ac5d3de7185c7e26fc95b35192a40 *R/reconstruct.R
@@ -145,7 +146,7 @@ e2d1339025ed901009bfed58dc6505ff *data/mat5M3ID.RData
828290996b613493f96e0fab024beebb *inst/CITATION
3f54f3775bcf382e25df2a12228894f6 *inst/doc/MoranI.R
2277b0efdb2f70dfdc5df8278702c6b6 *inst/doc/MoranI.Rnw
-31ca3eafda127eab7c4799764f8fa68c *inst/doc/MoranI.pdf
+49be1aafbc668a68960f23c66d6cf61a *inst/doc/MoranI.pdf
db2ed5302bddcfac232dcd6fabe4a651 *man/AAbin.Rd
e6876b193a0df06697c788a8e48cf4bc *man/CADM.global.Rd
dfa15a3e3bb57c9b21a30b8d5d790c62 *man/DNAbin.Rd
@@ -164,8 +165,9 @@ f0b85d459edb5ead2e6d39a86cbb350c *man/all.equal.phylo.Rd
c3bd9e7b02831d6086fabd5de6db917c *man/alview.Rd
42489c977f59031dc17a0ba7e55663b6 *man/ape-internal.Rd
f84e59c805efa00000d775b81bf9f8d0 *man/ape-package.Rd
+d051129f1de68268f43d1360601b9c25 *man/apetools.Rd
5bba4ae4bfc66b613855cfc182d9b1bc *man/as.alignment.Rd
-2771f272ad46b279022892d9f5d21eb2 *man/as.bitsplits.Rd
+4eed5aecc55738e6bdda6839b91cc635 *man/as.bitsplits.Rd
4f014cf2923e2eab6188acd48e8096fa *man/as.matching.Rd
78ff1c24cc3b91de1bf33ca5f45c7af5 *man/as.phylo.Rd
219cff7b94e167e18f4f9cd99cc8efc3 *man/as.phylo.formula.Rd
@@ -182,7 +184,7 @@ f929bc1b6391c57a6b0099c4561fd7be *man/binaryPGLMM.Rd
ef1c15d5d93410c21179997431112209 *man/birthdeath.Rd
87fafa0226908ea87c6d11aecc9bf72c *man/boot.phylo.Rd
5a64b90d3a6c7a8204946b00f45f4cfc *man/branching.times.Rd
-c808768a942998429ecc9fc925f3f3a3 *man/c.phylo.Rd
+99ffa532ab4397c374eaddd0f2ff8469 *man/c.phylo.Rd
9d3f9614732671a610cc8a37454885e2 *man/carnivora.Rd
953b3c21c60a38e8e53ec6894170dd36 *man/checkAlignment.Rd
5ff8c7e8fad519d978f166948c03059c *man/checkValidPhylo.Rd
@@ -234,7 +236,7 @@ eea313e8ee32597b4cec120d23113642 *man/gammaStat.Rd
86c49d080fdffd614d8056021e91cc55 *man/identify.phylo.Rd
dc00b2b43a44f4ba7ea76ea23a12c615 *man/image.DNAbin.Rd
0ede757d7bed82216980d3c4fd592cb6 *man/is.binary.tree.Rd
-80a5a228a43679edf75903590253a875 *man/is.compatible.Rd
+ab3d6e53acd7fe19bdf48a12912f005f *man/is.compatible.Rd
d2de8fd9549ef01a1dddeb726dd77fcf *man/is.monophyletic.Rd
ad9e7316219c3238b44b02d55c44b4d3 *man/is.ultrametric.Rd
3f6ff8340b6c9770a1f4d2fed72a295d *man/kronoviz.Rd
@@ -259,7 +261,7 @@ b74a379a0a9b80941118922ca7843b02 *man/multi2di.Rd
00fb7ade93c2dd887be27e3dad8e2115 *man/mvr.Rd
3df9e16b8a09df3f1dba5c4327a635fc *man/nj.Rd
9ea7d5899a190171de4165e3240de62e *man/njs.Rd
-8a9189a62be3d0e871deb20be5bd6e8a *man/node.dating.Rd
+1f2bbc4477c7fcfad397a6bdfa11faee *man/node.dating.Rd
a589b4cc04505185dc9ef1817c7ae304 *man/node.depth.Rd
ba754eba67181a1c39d417ffd9d23b76 *man/nodelabels.Rd
447ae03684ff56a4a24932aec182acf3 *man/nodepath.Rd
@@ -280,9 +282,9 @@ b24438c42cea969302ec6ba61002426e *man/print.phylo.Rd
81f756fdf2ec4c968095595dece8338f *man/rTraitMult.Rd
028b11582b3493cdefea60ccb78ad429 *man/read.GenBank.Rd
9e31bccfa08c1d2da2180cd8f50ab250 *man/read.caic.Rd
-2b32811be5bed71ffeca3801c4408e5b *man/read.dna.Rd
+87c739eafaaa091de24d8dc0c47248b8 *man/read.dna.Rd
51670875393d70cb7ef84b4b0388e8e4 *man/read.gff.Rd
-5c9aed0b6f0a7e0f0fac9b5a59ac732f *man/read.nexus.Rd
+fdc2393e70328ae27034c54bf8a847c7 *man/read.nexus.Rd
bc02e36c51d67074e661468993ed359b *man/read.nexus.data.Rd
3c8b53f5af4fdd36e7c0b449412d5a77 *man/read.tree.Rd
9224f164351fea4be9f7ce9d827f8e58 *man/reconstruct.Rd
@@ -352,7 +354,7 @@ d3364013ad597425d01ac50ae4ee0b34 *src/nj.c
7de13de6d828d870915c2efb6eb47ab6 *src/pic.c
28ef72b2ccede817c0431afcb04e9c8c *src/plot_phylo.c
aa8d9966da3b7e970879a49a40c05a07 *src/rTrait.c
-23d15137e1f3d9b0d365d04c1f4d9b05 *src/read_dna.c
+3789c732cdd750e3c5e813884968cf03 *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 5f76a37..03219a6 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -23,7 +23,7 @@ importFrom(stats, AIC, as.dist, as.hclust, biplot, coef, complete.cases,
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(utils, download.file, edit, read.table, str)
importFrom(parallel, mclapply)
## Methods for the classes defined in ape, including for the generics
diff --git a/R/DNA.R b/R/DNA.R
index 0757da6..010f61d 100644
--- a/R/DNA.R
+++ b/R/DNA.R
@@ -1,8 +1,8 @@
-## DNA.R (2016-11-26)
+## DNA.R (2017-01-12)
## Manipulations and Comparisons of DNA and AA Sequences
-## Copyright 2002-2016 Emmanuel Paradis, 2015 Klaus Schliep
+## Copyright 2002-2017 Emmanuel Paradis, 2015 Klaus Schliep
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
@@ -213,15 +213,13 @@ print.DNAbin <- function(x, printlen = 6, digits = 3, ...)
if (n == 1) {
cat("1 DNA sequence in binary format stored in a list.\n\n")
nTot <- length(x[[1]])
- cat("Sequence length:", nTot, "\n\n")
- cat("Label:", nms, "\n\n")
+ cat("Sequence length:", nTot, "\n")
} else {
cat(n, "DNA sequences in binary format stored in a list.\n\n")
tmp <- lengths(x, use.names = FALSE)
nTot <- sum(as.numeric(tmp))
- mini <- range(tmp)
- maxi <- mini[2]
- mini <- mini[1]
+ mini <- min(tmp)
+ maxi <- max(tmp)
if (mini == maxi)
cat("All sequences of same length:", maxi, "\n")
else {
@@ -229,35 +227,40 @@ print.DNAbin <- function(x, printlen = 6, digits = 3, ...)
cat(" Shortest sequence:", mini, "\n")
cat(" Longest sequence:", maxi, "\n")
}
- TAIL <- "\n\n"
- if (printlen < n) {
- nms <- nms[1:printlen]
- TAIL <- "...\n\n"
- }
- cat("\nLabels:", paste(nms, collapse = " "), TAIL)
}
} else {
nTot <- length(x)
if (is.matrix(x)) {
nd <- dim(x)
+ n <- nd[1]
nms <- rownames(x)
- cat(nd[1], "DNA sequences in binary format stored in a matrix.\n\n")
- cat("All sequences of same length:", nd[2], "\n")
- TAIL <- "\n\n"
- if (printlen < nd[1]) {
- nms <- nms[1:printlen]
- TAIL <- "...\n\n"
+ if (n == 1) {
+ cat("1 DNA sequence in binary format stored in a matrix.\n\n")
+ cat("Sequence length:", nd[2], "\n")
+ } else {
+ cat(n, "DNA sequences in binary format stored in a matrix.\n\n")
+ cat("All sequences of same length:", nd[2], "\n")
}
- cat("\nLabels:", paste(nms, collapse = " "), TAIL)
} else {
cat("1 DNA sequence in binary format stored in a vector.\n\n")
cat("Sequence length:", nTot, "\n\n")
}
}
- if (nTot <= 1e6) {
+ if (exists("nms")) {
+ HEAD <- if (n == 1) "\nLabel:" else "\nLabels:"
+ TAIL <- ""
+ if (printlen < n) {
+ nms <- nms[1:printlen]
+ TAIL <- "...\n"
+ }
+ if (any(longs <- nchar(nms) > 60))
+ nms[longs] <- paste0(substr(nms[longs], 1, 60), "...")
+ cat(HEAD, nms, TAIL, sep = "\n")
+ }
+ if (nTot <= 1e7) {
cat("Base composition:\n")
print(round(base.freq(x), digits))
- } else cat("More than 1 million nucleotides: not printing base composition\n")
+ } else cat("More than 10 million nucleotides: not printing base composition\n")
}
as.DNAbin <- function(x, ...) UseMethod("as.DNAbin")
diff --git a/R/apetools.R b/R/apetools.R
new file mode 100644
index 0000000..090498f
--- /dev/null
+++ b/R/apetools.R
@@ -0,0 +1,83 @@
+## apetools.R (2017-02-03)
+
+## APE Tools
+
+## Copyright 2005-2017 Emmanuel Paradis
+
+## This file is part of the R-package `ape'.
+## See the file ../COPYING for licensing issues.
+
+.file.extensions <-
+ list(clustal = "aln", fasta = c("fas", "fasta"),
+ fastq = c("fq", "fastq"), newick = c("nwk", "newick", "tre", "tree"),
+ nexus = c("nex", "nexus"), phylip = "phy")
+
+Xplorefiles <- function(from = "HOME", recursive = TRUE, ignore.case = TRUE)
+{
+ if (from == "HOME") from <- Sys.getenv("HOME")
+ FILES <- list.files(path = from, recursive = recursive, full.names = TRUE)
+ ext <- if (exists(".file.extensions", envir = .PlotPhyloEnv))
+ get(".file.extensions", envir = .PlotPhyloEnv)
+ else .file.extensions
+ res <- vector("list", length(ext))
+ names(res) <- names(ext)
+ for (i in seq_along(res)) {
+ e <- paste0("\\.", ext[[i]], "$")
+ if (length(e) > 1) e <- paste(e, collapse = "|")
+ x <- grep(e, FILES, ignore.case = ignore.case, value = TRUE)
+ res[[i]] <- data.frame(File = x, Size = file.size(x),
+ stringsAsFactors = FALSE)
+ }
+ res
+}
+
+editFileExtensions <- function()
+{
+ foo <- function(x) {
+ n <- length(x)
+ if (n < m) x[(n + 1):m] <- NA
+ x
+ }
+ res <- if (exists(".file.extensions", envir = .PlotPhyloEnv))
+ get(".file.extensions", envir = .PlotPhyloEnv)
+ else .file.extensions
+ m <- max(lengths(res, FALSE))
+ res <- lapply(res, foo)
+ res <- as.data.frame(res, stringsAsFactors = FALSE)
+ res <- edit(res)
+ res <- lapply(res, function(x) x[!is.na(x)])
+ assign(".file.extensions", res, envir = .PlotPhyloEnv)
+}
+
+bydir <- function(x)
+{
+ nofile <- which(sapply(x, nrow) == 0)
+ if (length(nofile)) x <- x[-nofile]
+ if (!length(x)) {
+ cat("No file\n")
+ return(invisible(NULL))
+ }
+ for (i in seq_along(x)) x[[i]]$Type <- names(x)[i]
+ x <- do.call(rbind, x)
+ x <- x[order(x$File), ]
+ SPLIT <- strsplit(x$File, "/")
+ LL <- lengths(SPLIT)
+ foo <- function(i, PATH) {
+ K <- grep(paste0("^", PATH, "/"), x$File)
+ sel <- intersect(K, which(LL == i + 1L))
+ if (length(sel)) {
+ y <- x[sel, ]
+ y$File <- gsub(".*/", "", y$File)
+ cat("\n", PATH, "/\n", sep = "")
+ print(y, row.names = FALSE)
+ }
+ if (length(sel) < length(K)) {
+ d <- setdiff(K, sel)
+ subdir <- unlist(lapply(SPLIT[d], "[", i + 1L))
+ for (z in unique(subdir))
+ foo(i + 1L, paste(PATH, z, sep = "/"))
+ }
+ }
+ top <- unlist(lapply(SPLIT, "[", 1L))
+ for (z in unique(top)) foo(1L, z)
+}
diff --git a/R/chronos.R b/R/chronos.R
index a7f63a9..3483bd2 100644
--- a/R/chronos.R
+++ b/R/chronos.R
@@ -1,8 +1,8 @@
-## chronos.R (2014-05-15)
+## chronos.R (2016-12-7)
## Molecular Dating With Penalized and Maximum Likelihood
-## 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.
@@ -79,6 +79,7 @@ chronos <-
ROOT <- n + 1L
m <- phy$Nnode
el <- phy$edge.length
+ if (is.null(el)) stop("the tree has no branch lengths")
if (any(el < 0)) stop("some branch lengths are negative")
e1 <- phy$edge[, 1L]
e2 <- phy$edge[, 2L]
@@ -437,8 +438,6 @@ maybe you need to adjust the calibration dates")
if (!quiet) cat("\nDone.\n")
-# browser()
-
if (model == "discrete") {
rate.freq <-
if (Nb.rates == 1) current.rates
diff --git a/R/drop.tip.R b/R/drop.tip.R
index 62bd65c..5187919 100644
--- a/R/drop.tip.R
+++ b/R/drop.tip.R
@@ -1,19 +1,19 @@
-## drop.tip.R (2015-07-28)
+## drop.tip.R (2017-02-14)
## Remove Tips in a Phylogenetic Tree
-## Copyright 2003-2015 Emmanuel Paradis
+## Copyright 2003-2017 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
extract.clade <- function(phy, node, root.edge = 0, interactive = FALSE)
{
- Ntip <- length(phy$tip.label)
- ROOT <- Ntip + 1
- Nedge <- dim(phy$edge)[1]
- wbl <- !is.null(phy$edge.length)
- if (interactive) node <- identify(phy)$nodes else {
+ n <- length(phy$tip.label)
+ if (interactive) {
+ cat("Click close to the node...\n")
+ node <- identify(phy)$nodes
+ } else {
if (length(node) > 1) {
node <- node[1]
warning("only the first value of 'node' has been considered")
@@ -21,60 +21,15 @@ extract.clade <- function(phy, node, root.edge = 0, interactive = FALSE)
if (is.character(node)) {
if (is.null(phy$node.label))
stop("the tree has no node labels")
- node <- which(phy$node.label %in% node) + Ntip
+ node <- match(node, phy$node.label) + n
+ if (is.na(node)) stop("'node' not among the node labels.")
}
- if (node <= Ntip)
+ if (node <= n)
stop("node number must be greater than the number of tips")
}
- if (node == ROOT) return(phy)
- phy <- reorder(phy) # insure it is in cladewise order
- root.node <- which(phy$edge[, 2] == node)
- start <- root.node + 1 # start of the clade looked for
- anc <- phy$edge[root.node, 1] # the ancestor of 'node'
- next.anc <- which(phy$edge[-(1:start), 1] <= anc) # find the next occurence of 'anc' or an 'older' node
-
- keep <- if (length(next.anc)) start + 0:(next.anc[1] - 1) else start:Nedge
-
- if (root.edge) {
- NewRootEdge <- phy$edge.length[root.node]
- root.edge <- root.edge - 1
- while (root.edge) {
- if (anc == ROOT) break
- i <- which(phy$edge[, 2] == anc)
- NewRootEdge <- NewRootEdge + phy$edge.length[i]
- root.edge <- root.edge - 1
- anc <- phy$edge[i, 1]
- }
- if (root.edge && !is.null(phy$root.edge))
- NewRootEdge <- NewRootEdge + phy$root.edge
- phy$root.edge <- NewRootEdge
- }
-
- phy$edge <- phy$edge[keep, ]
- if (wbl) phy$edge.length <- phy$edge.length[keep]
- TIPS <- phy$edge[, 2] <= Ntip
- tip <- phy$edge[TIPS, 2]
- ## Fix by Ludovic Mallet and Mahendra Mariadassou (2011-11-21):
- name <- vector("character", length(tip))
- name[order(tip)] <- phy$tip.label[tip]
- phy$tip.label <- name
- ## End of fix
- ## keep the ordering so no need to reorder tip.label:
- phy$edge[TIPS, 2] <- order(tip)
- if (!is.null(phy$node.label))
- phy$node.label <- phy$node.label[sort(unique(phy$edge[, 1])) - Ntip]
- Ntip <- length(phy$tip.label)
- phy$Nnode <- dim(phy$edge)[1] - Ntip + 1L
- ## The block below renumbers the nodes so that they conform
- ## to the "phylo" format -- same as in root()
- newNb <- integer(Ntip + phy$Nnode)
- newNb[node] <- Ntip + 1L
- sndcol <- phy$edge[, 2] > Ntip
- ## executed from right to left, so newNb is modified before phy$edge:
- phy$edge[sndcol, 2] <- newNb[phy$edge[sndcol, 2]] <-
- (Ntip + 2):(Ntip + phy$Nnode)
- phy$edge[, 1] <- newNb[phy$edge[, 1]]
- phy
+ if (node == n + 1L) return(phy)
+ keep <- prop.part(phy)[[node - n]]
+ drop.tip(phy, (1:n)[-keep], root.edge = root.edge, rooted = TRUE)
}
drop.tip <-
@@ -116,6 +71,18 @@ drop.tip <-
wbl <- !is.null(phy$edge.length)
+ if (length(tip) == Ntip - 1 && trim.internal) {
+ i <- which(phy$edge[, 2] == (1:Ntip)[-tip])
+ res <- list(edge = matrix(2:1, 1, 2),
+ tip.label = phy$tip.label[phy$edge[i, 2]],
+ Nnode = 1L)
+ class(res) <- "phylo"
+ if (wbl) res$edge.length <- phy$edge.length[i]
+ if (!is.null(phy$node.label))
+ res$node.label <- phy$node.label[phy$edge[i, 1] - Ntip]
+ return(res)
+ }
+
if (!rooted && subtree) {
phy <- root(phy, (1:Ntip)[-tip][1])
root.edge <- 0
diff --git a/R/multi2di.R b/R/multi2di.R
index 364b0ea..2e79eb9 100644
--- a/R/multi2di.R
+++ b/R/multi2di.R
@@ -1,8 +1,8 @@
-## multi2di.R (2016-10-16)
+## multi2di.R (2017-01-16)
## Collapse or Resolve Multichotomies
-## Copyright 2005-2016 Emmanuel Paradis
+## Copyright 2005-2017 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
@@ -69,7 +69,7 @@ multi2di <- function(phy, ...) UseMethod("multi2di")
if (!is.null(phy$node.label))
phy$node.label <-
c(phy$node.label, rep("", phy$Nnode - length(phy$node.label)))
- phy <- reorder(phy)
+ phy <- .reorder_ape(phy, "cladewise", FALSE, n, 1L) # fix by Klaus (2017-01-16)
## the node numbers are not in increasing order in edge[, 2]: this
## will confuse drop.tip and other functions (root), so renumber them
diff --git a/R/nj.R b/R/nj.R
index 3a49f45..7efc4ca 100644
--- a/R/nj.R
+++ b/R/nj.R
@@ -1,8 +1,8 @@
-## nj.R (2009-11-23)
+## nj.R (2017-01-24)
## Neighbor-Joining Tree Estimation
-## Copyright 2004-2009 Emmanuel Paradis
+## Copyright 2004-2017 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
@@ -12,6 +12,8 @@ nj <- function(X)
if (is.matrix(X)) X <- as.dist(X)
if (any(is.na(X)))
stop("missing values are not allowed in the distance matrix\nConsider using njs()")
+ if (any(is.infinite(X)))
+ stop("infinite values are not allowed in the distance matrix")
N <- attr(X, "Size")
labels <- attr(X, "Labels")
if (is.null(labels)) labels <- as.character(1:N)
diff --git a/R/node.dating.R b/R/node.dating.R
index bafe3cf..4152263 100644
--- a/R/node.dating.R
+++ b/R/node.dating.R
@@ -1,6 +1,6 @@
## node.dating.R (2016-06-21)
## This file is part of the R-package `ape'.
-## See the file ../COPYING for licensing issues.
+## See the file COPYING in the package ape available at cran.r-project.org for licensing issues.
# Copyright (c) 2016, Bradley R. Jones, BC Centre for Excellence in HIV/AIDS
# All rights reserved.
@@ -38,8 +38,7 @@
# 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.
+# clock.
#
# t: rooted tree with edge lengths equal to genetic distance
#
@@ -49,16 +48,16 @@
# 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) {
+estimate.mu <- function(t, node.dates, p.tol=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=""))
+ if ((1 - pchisq(AIC(null.g) - AIC(g) + 2, df=1)) > p.tol) {
+ warning(paste("Cannot reject null hypothesis (p=", (1 - pchisq(AIC(null.g) - AIC(g))), ")"))
}
-
+
coef(g)[[2]]
}
@@ -69,7 +68,8 @@ estimate.mu <- function(t, node.dates, p = 0.05) {
# 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
+# mu: mutation rate, either a vector of size one for a strict molecular clock
+# or a vector with a local molecular clock along each edge
#
# min.date: the minimum date that a node can have (needed for optimize()). The
# default is -.Machine$double.xmax
@@ -94,121 +94,167 @@ estimate.mu <- function(t, node.dates, p = 0.05) {
# 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)
+# returns a vector of the estimated dates of the tips and internal 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 = 0, is.binary = is.binary.phylo(t)) {
+
+ # check parameters
+ if (any(mu < 0))
stop(paste("mu (", mu, ") less than 0", sep=""))
-
- # init vars
+
+ # init vars
+ mu <- if (length(mu) == 1) rep(mu, length(t$edge.length)) else mu
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
}
+
+ lik.sens <- if (lik.tol == 0) opt.tol else lik.tol
# Don't count initial step if all values are seeded
- iter.step <- if (any(is.na(dates))) {
- 0
- } else {
- 1
- }
+ iter.step <- if (any(is.na(dates))) 0 else 1
children <- lapply(1:t$Nnode,
function(x) {
- t$edge[,1] == x + n.tips
+ which(t$edge[,1] == x + n.tips)
})
parent <- lapply(1:t$Nnode,
function(x) {
- t$edge[,2] == x + n.tips
+ which(t$edge[,2] == x + n.tips)
})
# to process children before parents
nodes <- c(1)
- i <- 1
+ for (i in 1:t$Nnode) {
+ to.add <- t$edge[children[[nodes[i]]], 2] - n.tips
- while (i <= length(nodes)) {
- nodes <- c(nodes, t$edge[t$edge[,1] == nodes[i] + n.tips, 2] - n.tips)
+ nodes <- c(nodes, to.add[to.add > 0])
i <- i + 1
}
-
- nodes <- nodes[nodes > 0]
nodes <- rev(nodes)
- # calculate likelihood
+ # calculate likelihood functions
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
+
+ t$edge.length[ch.edge]*log(tim)-mu[ch.edge]*tim
}
-
- opt.fun <- function(x, ch, p, ch.edge.length, p.edge.length, use.parent=T) {
+
+ opt.fun <- function(x, ch, p, ch.edge, p.edge, use.parent=T) {
sum(if (!use.parent || length(dates[p]) == 0 || is.na(dates[p])) {
- calc.Like(dates[ch], ch.edge.length, x)
+ calc.Like(dates[ch], t$edge.length[ch.edge], x)
} else {
- calc.Like(c(dates[ch], x), c(ch.edge.length, p.edge.length), c(rep(x, length(dates[ch])), dates[p]))
+ calc.Like(c(dates[ch], x), c(t$edge.length[ch.edge], t$edge.length[p.edge]), 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]
+ solve.lin <- function(bounds, ch.times, ch.edge) {
+ y <- (mu[ch.edge] * ch.times - t$edge.length[ch.edge]) / mu[ch.edge]
+ x <- c(bounds[1] + opt.tol, bounds[2] - opt.tol)
+ if (bounds[1] < y && y < bounds[2])
+ x <- c(x, y)
+
+ x[which.max(unlist(lapply(x, function(y) sum(calc.Like(ch.times, ch.edge, y)))))]
+ }
+
+ solve.poly2 <- function(bounds, a, b, c.0) {
+ x <- c(bounds[1] + opt.tol, bounds[2] - opt.tol)
+
+ if (b ^ 2 - 4 * a * c.0 >= 0) {
+ if (a == 0) {
+ y <- -c.0 / b
+
+ if (bounds[1] < y && y < bounds[2])
+ x <- c(x, y)
+ } 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)
+
+ 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)
+ }
+ }
+
+ x
+ }
+
+ solve.bin <- function(bounds, ch.times, ch.edge) {
+ ch.edge.length <- t$edge.length[ch.edge]
+ a <- sum(mu[ch.edge])
+ b <- ch.edge.length[1] + ch.edge.length[2] - a * (ch.times[1] + ch.times[2])
+ c.0 <- a*ch.times[1] * ch.times[2] - ch.times[1] * ch.edge.length[2] - ch.times[2] * ch.edge.length[1]
+
+ x <- solve.poly2(bounds, a, b, c.0)
+
+ x[which.max(unlist(lapply(x, function(y) sum(calc.Like(ch.times, ch.edge, y)))))]
+ }
+
+
+ solve.bin2 <- function(bounds, ch.times, ch.edge, par.time, par.edge) {
+ ch.edge.length <- t$edge.length[ch.edge]
+ par.edge.length <- t$edge.length[par.edge]
+ a <- mu[ch.edge] - mu[par.edge]
+ b <- ch.edge.length + par.edge.length - a * (ch.times + par.time)
+ c.0 <- a*ch.times * par.time - ch.times * par.edge.length - par.time * ch.edge.length
- b ^ 2 - 4 * a * c.0 < 0
+ cat(sprintf("a: %f, b: %f, c: %f\n", a, b, c.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 <- solve.poly2(bounds, a, b, c.0)
- x <- c(bounds[1] + opt.tol, bounds[2] - opt.tol)
- if (bounds[1] < x.1 && x.1 < bounds[2])
+ x[which.max(unlist(lapply(x, function(y) sum(calc.Like(c(ch.times, y), c(ch.edge, par.edge), c(y, par.time))))))]
+ }
+
+ solve.poly3 <- function(bounds, a, b, c.0, d) {
+ x <- c(bounds[1] + opt.tol, bounds[2] - opt.tol)
+
+ if (a == 0)
+ x <- c(x, solve.poly2(bounds, b, c.0, d))
+ else {
+ 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)))
+
+ 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)
-
- return(x[which.max(unlist(lapply(x, function(y) sum(calc.Like(ch.times, ch.edge.length, y)))))])
+ if (bounds[1] < x.3 && x.3 < bounds[2])
+ x <- c(x, x.3)
}
+
+ x
}
- 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)))
+ solve.cube <- function(bounds, ch.times, ch.edge, par.time, par.edge) {
+ ch.edge.length <- t$edge.length[ch.edge]
+ par.edge.length <- t$edge.length[par.edge]
+
+ a <- sum(mu[ch.edge]) - mu[par.edge]
+ b <- sum(ch.edge.length) + par.edge.length - a * (sum(ch.times) + par.time)
+ c.0 <- a * (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] - a * prod(ch.times) * par.time
- 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)
+ x <- solve.poly3(bounds, a, b, c.0, d)
- 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))))))])
+ x[which.max(unlist(lapply(x, function(y) sum(calc.Like(c(ch.times, y), c(ch.edge, par.edge), 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]]]
+ ch.edge <- children[[node]]
+ ch <- t$edge[ch.edge, 2]
+
+ p.edge <- parent[[node]]
+ p <- t$edge[p.edge, 1]
m <- if (length(p) == 0 || is.na(dates[p])) {
min.date
@@ -217,13 +263,19 @@ estimate.dates <- function(t, node.dates, mu = estimate.mu(t, node.dates), min.d
}
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)
+ if (length(dates[p]) == 0 || is.na(dates[p])) {
+ if (length(ch.edge) == 2)
+ solve.bin(c(m, min(dates[ch])), dates[ch], ch.edge)
+ else
+ solve.lin(c(m, min(dates[ch])), dates[ch], ch.edge)
+ } else {
+ if (length(ch.edge) == 2)
+ solve.cube(c(m, min(dates[ch])), dates[ch], ch.edge, dates[p], p.edge)
+ else
+ solve.bin2(c(m, min(dates[ch])), dates[ch], ch.edge, dates[p], p.edge)
+ }
+ } else {
+ res <- optimize(opt.fun, c(m, min(dates[ch])), ch, p, ch.edge, p.edge, maximum=T)
res$maximum
}
@@ -238,15 +290,8 @@ estimate.dates <- function(t, node.dates, mu = estimate.mu(t, node.dates), min.d
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
+ all.lik <- calc.Like(dates[t$edge[,2]], 1:length(t$edge.length), dates[t$edge[,1]]) + scale.lik
+ new.lik <- sum(all.lik)
if (show.steps > 0 && ((iter.step %% show.steps) == 0)) {
cat(paste("Step: ", iter.step, ", Likelihood: ", new.lik, "\n", sep=""))
@@ -256,7 +301,7 @@ estimate.dates <- function(t, node.dates, mu = estimate.mu(t, node.dates), min.d
if (is.infinite(lik) || is.infinite(new.lik)) {
warning("Likelihood infinite")
}
- else if (!is.na(lik) && new.lik < lik - lik.tol) {
+ else if (!is.na(lik) && new.lik + lik.sens < lik) {
warning("Likelihood less than previous estimate")
}
@@ -274,4 +319,3 @@ estimate.dates <- function(t, node.dates, mu = estimate.mu(t, node.dates), min.d
dates
}
-
diff --git a/R/read.GenBank.R b/R/read.GenBank.R
index d2a44c2..8115d7a 100644
--- a/R/read.GenBank.R
+++ b/R/read.GenBank.R
@@ -1,8 +1,8 @@
-## read.GenBank.R (2016-07-27)
+## read.GenBank.R (2017-01-18)
## Read DNA Sequences from GenBank via Internet
-## Copyright 2002-2016 Emmanuel Paradis
+## Copyright 2002-2017 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
@@ -26,8 +26,13 @@ read.GenBank <- function(access.nb, seq.names = access.nb, species.names = TRUE,
if (b > N) b <- N
}
res <- read.FASTA(fl)
+ if (is.null(res)) return(NULL)
attr(res, "description") <- names(res)
- names(res) <- access.nb
+ if (length(access.nb) != length(res)) {
+ names(res) <- gsub("\\..*$", "", names(res))
+ failed <- paste(access.nb[! access.nb %in% names(res)], collapse = ", ")
+ warning(paste0("cannot get the following sequences:\n", failed))
+ } else names(res) <- access.nb
if (as.character) res <- as.character(res)
diff --git a/R/read.dna.R b/R/read.dna.R
index 79ce558..9ec231c 100644
--- a/R/read.dna.R
+++ b/R/read.dna.R
@@ -1,28 +1,42 @@
-## read.dna.R (2016-06-21)
+## read.dna.R (2017-02-03)
## Read DNA Sequences in a File
-## Copyright 2003-2016 Emmanuel Paradis
+## Copyright 2003-2017 Emmanuel Paradis, 2017 RJ Ewing
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
read.FASTA <- function(file)
{
- if (length(grep("^(ht|f)tp:", file))) {
+ if (length(grep("^(ht|f)tp(s|):", file))) {
url <- file
file <- tempfile()
download.file(url, file)
}
- sz <- file.size(file)
- x <- readBin(file, "raw", sz)
- ## if the file is larger than 2 Gb we assume that it is UNIX-encoded
- ## and skip the search-replacement of carriage returns
+ if (inherits(file, "connection")) {
+ if (!isOpen(file, "rt")) {
+ open(file, "rt")
+ on.exit(close(file))
+ }
+ x <- scan(file, what = character(), sep = "\n", quiet = TRUE)
+ x <- charToRaw(paste(x, collapse = "\n"))
+ sz <- length(x)
+ } else {
+ sz <- file.size(file)
+ x <- readBin(file, "raw", sz)
+ }
+ ## if the file is larger than 1 Gb we assume that it is
+ ## UNIX-encoded and skip the search-replace of carriage returns
if (sz < 1e9) {
icr <- which(x == as.raw(0x0d)) # CR
if (length(icr)) x <- x[-icr]
}
- res <- .Call(rawStreamToDNAbin, x)
+ res <- .Call("rawStreamToDNAbin", x)
+ if (identical(res, 0L)) {
+ warning("failed to read sequences, returns NULL")
+ return(NULL)
+ }
names(res) <- sub("^ +", "", names(res)) # to permit phylosim
class(res) <- "DNAbin"
res
@@ -139,3 +153,22 @@ read.dna <- function(file, format = "interleaved", skip = 0,
}
obj
}
+
+read.fastq <- function(file, offset = -33)
+{
+ Z <- scan(file, "", sep="\n", quiet = TRUE)
+ tmp <- Z[c(TRUE, TRUE, FALSE, FALSE)]
+ sel <- c(TRUE, FALSE)
+ tmp[sel] <- gsub("^@", ">", tmp[sel])
+ fl <- tempfile()
+ cat(tmp, file = fl, sep = "\n")
+ DNA <- read.FASTA(fl)
+
+ ## get the qualities:
+ tmp <- Z[c(FALSE, FALSE, FALSE, TRUE)]
+ QUAL <- lapply(tmp, function(x) as.integer(charToRaw(x)))
+ if (offset) QUAL <- lapply(QUAL, "+", offset)
+ names(QUAL) <- names(DNA)
+ attr(DNA, "QUAL") <- QUAL
+ DNA
+}
diff --git a/R/read.nexus.R b/R/read.nexus.R
index 2c9412c..b34caa2 100644
--- a/R/read.nexus.R
+++ b/R/read.nexus.R
@@ -1,8 +1,8 @@
-## read.nexus.R (2014-10-20)
+## read.nexus.R (2017-01-10)
## Read Tree File in Nexus Format
-## Copyright 2003-2014 Emmanuel Paradis and 2010 Klaus Schliep
+## Copyright 2003-2017 Emmanuel Paradis and 2010 Klaus Schliep
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
@@ -102,7 +102,7 @@ clado.build <- function(tp)
obj
}
-read.nexus <- function(file, tree.names = NULL)
+read.nexus <- function(file, tree.names = NULL, force.multi = FALSE)
{
X <- scan(file = file, what = "", sep = "\n", quiet = TRUE)
## remove all comments
@@ -227,7 +227,7 @@ read.nexus <- function(file, tree.names = NULL)
stop(paste("The tree has apparently singleton node(s): cannot read tree file.\n Reading NEXUS file aborted at tree no.", i, sep = ""))
}
}
- if (Ntree == 1) {
+ if (Ntree == 1 && !force.multi) {
trees <- trees[[1]]
if (translation) {
trees$tip.label <-
diff --git a/inst/doc/MoranI.pdf b/inst/doc/MoranI.pdf
index aa82ee1..4ebc3fa 100644
Binary files a/inst/doc/MoranI.pdf and b/inst/doc/MoranI.pdf differ
diff --git a/man/apetools.Rd b/man/apetools.Rd
new file mode 100644
index 0000000..c231adf
--- /dev/null
+++ b/man/apetools.Rd
@@ -0,0 +1,47 @@
+\name{apetools}
+\alias{apetools}
+\alias{Xplorefiles}
+\alias{editFileExtensions}
+\alias{bydir}
+\title{Tools to Explore Files}
+\description{
+ These functions help to find files on the local disk.
+}
+\usage{
+Xplorefiles(from = "HOME", recursive = TRUE, ignore.case = TRUE)
+editFileExtensions()
+bydir(x)
+}
+\arguments{
+ \item{from}{the directory where to start the file search; by default,
+ the `HOME' directory. Use \code{from = getwd()} to start from the
+ current working directory.}
+ \item{recursive}{whether to search the subdirectories; \code{TRUE} by
+ default.}
+ \item{ignore.case}{whether to ignore the case of the file extensions;
+ \code{TRUE} by default.}
+ \item{x}{a list returned by \code{Xplorefiles}.}
+}
+\details{
+ \code{Xplorefiles} looks for all files with a specified extension in
+ their names. The default is to look for the following file types:
+ CLUSTAL (.aln), FASTA (.fas, .fasta), FASTQ (.fq, .fastq), NEWICK
+ (.nwk, .newick, .tre, .tree), NEXUS (.nex, .nexus), and PHYLIP
+ (.phy). This list can be modified with \code{editFileExtensions}.
+
+ \code{bydir} sorts the list of files by directories.
+}
+\value{
+ \code{Xplorefiles} returns a list. \code{bydir} prints the file
+ listings on the console.
+}
+\author{Emmanuel Paradis}
+\examples{
+\dontrun{
+x <- Xplorefiles()
+x # all data files on your disk
+bydir(x) # sorted by directories
+bydir(x["fasta"]) # only the FASTA files
+Xplorefiles(getwd(), recursive = FALSE) # look only in current dir
+}}
+\keyword{manip}
diff --git a/man/as.bitsplits.Rd b/man/as.bitsplits.Rd
index d9cb25e..a9f58e4 100644
--- a/man/as.bitsplits.Rd
+++ b/man/as.bitsplits.Rd
@@ -52,7 +52,7 @@ as.prop.part(x)
\code{as.prop.part} returns an object of class \code{"prop.part"}.
}
\author{Emmanuel Paradis}
-\seealso{\code{\link{prop.part}}}
+\seealso{\code{\link{prop.part}}, \code{\link{is.compatible}}}
\examples{
tr <- rtree(20)
pp <- prop.part(tr)
diff --git a/man/c.phylo.Rd b/man/c.phylo.Rd
index 283fffa..1a5bfaf 100644
--- a/man/c.phylo.Rd
+++ b/man/c.phylo.Rd
@@ -16,7 +16,7 @@
\arguments{
\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{recursive}{see details.}
\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.}
diff --git a/man/is.compatible.Rd b/man/is.compatible.Rd
index e8f2af7..454ea19 100644
--- a/man/is.compatible.Rd
+++ b/man/is.compatible.Rd
@@ -22,4 +22,5 @@ arecompatible(x, y, n)
\code{TRUE} if the splits are compatible, \code{FALSE} otherwise.
}
\author{Andrei Popescu \email{niteloserpopescu at gmail.com}}
+\seealso{\code{\link{as.bitsplits}}}
\keyword{manip}
diff --git a/man/node.dating.Rd b/man/node.dating.Rd
index 0667b26..d2e10d9 100644
--- a/man/node.dating.Rd
+++ b/man/node.dating.Rd
@@ -1,24 +1,25 @@
\name{node.dating}
\alias{estimate.mu}
\alias{estimate.dates}
+\alias{node.dating}
\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.mu(t, node.dates, p.tol = 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,
+ opt.tol = 1e-8, nsteps = 1000, lik.tol = 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{p.tol}{p-value cutoff for failed regression.}
+ \item{mu}{mutation rate. Either a numeric vector of size 1 for a strict rate or a
+ numeric vector containing the mutation rate along each edge for a relaxed 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.}
@@ -52,6 +53,10 @@ estimate.dates(t, node.dates, mu = estimate.mu(t, node.dates),
priors in the inital step. If all of the dates for nodes are given, then
'estimate.dates' will not run the inital step.
+ 'mu' can be either a vector of size 1 to use a strict molecular clock or a vector
+ with a mutation rate for each edge, in the order of t$edge, for a local molecular
+ clock on each edge.
+
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
diff --git a/man/read.dna.Rd b/man/read.dna.Rd
index 00fee97..d845746 100644
--- a/man/read.dna.Rd
+++ b/man/read.dna.Rd
@@ -1,23 +1,28 @@
\name{read.dna}
\alias{read.dna}
\alias{read.FASTA}
+\alias{read.fastq}
\title{Read DNA Sequences in a File}
\description{
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
- case.
+ in binary format, otherwise (if \code{as.character = TRUE}) in
+ lowercase.
}
\usage{
read.dna(file, format = "interleaved", skip = 0,
nlines = 0, comment.char = "#",
as.character = FALSE, as.matrix = NULL)
read.FASTA(file)
+read.fastq(file, offset = -33)
}
\arguments{
\item{file}{a file name specified by either a variable of mode character,
- or a double-quoted string.}
+ or a double-quoted string. Can also be a \link{connection} (which
+ will be opened for reading if necessary, and if so
+ \code{\link{close}}d (and hence destroyed) at the end of the
+ function call).}
\item{format}{a character string specifying the format of the DNA
sequences. Four choices are possible: \code{"interleaved"},
\code{"sequential"}, \code{"clustal"}, or \code{"fasta"}, or any
@@ -36,6 +41,9 @@ read.FASTA(file)
returns the sequences in a matrix, or stops with an error if they
are of different lengths; (iii) \code{FALSE}: always returns the
sequences in a list.}
+ \item{offset}{the value to be added to the quality scores (the default
+ applies to the Sanger format and should work for most recent FASTQ
+ files).}
}
\details{
\code{read.dna} follows the interleaved and sequential formats defined
@@ -44,8 +52,8 @@ read.FASTA(file)
formats, the first line of the file must contain the dimensions of the
data (the numbers of taxa and the numbers of nucleotides); the
sequences are considered as aligned and thus must be of the same
- lengths for all taxa. For the FASTA format, the conventions defined in
- the URL below (see References) are followed; the sequences are taken as
+ lengths for all taxa. For the FASTA and FASTQ formats, the conventions
+ defined in the references are followed; the sequences are taken as
non-aligned. For all formats, the nucleotides can be arranged in any
way with blanks and line-breaks inside (with the restriction that the
first ten nucleotides must be contiguous for the interleaved and
@@ -56,9 +64,9 @@ read.FASTA(file)
\item{Interleaved:}{the function starts to read the sequences after it
finds one or more spaces (or tabulations). All characters before the
sequences are taken as the taxa names after removing the leading and
- trailing spaces (so spaces in taxa names are allowed). It is assumed
- that the taxa names are not repeated in the subsequent blocks of
- nucleotides.}
+ trailing spaces (so spaces in taxa names are not allowed). It is
+ assumed that the taxa names are not repeated in the subsequent
+ blocks of nucleotides.}
\item{Sequential:}{the same criterion than for the interleaved format
is used to start reading the sequences and the taxa names; the
@@ -66,26 +74,38 @@ read.FASTA(file)
the first line of the file is reached. This is repeated for each taxa.}
\item{Clustal:}{this is the format output by the Clustal programs
- (.aln). It is somehow similar to the interleaved format: the
- differences being that the dimensions of the data are not indicated
- in the file, and the names of the sequences are repeated in each block.}
+ (.aln). It is close to the interleaved format: the differences are
+ that the dimensions of the data are not indicated in the file, and
+ the names of the sequences are repeated in each block.}
\item{FASTA:}{This looks like the sequential format but the taxa names
- (or rather a description of the sequence) are on separate lines
- beginning with a `greater than' character `>' (there may be
- leading spaces before this character). These lines are taken as taxa
- names after removing the `>' and the possible leading and trailing
- spaces. All the data in the file before the first sequence is ignored.}
-}}
+ (or a description of the sequence) are on separate lines beginning
+ with a `greater than' character `>' (there may be leading spaces
+ before this character). These lines are taken as taxa names after
+ removing the `>' and the possible leading and trailing spaces. All
+ the data in the file before the first sequence are ignored.}
+}
+
+The FASTQ format is explained in the references.
+
+Compressed files must be read through connections (see examples).
+\code{read.fastq} can read compressed files directly (see
+examples).
+}
\value{
a matrix or a list (if \code{format = "fasta"}) of DNA sequences
stored in binary format, or of mode character (if \code{as.character =
"TRUE"}).
\code{read.FASTA} always returns a list of class \code{"DNAbin"}.
+
+ \code{read.fastq} returns a list of class \code{"DNAbin"} with an
+ atrribute \code{"QUAL"} (see examples).
}
\references{
- Anonymous. FASTA format. \url{http://en.wikipedia.org/wiki/FASTA_format}
+ Anonymous. FASTA format. \url{https://en.wikipedia.org/wiki/FASTA_format}
+
+ Anonymous. FASTQ format. \url{https://en.wikipedia.org/wiki/FASTQ_format}
Felsenstein, J. (1993) Phylip (Phylogeny Inference Package) version
3.5c. Department of Genetics, University of Washington.
@@ -95,9 +115,9 @@ read.FASTA(file)
\code{\link{read.GenBank}}, \code{\link{write.dna}},
\code{\link{DNAbin}}, \code{\link{dist.dna}}, \code{\link{woodmouse}}
}
-\author{Emmanuel Paradis}
+\author{Emmanuel Paradis and RJ Ewing}
\examples{
-### a small extract from `data(woddmouse)'
+## a small extract from data(woddmouse) in sequential format:
cat("3 40",
"No305 NTTCGAAAAACACACCCACTACTAAAANTTATCAGTCACT",
"No304 ATTCGAAAAACACACCCACTACTAAAAATTATCAACCACT",
@@ -106,7 +126,7 @@ file = "exdna.txt", sep = "\n")
ex.dna <- read.dna("exdna.txt", format = "sequential")
str(ex.dna)
ex.dna
-### the same data in interleaved format...
+## the same data in interleaved format...
cat("3 40",
"No305 NTTCGAAAAA CACACCCACT",
"No304 ATTCGAAAAA CACACCCACT",
@@ -116,7 +136,7 @@ cat("3 40",
" ACTAAAAATT ATCAATCACT",
file = "exdna.txt", sep = "\n")
ex.dna2 <- read.dna("exdna.txt")
-### ... in clustal format...
+## ... in clustal format...
cat("CLUSTAL (ape) multiple sequence alignment", "",
"No305 NTTCGAAAAACACACCCACTACTAAAANTTATCAGTCACT",
"No304 ATTCGAAAAACACACCCACTACTAAAAATTATCAACCACT",
@@ -124,19 +144,56 @@ cat("CLUSTAL (ape) multiple sequence alignment", "",
" ************************** ****** ****",
file = "exdna.txt", sep = "\n")
ex.dna3 <- read.dna("exdna.txt", format = "clustal")
-### ... and in FASTA format
+## ... and in FASTA format
cat(">No305",
"NTTCGAAAAACACACCCACTACTAAAANTTATCAGTCACT",
">No304",
"ATTCGAAAAACACACCCACTACTAAAAATTATCAACCACT",
">No306",
"ATTCGAAAAACACACCCACTACTAAAAATTATCAATCACT",
-file = "exdna.txt", sep = "\n")
-ex.dna4 <- read.dna("exdna.txt", format = "fasta")
-### They are the same!
+file = "exdna.fas", sep = "\n")
+ex.dna4 <- read.dna("exdna.fas", format = "fasta")
+## They are the same:
identical(ex.dna, ex.dna2)
identical(ex.dna, ex.dna3)
identical(ex.dna, ex.dna4)
+
+## How to read compressed files:
+
+## create the ZIP file:
+zip("exdna.fas.zip", "exdna.fas")
+## create the GZ file with a connection:
+con <- gzfile("exdna.fas.gz", "wt")
+cat(">No305", "NTTCGAAAAACACACCCACTACTAAAANTTATCAGTCACT",
+ ">No304", "ATTCGAAAAACACACCCACTACTAAAAATTATCAACCACT",
+ ">No306", "ATTCGAAAAACACACCCACTACTAAAAATTATCAATCACT",
+ file = con, sep = "\n")
+close(con)
+
+ex.dna5 <- read.dna(unz("exdna.fas.zip", "exdna.fas"), "fasta")
+ex.dna6 <- read.dna(gzfile("exdna.fas.gz"), "fasta")
+identical(ex.dna5, ex.dna4)
+identical(ex.dna6, ex.dna4)
+
unlink("exdna.txt") # clean-up
-}
+unlink("exdna.fas") # clean-up
+unlink("exdna.fas.zip")
+unlink("exdna.fas.gz")
+
+## read a FASTQ file from 1000 Genomes:
+\dontrun{
+a <- "ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG00096/sequence_read/"
+b <- "SRR062641.filt.fastq.gz"
+URL <- paste0(a, b)
+download.file(URL, b)
+X <- read.fastq(b)
+X # 109,811 sequences
+## get the qualities of the first sequence:
+(qual1 <- attr(X, "QUAL")[[1]])
+## the corresponding probabilities:
+10^(-qual1/10)
+## get the mean quality for each sequence:
+mean.qual <- sapply(attr(X, "Q"), mean)
+## can do the same for var, sd, ...
+}}
\keyword{IO}
diff --git a/man/read.nexus.Rd b/man/read.nexus.Rd
index 5594e0e..5aa9e3c 100644
--- a/man/read.nexus.Rd
+++ b/man/read.nexus.Rd
@@ -2,7 +2,7 @@
\alias{read.nexus}
\title{Read Tree File in Nexus Format}
\usage{
-read.nexus(file, tree.names = NULL)
+read.nexus(file, tree.names = NULL, force.multi = FALSE)
}
\arguments{
\item{file}{a file name specified by either a variable of mode character,
@@ -10,6 +10,9 @@ read.nexus(file, tree.names = NULL)
\item{tree.names}{if there are several trees to be read, a vector of
mode character giving names to the individual trees (by default,
this uses the labels in the NEXUS file if these are present).}
+ \item{force.multi}{a logical value; if \code{TRUE}, an object of class
+ \code{"multiPhylo"} is always returned even if the file contains a
+ single tree (see details).}
}
\description{
This function reads one or several trees in a NEXUS file.
@@ -30,6 +33,9 @@ read.nexus(file, tree.names = NULL)
consistent with \code{\link{write.nexus}} which translates only the
tip labels.
+ Using \code{force.multi = TRUE} when the file contains a single tree
+ makes possible to keep the tree name (as names of the list).
+
`read.nexus' tries to represent correctly trees with a badly
represented root edge (i.e. with an extra pair of parentheses). For
instance, the tree "((A:1,B:1):10);" will be read like "(A:1,B:1):10;"
@@ -39,7 +45,7 @@ read.nexus(file, tree.names = NULL)
message is issued.
}
\value{
- an object of class \code{"phylo"} or of class \code{"multiPhylo"}.
+ an object of class \code{"phylo"} or \code{"multiPhylo"}.
}
\references{
Maddison, D. R., Swofford, D. L. and Maddison, W. P. (1997) NEXUS: an
diff --git a/src/read_dna.c b/src/read_dna.c
index c1ef4f0..4427503 100644
--- a/src/read_dna.c
+++ b/src/read_dna.c
@@ -1,6 +1,6 @@
-/* read_dna.c 2016-06-21 */
+/* read_dna.c 2017-01-04 */
-/* Copyright 2013-2016 Emmanuel Paradis */
+/* Copyright 2013-2017 Emmanuel Paradis */
/* This file is part of the R-package `ape'. */
/* See the file ../COPYING for licensing issues. */
@@ -112,6 +112,13 @@ SEXP rawStreamToDNAbin(SEXP x)
}
}
+ if (n == 0) {
+ PROTECT(obj = allocVector(INTSXP, 1));
+ INTEGER(obj)[0] = 0;
+ UNPROTECT(2);
+ return obj;
+ }
+
PROTECT(obj = allocVector(VECSXP, n));
PROTECT(nms = allocVector(STRSXP, n));
--
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