API tools faq
paste
Advertisement
Guest User

my.plot.phylo.r

a guest
Aug 23rd, 2017
345
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 19.89 KB | None | 0 0
raw download clone embed print report
  1. my.plot.phylo <- function (x, type = "phylogram", use.edge.length = TRUE, node.pos = NULL,
  2. show.tip.label = TRUE, show.node.label = FALSE, edge.color = "black",
  3. edge.width = 1, edge.lty = 1, font = 3, cex = par("cex"),
  4. adj = NULL, srt = 0, no.margin = FALSE, root.edge = FALSE,
  5. label.offset = 0, underscore = FALSE, x.lim = NULL, y.lim = NULL,
  6. direction = "rightwards", lab4ut = NULL, tip.color = "black",
  7. plot = TRUE, rotate.tree = 0, open.angle = 0, node.depth = 1,
  8. align.tip.label = FALSE, ...)
  9. {
  10. Ntip <- length(x$tip.label)
  11. if (Ntip < 2) {
  12. warning("found less than 2 tips in the tree")
  13. return(NULL)
  14. }
  15. if (any(tabulate(x$edge[, 1]) == 1))
  16. stop("there are single (non-splitting) nodes in your tree; you may need to use collapse.singles()")
  17. .nodeHeight <- function(Ntip, Nnode, edge, Nedge, yy) .C(node_height,
  18. as.integer(Ntip), as.integer(Nnode), as.integer(edge[,
  19. 1]), as.integer(edge[, 2]), as.integer(Nedge), as.double(yy))[[6]]
  20. .nodeDepth <- function(Ntip, Nnode, edge, Nedge, node.depth) .C(node_depth,
  21. as.integer(Ntip), as.integer(Nnode), as.integer(edge[,
  22. 1]), as.integer(edge[, 2]), as.integer(Nedge), double(Ntip +
  23. Nnode), as.integer(node.depth))[[6]]
  24. .nodeDepthEdgelength <- function(Ntip, Nnode, edge, Nedge,
  25. edge.length) .C(node_depth_edgelength, as.integer(Ntip),
  26. as.integer(Nnode), as.integer(edge[, 1]), as.integer(edge[,
  27. 2]), as.integer(Nedge), as.double(edge.length), double(Ntip +
  28. Nnode))[[7]]
  29. Nedge <- dim(x$edge)[1]
  30. Nnode <- x$Nnode
  31. if (any(x$edge < 1) || any(x$edge > Ntip + Nnode))
  32. stop("tree badly conformed; cannot plot. Check the edge matrix.")
  33. ROOT <- Ntip + 1
  34. type <- match.arg(type, c("phylogram", "cladogram", "fan",
  35. "unrooted", "radial"))
  36. direction <- match.arg(direction, c("rightwards", "leftwards",
  37. "upwards", "downwards"))
  38. if (is.null(x$edge.length)) {
  39. use.edge.length <- FALSE
  40. }
  41. else {
  42. if (use.edge.length && type != "radial") {
  43. tmp <- sum(is.na(x$edge.length))
  44. if (tmp) {
  45. warning(paste(tmp, "branch length(s) NA(s): branch lengths ignored in the plot"))
  46. use.edge.length <- FALSE
  47. }
  48. }
  49. }
  50. if (is.numeric(align.tip.label)) {
  51. align.tip.label.lty <- align.tip.label
  52. align.tip.label <- TRUE
  53. }
  54. else {
  55. if (align.tip.label)
  56. align.tip.label.lty <- 1
  57. }
  58. if (align.tip.label) {
  59. if (type %in% c("unrooted", "radial") || !use.edge.length ||
  60. is.ultrametric(x))
  61. align.tip.label <- FALSE
  62. }
  63. if (type %in% c("unrooted", "radial") || !use.edge.length ||
  64. is.null(x$root.edge) || !x$root.edge)
  65. root.edge <- FALSE
  66. phyloORclado <- type %in% c("phylogram", "cladogram")
  67. horizontal <- direction %in% c("rightwards", "leftwards")
  68. xe <- x$edge
  69. if (phyloORclado) {
  70. phyOrder <- attr(x, "order")
  71. if (is.null(phyOrder) || phyOrder != "cladewise") {
  72. x <- reorder(x)
  73. if (!identical(x$edge, xe)) {
  74. ereorder <- match(x$edge[, 2], xe[, 2])
  75. if (length(edge.color) > 1) {
  76. edge.color <- rep(edge.color, length.out = Nedge)
  77. edge.color <- edge.color[ereorder]
  78. }
  79. if (length(edge.width) > 1) {
  80. edge.width <- rep(edge.width, length.out = Nedge)
  81. edge.width <- edge.width[ereorder]
  82. }
  83. if (length(edge.lty) > 1) {
  84. edge.lty <- rep(edge.lty, length.out = Nedge)
  85. edge.lty <- edge.lty[ereorder]
  86. }
  87. }
  88. }
  89. yy <- numeric(Ntip + Nnode)
  90. TIPS <- x$edge[x$edge[, 2] <= Ntip, 2]
  91. yy[TIPS] <- 1:Ntip
  92. }
  93. z <- reorder(x, order = "postorder")
  94. if (phyloORclado) {
  95. if (is.null(node.pos))
  96. node.pos <- if (type == "cladogram" && !use.edge.length)
  97. 2
  98. else 1
  99. if (node.pos == 1)
  100. yy <- .nodeHeight(Ntip, Nnode, z$edge, Nedge, yy)
  101. else {
  102. ans <- .C(node_height_clado, as.integer(Ntip), as.integer(Nnode),
  103. as.integer(z$edge[, 1]), as.integer(z$edge[,
  104. 2]), as.integer(Nedge), double(Ntip + Nnode),
  105. as.double(yy))
  106. xx <- ans[[6]] - 1
  107. yy <- ans[[7]]
  108. }
  109. if (!use.edge.length) {
  110. if (node.pos != 2)
  111. xx <- .nodeDepth(Ntip, Nnode, z$edge, Nedge,
  112. node.depth) - 1
  113. xx <- max(xx) - xx
  114. }
  115. else {
  116. xx <- .nodeDepthEdgelength(Ntip, Nnode, z$edge, Nedge,
  117. z$edge.length)
  118. }
  119. }
  120. else {
  121. twopi <- 2 * pi
  122. rotate.tree <- twopi * rotate.tree/360
  123. if (type != "unrooted") {
  124. TIPS <- x$edge[which(x$edge[, 2] <= Ntip), 2]
  125. xx <- seq(0, twopi * (1 - 1/Ntip) - twopi * open.angle/360,
  126. length.out = Ntip)
  127. theta <- double(Ntip)
  128. theta[TIPS] <- xx
  129. theta <- c(theta, numeric(Nnode))
  130. }
  131. switch(type, fan = {
  132. theta <- .nodeHeight(Ntip, Nnode, z$edge, Nedge,
  133. theta)
  134. if (use.edge.length) {
  135. r <- .nodeDepthEdgelength(Ntip, Nnode, z$edge,
  136. Nedge, z$edge.length)
  137. } else {
  138. r <- .nodeDepth(Ntip, Nnode, z$edge, Nedge, node.depth)
  139. r <- 1/r
  140. }
  141. theta <- theta + rotate.tree
  142. if (root.edge) r <- r + x$root.edge
  143. xx <- r * cos(theta)
  144. yy <- r * sin(theta)
  145. }, unrooted = {
  146. nb.sp <- .nodeDepth(Ntip, Nnode, z$edge, Nedge, node.depth)
  147. XY <- if (use.edge.length) unrooted.xy(Ntip, Nnode,
  148. z$edge, z$edge.length, nb.sp, rotate.tree) else unrooted.xy(Ntip,
  149. Nnode, z$edge, rep(1, Nedge), nb.sp, rotate.tree)
  150. xx <- XY$M[, 1] - min(XY$M[, 1])
  151. yy <- XY$M[, 2] - min(XY$M[, 2])
  152. }, radial = {
  153. r <- .nodeDepth(Ntip, Nnode, z$edge, Nedge, node.depth)
  154. r[r == 1] <- 0
  155. r <- 1 - r/Ntip
  156. theta <- .nodeHeight(Ntip, Nnode, z$edge, Nedge,
  157. theta) + rotate.tree
  158. xx <- r * cos(theta)
  159. yy <- r * sin(theta)
  160. })
  161. }
  162. if (phyloORclado) {
  163. if (!horizontal) {
  164. tmp <- yy
  165. yy <- xx
  166. xx <- tmp - min(tmp) + 1
  167. }
  168. if (root.edge) {
  169. if (direction == "rightwards")
  170. xx <- xx + x$root.edge
  171. if (direction == "upwards")
  172. yy <- yy + x$root.edge
  173. }
  174. }
  175. if (no.margin)
  176. par(mai = rep(0, 4))
  177. if (show.tip.label)
  178. nchar.tip.label <- nchar(x$tip.label)
  179. max.yy <- max(yy)
  180. if (is.null(x.lim)) {
  181. if (phyloORclado) {
  182. if (horizontal) {
  183. x.lim <- c(0, NA)
  184. pin1 <- par("pin")[1]
  185. strWi <- strwidth(x$tip.label, "inches", cex = cex)
  186. xx.tips <- xx[1:Ntip] * 1.04
  187. alp <- try(uniroot(function(a) max(a * xx.tips +
  188. strWi) - pin1, c(0, 1e+06))$root, silent = TRUE)
  189. if (is.character(alp)) {
  190. tmp <- max(xx.tips)
  191. if (show.tip.label)
  192. tmp <- tmp * 1.5
  193. }
  194. else {
  195. tmp <- if (show.tip.label)
  196. max(xx.tips + strWi/alp)
  197. else max(xx.tips)
  198. }
  199. if (show.tip.label)
  200. tmp <- tmp + label.offset
  201. x.lim[2] <- tmp
  202. }
  203. else x.lim <- c(1, Ntip)
  204. }
  205. else switch(type, fan = {
  206. if (show.tip.label) {
  207. offset <- max(nchar.tip.label * 0.018 * max.yy *
  208. cex)
  209. x.lim <- range(xx) + c(-offset, offset)
  210. } else x.lim <- range(xx)
  211. }, unrooted = {
  212. if (show.tip.label) {
  213. offset <- max(nchar.tip.label * 0.018 * max.yy *
  214. cex)
  215. x.lim <- c(0 - offset, max(xx) + offset)
  216. } else x.lim <- c(0, max(xx))
  217. }, radial = {
  218. if (show.tip.label) {
  219. offset <- max(nchar.tip.label * 0.03 * cex)
  220. x.lim <- c(-1 - offset, 1 + offset)
  221. } else x.lim <- c(-1, 1)
  222. })
  223. }
  224. else if (length(x.lim) == 1) {
  225. x.lim <- c(0, x.lim)
  226. if (phyloORclado && !horizontal)
  227. x.lim[1] <- 1
  228. if (type %in% c("fan", "unrooted") && show.tip.label)
  229. x.lim[1] <- -max(nchar.tip.label * 0.018 * max.yy *
  230. cex)
  231. if (type == "radial")
  232. x.lim[1] <- if (show.tip.label)
  233. -1 - max(nchar.tip.label * 0.03 * cex)
  234. else -1
  235. }
  236. if (phyloORclado && direction == "leftwards")
  237. xx <- x.lim[2] - xx
  238. if (is.null(y.lim)) {
  239. if (phyloORclado) {
  240. if (horizontal)
  241. y.lim <- c(1, Ntip)
  242. else {
  243. y.lim <- c(0, NA)
  244. pin2 <- par("pin")[2]
  245. strWi <- strwidth(x$tip.label, "inches", cex = cex)
  246. yy.tips <- yy[1:Ntip] * 1.04
  247. alp <- try(uniroot(function(a) max(a * yy.tips +
  248. strWi) - pin2, c(0, 1e+06))$root, silent = TRUE)
  249. if (is.character(alp)) {
  250. tmp <- max(yy.tips)
  251. if (show.tip.label)
  252. tmp <- tmp * 1.5
  253. }
  254. else {
  255. tmp <- if (show.tip.label)
  256. max(yy.tips + strWi/alp)
  257. else max(yy.tips)
  258. }
  259. if (show.tip.label)
  260. tmp <- tmp + label.offset
  261. y.lim[2] <- tmp
  262. }
  263. }
  264. else switch(type, fan = {
  265. if (show.tip.label) {
  266. offset <- max(nchar.tip.label * 0.018 * max.yy *
  267. cex)
  268. y.lim <- c(min(yy) - offset, max.yy + offset)
  269. } else y.lim <- c(min(yy), max.yy)
  270. }, unrooted = {
  271. if (show.tip.label) {
  272. offset <- max(nchar.tip.label * 0.018 * max.yy *
  273. cex)
  274. y.lim <- c(0 - offset, max.yy + offset)
  275. } else y.lim <- c(0, max.yy)
  276. }, radial = {
  277. if (show.tip.label) {
  278. offset <- max(nchar.tip.label * 0.03 * cex)
  279. y.lim <- c(-1 - offset, 1 + offset)
  280. } else y.lim <- c(-1, 1)
  281. })
  282. }
  283. else if (length(y.lim) == 1) {
  284. y.lim <- c(0, y.lim)
  285. if (phyloORclado && horizontal)
  286. y.lim[1] <- 1
  287. if (type %in% c("fan", "unrooted") && show.tip.label)
  288. y.lim[1] <- -max(nchar.tip.label * 0.018 * max.yy *
  289. cex)
  290. if (type == "radial")
  291. y.lim[1] <- if (show.tip.label)
  292. -1 - max(nchar.tip.label * 0.018 * max.yy * cex)
  293. else -1
  294. }
  295. if (phyloORclado && direction == "downwards")
  296. yy <- y.lim[2] - yy
  297. if (phyloORclado && root.edge) {
  298. if (direction == "leftwards")
  299. x.lim[2] <- x.lim[2] + x$root.edge
  300. if (direction == "downwards")
  301. y.lim[2] <- y.lim[2] + x$root.edge
  302. }
  303. asp <- if (type %in% c("fan", "radial", "unrooted"))
  304. 1
  305. else NA
  306. plot.default(0, type = "n", xlim = x.lim, ylim = y.lim, xlab = "",
  307. ylab = "", axes = FALSE, asp = asp, ...)
  308. if (plot) {
  309. if (is.null(adj))
  310. adj <- if (phyloORclado && direction == "leftwards")
  311. 1
  312. else 0
  313. if (phyloORclado && show.tip.label) {
  314. MAXSTRING <- max(strwidth(x$tip.label, cex = cex))
  315. loy <- 0
  316. if (direction == "rightwards") {
  317. lox <- label.offset + MAXSTRING * 1.05 * adj
  318. }
  319. if (direction == "leftwards") {
  320. lox <- -label.offset - MAXSTRING * 1.05 * (1 -
  321. adj)
  322. }
  323. if (!horizontal) {
  324. psr <- par("usr")
  325. MAXSTRING <- MAXSTRING * 1.09 * (psr[4] - psr[3])/(psr[2] -
  326. psr[1])
  327. loy <- label.offset + MAXSTRING * 1.05 * adj
  328. lox <- 0
  329. srt <- 90 + srt
  330. if (direction == "downwards") {
  331. loy <- -loy
  332. srt <- 180 + srt
  333. }
  334. }
  335. }
  336. if (type == "phylogram") {
  337. phylogram.plot(x$edge, Ntip, Nnode, xx, yy, horizontal,
  338. edge.color, edge.width, edge.lty)
  339. }
  340. else {
  341. if (type == "fan") {
  342. ereorder <- match(z$edge[, 2], x$edge[, 2])
  343. if (length(edge.color) > 1) {
  344. edge.color <- rep(edge.color, length.out = Nedge)
  345. edge.color <- edge.color[ereorder]
  346. }
  347. if (length(edge.width) > 1) {
  348. edge.width <- rep(edge.width, length.out = Nedge)
  349. edge.width <- edge.width[ereorder]
  350. }
  351. if (length(edge.lty) > 1) {
  352. edge.lty <- rep(edge.lty, length.out = Nedge)
  353. edge.lty <- edge.lty[ereorder]
  354. }
  355. circular.plot(z$edge, Ntip, Nnode, xx, yy, theta,
  356. r, edge.color, edge.width, edge.lty)
  357. }
  358. else cladogram.plot(x$edge, xx, yy, edge.color, edge.width,
  359. edge.lty)
  360. }
  361. if (root.edge) {
  362. rootcol <- if (length(edge.color) == 1)
  363. edge.color
  364. else "black"
  365. rootw <- if (length(edge.width) == 1)
  366. edge.width
  367. else 1
  368. rootlty <- if (length(edge.lty) == 1)
  369. edge.lty
  370. else 1
  371. if (type == "fan") {
  372. tmp <- polar2rect(x$root.edge, theta[ROOT])
  373. segments(0, 0, tmp$x, tmp$y, col = rootcol, lwd = rootw,
  374. lty = rootlty)
  375. }
  376. else {
  377. switch(direction, rightwards = segments(0, yy[ROOT],
  378. x$root.edge, yy[ROOT], col = rootcol, lwd = rootw,
  379. lty = rootlty), leftwards = segments(xx[ROOT],
  380. yy[ROOT], xx[ROOT] + x$root.edge, yy[ROOT],
  381. col = rootcol, lwd = rootw, lty = rootlty),
  382. upwards = segments(xx[ROOT], 0, xx[ROOT], x$root.edge,
  383. col = rootcol, lwd = rootw, lty = rootlty),
  384. downwards = segments(xx[ROOT], yy[ROOT], xx[ROOT],
  385. yy[ROOT] + x$root.edge, col = rootcol, lwd = rootw,
  386. lty = rootlty))
  387. }
  388. }
  389. if (show.tip.label) {
  390. if (is.expression(x$tip.label))
  391. underscore <- TRUE
  392. if (!underscore)
  393. x$tip.label <- gsub("_", " ", x$tip.label)
  394. if (phyloORclado) {
  395. if (align.tip.label) {
  396. xx.tmp <- switch(direction, rightwards = max(xx[1:Ntip]),
  397. leftwards = min(xx[1:Ntip]), upwards = xx[1:Ntip],
  398. downwards = xx[1:Ntip])
  399. yy.tmp <- switch(direction, rightwards = yy[1:Ntip],
  400. leftwards = yy[1:Ntip], upwards = max(yy[1:Ntip]),
  401. downwards = min(yy[1:Ntip]))
  402. segments(xx[1:Ntip], yy[1:Ntip], xx.tmp, yy.tmp,
  403. lty = align.tip.label.lty)
  404. }
  405. else {
  406. xx.tmp <- xx[1:Ntip]
  407. yy.tmp <- yy[1:Ntip]
  408. }
  409. text(xx.tmp + lox, yy.tmp + loy, x$tip.label,
  410. adj = adj, font = font, srt = srt, cex = cex,
  411. col = tip.color)
  412. }
  413. else {
  414. angle <- if (type == "unrooted")
  415. XY$axe
  416. else atan2(yy[1:Ntip], xx[1:Ntip])
  417. lab4ut <- if (is.null(lab4ut)) {
  418. if (type == "unrooted")
  419. "horizontal"
  420. else "axial"
  421. }
  422. else match.arg(lab4ut, c("horizontal", "axial"))
  423. xx.tips <- xx[1:Ntip]
  424. yy.tips <- yy[1:Ntip]
  425. if (label.offset) {
  426. xx.tips <- xx.tips + label.offset * cos(angle)
  427. yy.tips <- yy.tips + label.offset * sin(angle)
  428. }
  429. if (lab4ut == "horizontal") {
  430. y.adj <- x.adj <- numeric(Ntip)
  431. sel <- abs(angle) > 0.75 * pi
  432. x.adj[sel] <- -strwidth(x$tip.label)[sel] *
  433. 1.05
  434. sel <- abs(angle) > pi/4 & abs(angle) < 0.75 *
  435. pi
  436. x.adj[sel] <- -strwidth(x$tip.label)[sel] *
  437. (2 * abs(angle)[sel]/pi - 0.5)
  438. sel <- angle > pi/4 & angle < 0.75 * pi
  439. y.adj[sel] <- strheight(x$tip.label)[sel]/2
  440. sel <- angle < -pi/4 & angle > -0.75 * pi
  441. y.adj[sel] <- -strheight(x$tip.label)[sel] *
  442. 0.75
  443. text(xx.tips + x.adj * cex, yy.tips + y.adj *
  444. cex, x$tip.label, adj = c(adj, 0), font = font,
  445. srt = srt, cex = cex, col = tip.color)
  446. }
  447. else {
  448. if (align.tip.label) {
  449. POL <- rect2polar(xx.tips, yy.tips)
  450. POL$r[] <- max(POL$r)
  451. REC <- polar2rect(POL$r, POL$angle)
  452. xx.tips <- REC$x
  453. yy.tips <- REC$y
  454. segments(xx[1:Ntip], yy[1:Ntip], xx.tips,
  455. yy.tips, lty = align.tip.label.lty)
  456. }
  457. if (type == "unrooted") {
  458. adj <- abs(angle) > pi/2
  459. angle <- angle * 180/pi
  460. angle[adj] <- angle[adj] - 180
  461. adj <- as.numeric(adj)
  462. }
  463. else {
  464. s <- xx.tips < 0
  465. angle <- angle * 180/pi
  466. angle[s] <- angle[s] + 180
  467. adj <- as.numeric(s)
  468. }
  469. font <- rep(font, length.out = Ntip)
  470. tip.color <- rep(tip.color, length.out = Ntip)
  471. cex <- rep(cex, length.out = Ntip)
  472. for (i in 1:Ntip) text(xx.tips[i], yy.tips[i],
  473. x$tip.label[i], font = font[i], cex = cex[i],
  474. srt = angle[i], adj = adj[i], col = tip.color[i])
  475. }
  476. }
  477. }
  478. if (show.node.label)
  479. text(xx[ROOT:length(xx)] + label.offset, yy[ROOT:length(yy)],
  480. x$node.label, adj = adj, font = font, srt = srt,
  481. cex = cex)
  482. }
  483. L <- list(type = type, use.edge.length = use.edge.length,
  484. node.pos = node.pos, node.depth = node.depth, show.tip.label = show.tip.label,
  485. show.node.label = show.node.label, font = font, cex = cex,
  486. adj = adj, srt = srt, no.margin = no.margin, label.offset = label.offset,
  487. x.lim = x.lim, y.lim = y.lim, direction = direction,
  488. tip.color = tip.color, Ntip = Ntip, Nnode = Nnode, root.time = x$root.time,
  489. align.tip.label = align.tip.label)
  490. assign("last_plot.phylo", c(L, list(edge = xe, xx = xx, yy = yy)),
  491. envir = .PlotPhyloEnv)
  492. invisible(L)
  493. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!