用R在系统发育边缘上绘制 [英] Draw on a phylogeny edge with R

问题描述

我想在带R的系统发育学边缘的任意位置绘制一个符号(十字).

I want to draw a symbol (cross) anywhere along an edge of a phylogeny with R.

让我们带上那棵树:

tree.mod

(((((hg19:0.0295977,macFas5:0.0351997)hg19-macFas5:0.0796862,otoGar3:0.153768)hg19-otoGar3:0.0189364,(((speTri2:0.136745,(mm10:0.0836004,rn5:0.0894755)mm10-rn5:0.221692)speTri2-mm10:0.009583,cavPor3:0.230585)speTri2-cavPor3:0.0279417,oryCun2:0.212423)speTri2-oryCun2:0.0141334)hg19-speTri2:0.0224427,((susScr3:0.127502,(camFer1:0.111428((HLbalAcu1:0.0196864,(HLphyCat1:0.0216245,(HLlipVex1:0.0178214,(turTru2:0.00660779,orcOrc1:0.00557822)turTru2-orcOrc1:0.0113362)HLlipVex1-turTru2:0.00784223)HLphyCat1-HLlipVex1:0.00385294)HLbalAcu1-HLphyCat1:0.0417511(bosTau7:0.0320796,oviAri3:0.0370855)bosTau7-oviAri3:0.0963575)HLbalAcu1-bosTau7:0.0212413)camFer1-HLbalAcu1:0.00459573)susScr3-camFer1:0.0425276,((equCab2:0.107972,(felCat5:0.0893591,(canFam3:0.0908515,((ailMel1:0.0208923,HLursMar1:0.0194161)ailMel1-HLursMar1:0.041512,(odoRosDiv1:0.0259201,lepWed1:0.0238902)odoRosDiv1-lepWed1:0.0288135)ailMel1-odoRosDiv1:0.0213261)canFam3-ailMel1:0.0207311)felCat5-canFam3:0.0501969)equCab2-felCat5:0.00485202,myoLuc2:0.169981)equCab2-myoLuc2:0.00361502)susScr3-equCab2:0.0302374)hg19-susScr3:0.0217021,(((loxAfr3:0.0821245,triMan1:0.0662932)loxAfr3-triMan1:0.0296365,echTel2:0.245195)loxAfr3-echTel2:0.0492624,dasNov3:0.159632)loxAfr3-dasNov3:0.00825218)hg19-loxAfr3;

并在hg19提示处用符号进行绘制:

and plot it with a symbol at the hg19 tip:

R
library(ape)
tree = read.tree("tree.mod")
plot(tree)
points( 0.172365, 1, col="red", pch=20)

容易理解的技巧:'x'值为完整的分支长度，'y'为1,2,3 ...

that is easy for the tips: the 'x' value is the full branch length, and 'y' is 1,2,3...

但我看不到如何对内部节点(例如loxAfr3-triMan1)执行此操作.我有'x'但找不到'y'...

but I don't see how to do it for inner nodes like for instance loxAfr3-triMan1. I have the 'x' but cannot find the 'y'...

欢迎任何输入！

推荐答案

确定.我觉得我肯定想念一种更简单的方法，但是这些代码中的许多似乎都藏在了由程序包的plotting函数调用的C函数中.不过，R中的一些功能应该能够提取特定节点的x和y位置.

OK. I feel like I must be missing an easier way, but much of those code seems to be tucked away inside C functions called by the package's plotting function. Nevertheless, here are some functions in R that should be able to extract the x and y positions for particular nodes.

getphylo_x <- function(tree, node) {
    if(is.character(node)) {
        node <- which(c(tree$tip.label, tree$node.label)==node)
    }
    pi <- tree$edge[tree$edge[,2]==node, 1]
    if (length(pi)) {
        ei<-which(tree$edge[,1]==pi & tree$edge[,2]==node)
        tree$edge.length[ei] + Recall(tree, pi)
    } else {
        if(!is.null(tree$root.edge)) {
            tree$root.edge
        } else {
            0
        }
    }
}

getphylo_y <- function(tree, node) {
    if(is.character(node)) {
        node <- which(c(tree$tip.label, tree$node.label)==node)
    }
    ci <- tree$edge[tree$edge[,1]==node, 2]
    if (length(ci)==2) {
        mean(c(Recall(tree, ci[1]), Recall(tree, ci[2])))
    } else if (length(ci)==0) {
        Ntip <- length(tree$tip.label)
        which(tree$edge[tree$edge[, 2] <= Ntip, 2] == node)
    } else {
        stop(paste("error", length(ci)))
    }
}

要使用它们，只需传递树和节点名称即可.例如

To use them, you just pass in your tree and your node name. For example

node <- "loxAfr3-triMan1"
x <- getphylo_x(tree, node)
y <- getphylo_y(tree, node)


plot(tree)
points(x,y,pch=18, col="red", cex=2)

产生

这篇关于用R在系统发育边缘上绘制的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！