Nothing
R/phyloTreeAnno.R
In MesKit: A tool kit for dissecting cancer evolution from multi-region derived tumor biopsies via somatic alterations
Defines functions
labelBranch
getSigColors
calChildNodeNum
getNodeAngle
calMainTrunk
setPhyloTree
getTreeData
plotTree
plotTree <- function(phyloTree,
treeData = NULL,
branchCol = "mutType",
show.bootstrap = TRUE,
signaturesRef = "cosmic_v2",
min.mut.count = 15,
min.ratio = 1/20,
uncommon.col = "red",
compare.tree.name = "",
use.tumorSampleLabel = FALSE,
show.scale.bar = FALSE,
scale.bar.x = NULL,
scale.bar.y = NULL){
if(min.ratio <= 0|min.ratio > 1){
stop("min.ratio should be within (0,1]")
}
## adjust length of branches by min.ratio
tree <- getTree(phyloTree)
min.len <- max(tree$edge.length)*min.ratio
tree$edge.length[tree$edge.length < min.len] <- min.len
phyloTree <- new('phyloTree',
patientID = getPhyloTreePatient(phyloTree),
tree = tree,
binary.matrix = getBinaryMatrix(phyloTree),
ccf.matrix = getCCFMatrix(phyloTree),
mut.branches = getMutBranches(phyloTree),
branch.type = getBranchType(phyloTree),
ref.build = getPhyloTreeRef(phyloTree),
bootstrap.value = getBootstrapValue(phyloTree),
method = getTreeMethod(phyloTree),
tsb.label = getPhyloTreeTsbLabel(phyloTree))
patient <- getPhyloTreePatient(phyloTree)
if(is.null(treeData)){
treeData <- getTreeData(phyloTree = phyloTree,
branchCol = branchCol,
signaturesRef = signaturesRef,
min.mut.count = min.mut.count)
if(identical(treeData, NA)){
return(NA)
}
compare <- FALSE
}else{
compare <- TRUE
}
# set.seed(1234)
## get bootstrap value
boot_value <- getBootstrapValue(phyloTree)
rootLabel <- "NORMAL"
## plot phylotree
samplePointsSize <- 3
comparePointsSize <- 3.5
sampleTextSize <- 3.5
nodePointsSize <- 1.7
segmentSize <- 1.5
nodeStrokeSize <- 0.25
sampleStrokeSize <- 1.5
bootLabelSize <- 3
bootPaddingSize <- 0.35
bootLabelPaddingSize <- 0.15
legend.text.size <- 10
legend.title.size <- legend.text.size + 0.5
samplesLength <- nrow(treeData[treeData$sample != "internal node",])
if(samplesLength > 7){
samplePointsSize <- 1.5
comparePointsSize <- 2
sampleTextSize <- 3
segmentSize <- 0.8
nodePointsSize <- 0.8
nodeStrokeSize <- 0.15
sampleStrokeSize <- 0.8
bootLabelSize <- 2.5
bootPaddingSize <- 0.1
bootLabelPaddingSize <- 0.1
legend.text.size <- 9
legend.title.size <- legend.text.size + 0.5
}
rootNode <- treeData[treeData$sample == rootLabel,]$node
if(length(boot_value) == 1){
bootsData <- data.frame(x2 = 0, y2 = 0, node = rootNode, end_num = rootNode, boots = boot_value)
}else{
sub <- data.table::data.table(x2 = 0, y2 = 0,node = rootNode, end_num = rootNode)
bootsData <- rbind(treeData[treeData$sample == 'internal node',] %>%
dplyr::select("x2","y2","node","end_num"),sub)
boots <- c()
LN <- min(bootsData$node)-1
boots_list <- lapply(seq_len(nrow(bootsData)), function(i){
if(i == nrow(bootsData)){
boots <- append(boots, boot_value[rootNode - LN])
}else{
boots <- append(boots, boot_value[bootsData$end_num[i] - LN])
}
return(boots)
}) %>% unlist()
bootsData <- cbind(bootsData, boots = boots_list)
}
## initialize
x2 <- NULL
y2 <- NULL
x1 <- NULL
y1 <- NULL
## get the max value of X axis
x_max <- max(abs(treeData$x2))
# print(treeData)
p <- ggplot(data = treeData) +
## balance the space on the left and right sides
geom_segment(aes(x = 0, y = 0, xend = x_max, yend = 0),color = "white",size = 0.01)+
geom_segment(aes(x = 0, y = 0, xend = -x_max, yend = 0),color = "white",size = 0.01)
textAdjust <- mean(as.numeric(treeData$distance))
if(use.tumorSampleLabel){
tsb.label <- getPhyloTreeTsbLabel(phyloTree)
if(nrow(tsb.label) == 0){
stop("Tumor_Sample_Label was not found. Please check clinical data or let use.tumorSampleLabel be FALSE.")
}
ctsb <- treeData$sample[!treeData$sample %in% c("NORMAL", "internal node")]
treeData$sample[!treeData$sample %in% c("NORMAL", "internal node")] <- lapply(as.list(ctsb),
function(x){
tsb.label[which(tsb.label$Tumor_Sample_Barcode == x),]$Tumor_Sample_Label
}) %>%
unlist()
}
if(!is.null(branchCol)){
if(branchCol == "mutSig"){
color_scale <- getSigColors(as.character(unique(treeData$Signature)) )
sig_level <- levels(treeData$Signature)
## initialize
Signature <- NULL
p <- p + geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2, color = Signature), size=segmentSize)
p <- p + scale_color_manual(breaks = sig_level ,values = color_scale) +
theme(legend.title = element_text(size = legend.title.size))
}else{
## sort branch tumor type
all.types <- unique(treeData$Mutation_Type)
public <- all.types[grep("Public", all.types)]
shared <- all.types[grep("Shared", all.types)]
private <- all.types[grep("Private", all.types)]
type.level <- c(public, shared, private)
## get colors
type_all_colors <- c("#7fc97f","#fdc086", "#E64B35FF", "#82166E",
"#B77B42","#6349B7","#D5017D","#B77562",
"#88A4FF", "#439F18", "#971D37","#8C9F3C")
if(length(type.level) > length(type_all_colors)){
left_colors <- sample(colors(),
length(type.level)-length(type_all_colors),
replace = FALSE)
type_all_colors <- append(type_all_colors,left_colors)
}else{
type_all_colors <- type_all_colors[seq_len(length(type.level)) ]
}
names(type_all_colors) <- type.level
treeData$Mutation_Type <- factor(treeData$Mutation_Type, levels = type.level)
## initialize
Mutation_Type <- NULL
p <- p + geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2, color = factor(Mutation_Type)), size=segmentSize)
## remove legend title
p <- p + theme(legend.title = element_blank()) +
scale_color_manual(breaks = type.level, values = type_all_colors)
}
}else{
if(compare){
# p <- p + geom_point(aes(x = x1, y = y1),
# color = uncommon.col,
# data = treeData[treeData$is.match == "NO"],
# size = samplePointsSize)
p <- p + geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2),
color = "black",
data = treeData ,size = segmentSize)
}else{
p <- p + geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2), color = "black",
data = treeData,
size=segmentSize, show.legend = FALSE)
}
}
p <- p + theme(axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.line = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
panel.border = element_blank(),
legend.text = element_text(size = legend.text.size),
legend.position = 'right') +
scale_x_discrete(expand = expansion(add = mean(treeData$distance)))+
coord_fixed(ratio= 1)
if(compare & is.null(branchCol)){
unmatch_dat <- treeData[treeData$is.match == "NO"&treeData$sample == 'internal node',]
if(nrow(unmatch_dat) > 1){
unmatch_label <- paste0("Clades absent in ", compare.tree.name)
}else if(nrow(unmatch_dat) == 1){
unmatch_label <- paste0("Clade absent in ", compare.tree.name)
}
if(nrow(unmatch_dat) > 0){
rootLabel <- "NORMAL"
numNORMAL <- which(tree$tip.label == rootLabel)
tree <- ape::root(tree, tree$tip.label[numNORMAL])
treeEdge <- data.table::data.table(node = tree$edge[,1], endNum = tree$edge[,2])
rootRow <- which(treeEdge$endNum == numNORMAL)
rootNode <- treeEdge$node[rootRow]
node_point_list <- c()
for(i in seq_len(nrow(unmatch_dat))){
point <- unmatch_dat$end_num[i]
for(i1 in seq_len(length(tree$tip.label))){
tip <- tree$tip.label[i1]
if(tip == rootLabel){
next
}
start <- i1
tip_point_list <- c(i1)
while(TRUE){
end <- as.numeric(treeEdge[treeEdge$endNum == start,]$node)
tip_point_list <- c(tip_point_list, end)
if(end == point){
node_point_list <- unique(c(node_point_list, tip_point_list))
node_point_list <- node_point_list[node_point_list!=point]
break
}
if(end == rootNode){
break
}
start <- end
}
}
}
p <- p + geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2),
color = uncommon.col,
data = treeData[treeData$end_num %in% node_point_list,],
size = segmentSize)
}
}
p <- p + geom_point(aes(x = x2,y = y2),
data = treeData[treeData$sample == 'internal node',],
size = nodePointsSize, color = "#8c510a", fill = "white", shape = 21,
stroke = nodeStrokeSize)
p <- p + geom_point(aes(x = x2, y = y2),
data = treeData[treeData$sample != 'internal node',],
size = samplePointsSize,color = "#67001F", fill = 'white', shape = 21, stroke = sampleStrokeSize)
if(compare & is.null(branchCol)){
if(nrow(unmatch_dat) > 0){
p <- p + geom_point(aes(x = x2, y = y2, color = is.match),
# color = uncommon.col,
data = unmatch_dat,
size = comparePointsSize) +
scale_color_manual(label = unmatch_label,values = uncommon.col) +
theme(legend.title = element_blank())
}
}
Nd <- treeData[sample == rootLabel,]$distance
if(length(Nd)!=0){
if(Nd != 0){
p <- p + geom_point(aes(x =0 , y = 0), size = nodePointsSize, color = "#8c510a",
fill = 'white', shape = 21, stroke = nodeStrokeSize)
}
}
p <- p + geom_text_repel(aes(x = x2, y = y2, label = sample),
nudge_y = textAdjust/10,
nudge_x = textAdjust/10,
segment.color = "grey",
segment.size = 0.25,
max.overlaps = 200,
data = treeData[(!treeData$sample %in% c("internal node",rootLabel)) &
treeData$x2 >= 0, ],
size = sampleTextSize ,force = 10)
p <- p + geom_text_repel(aes(x = x2, y = y2, label = sample),
nudge_y = textAdjust/10,
nudge_x = -textAdjust/10,
segment.color = "grey",
segment.size = 0.25,
max.overlaps = 200,
data = treeData[(!treeData$sample %in% c("internal node",rootLabel)) &
treeData$x2 < 0, ],
size = sampleTextSize ,force = 10)
## label NORMAL
p <- p + geom_text(aes(x = x2,y = y2-textAdjust/5),
label = rootLabel,
data = treeData[treeData$sample == rootLabel,],
size = sampleTextSize)
if(show.bootstrap){
p <- p + geom_label_repel(aes(x = x2, y = y2,label = boots),
data = bootsData,
# nudge_y = textAdjust/6,
# fontface = 'bold',
size = bootLabelSize,
box.padding = unit(bootPaddingSize, "lines"),
label.padding = bootLabelPaddingSize,
segment.colour = "grey50",
segment.size = 0.25,
max.overlaps = 200,
force = 5)
}
# if(compare){
# is.match <- NULL
# p <- p + geom_label_repel(aes(x = x1 + (x2-x1)/2 , y = y1 + (y2 - y1)/2, label = is.match),
# data = treeData[treeData$is.match != "NO",],
# fontface = 'bold',
# size = bootLabelSize, box.padding = unit(bootPaddingSize, "lines"), point.padding = unit(0.5, "lines"),
# segment.colour = "grey50", segment.size = 0.5, force = 5)
# }
#
tree.title <- patient
if(show.scale.bar){
mean_len <- ceiling(mean(treeData$distance))
if(!is.null(scale.bar.x)){
x_bar <- as.numeric(scale.bar.x)
}else{
x_bar <- min(treeData$x2) + min(treeData$x2)/2
}
if(!is.null(scale.bar.y)){
y_bar <- as.numeric(scale.bar.y)
}else{
y_bar <- (max(treeData$y2) - min(treeData$y2))/2
}
df_bar <- data.frame(
x = x_bar,
xend = x_bar,
y = y_bar - mean_len/2,
yend = y_bar + mean_len/2
)
p <- p +
geom_segment(
aes(x = x, y = y, xend = xend, yend = yend),
data = df_bar,
size = 0.6,
color = "black"
) +
geom_text(
aes(x = x, y = yend + mean_len/10),
data = df_bar,
label = as.character(mean_len),
size = 4,
color = "black"
)
}
if(!compare){
p <- p +
ggtitle(tree.title)+
theme(plot.title = element_text(face = "bold",colour = "black", hjust = 0.5,size = 13.5))
}
return(p)
}
getTreeData <- function(phyloTree = NULL,
branchCol = "mutSig",
signaturesRef = "cosmic_v2",
min.mut.count = 15,
compare = FALSE){
tree <- getTree(phyloTree)
rootLabel <- "NORMAL"
tree <- ape::root(tree, tree$tip.label[which(tree$tip.label == rootLabel)])
treeEdge <- data.table::data.table(node = tree$edge[,1], endNum = tree$edge[,2], length = tree$edge.length)
numNORMAL <- which(tree$tip.label == rootLabel)
#the position of NORMAL in edge
rootRow <- which(treeEdge$endNum == numNORMAL)
# the Node connected to NORMAL
rootNode <- treeEdge$node[rootRow]
rootEdge <- treeEdge$length[rootRow]
lb <- labelBranch(tree)
branchLabel <- lb[[1]]
subnumlist <- lb[[2]]
mainTrunk <- calMainTrunk(tree, treeEdge, rootNode = rootNode, rootLabel = rootLabel)
ccn <- calChildNodeNum(tree, treeEdge = treeEdge, mainTrunk, rootNode = rootNode)
numList <- ccn[[1]]
pointsList <- ccn[[2]]
nodeNoOnTree <- ccn[[3]]
nodeOnTree <- ccn[[4]]
if(length(mainTrunk) != 0){
gta <- getNodeAngle(tree = tree, treeEdge = treeEdge,
mainTrunk = mainTrunk, numList = numList,
pointsList = pointsList, rootNode = rootNode,
horizon = pi/2, W = pi, numNORMAL = numNORMAL)
adjacentWs <- gta[[1]]
adjacentAngles <- gta[[2]]
adjacentPoints <- gta[[3]]
}
#X1, y1 are the starting point, x2 and y2 are the end point.Horizon stands for the Angle between the line and the positive x axis;W stands for the Angle occupied by branches (explained in the paper)
#Distance means distance;Node is the internal node with the start;End_num is the number of each node in edge;
treeData <- data.table::data.table('x1'= 0, 'y1' = 0,
'x2' = 0, 'y2' = 0,
'horizon' = pi/2, 'W' = pi,
'angle' = 0,'distance' = 0,
'node' = rootNode, 'end_num' = rootNode)
treeData <- setPhyloTree(tree = tree, treeEdge = treeEdge, treeData = treeData, rootNode = rootNode,
mainTrunk = mainTrunk, adjacentPoints = adjacentPoints,
adjacentWs = adjacentWs, adjacentAngles = adjacentAngles,
horizon = pi/2, W = pi)
# subTrees <- ape::subtrees(tree)
# subroot <- unlist(lapply(subTrees, function(x){return(x$name)}))
if(length(nodeNoOnTree)!=0){
t <- 1
while(TRUE){
node <- nodeNoOnTree[t]
subnodes <- subnumlist[[node]]
subEdge <- treeEdge[treeEdge$node %in% subnodes, ]
rootNode <- node
numNORMAL <- NULL
mainTrunk <- calMainTrunk(tree, subEdge, rootNode = rootNode, rootLabel = "")
ccn <- calChildNodeNum(tree, subEdge, mainTrunk, rootNode = rootNode, ft = TRUE)
numList <- ccn[[1]]
pointsList <- ccn[[2]]
nodeNoOnTree <- append(nodeNoOnTree,ccn[[3]])
nodeOnTree <- append(nodeOnTree,ccn[[4]])
horizon <- treeData[treeData$end_num == rootNode,]$angle
W <- treeData[treeData$end_num == rootNode,]$W
gta <- getNodeAngle(tree = tree, treeEdge = subEdge,
mainTrunk = mainTrunk, numList = numList,
pointsList = pointsList, rootNode = rootNode,
horizon = horizon, W = W, numNORMAL = NULL)
adjacentWs <- gta[[1]]
adjacentAngles <- gta[[2]]
adjacentPoints <- gta[[3]]
treeData <- setPhyloTree(tree = tree, treeEdge = subEdge, treeData = treeData, rootNode = rootNode,
mainTrunk = mainTrunk, adjacentPoints = adjacentPoints,
adjacentWs = adjacentWs, adjacentAngles = adjacentAngles,
horizon = horizon, W = W)
if(length(nodeOnTree) == tree$Nnode){
break
}
t <- t + 1
}
}
## bind NORMAL
rootNode <- treeData[1,]$node
treeData[1,]$end_num <- which(tree$tip.label == "NORMAL")
treeData[1,]$y2 <- -rootEdge
treeData[1,]$distance <- rootEdge
treeData$sample <- ""
sample_list <- vapply(treeData$end_num, function(e){
if(e > length(tree$tip.label)){
return("internal node")
}else{
return(tree$tip.label[e])
}
}, FUN.VALUE = character(1))
treeData$sample <- sample_list
# for(i in seq_len(nrow(treeData))){
# if(treeData$end_num[i] > length(tree$tip.label)){
# treeData$sample[i] <- "internal node"
# }
# else{
# pos <- treeData$end_num[i]
# treeData$sample[i] <- tree$tip.label[pos]
# }
# }
if(nrow(treeData) == 3){
angleList <- c(pi/3, 2*pi/3)
rows <- which(treeData$sample != "NORMAL")
treeData[rows, ]$angle <- angleList
treeData[rows, ]$x2 <- treeData[rows, ]$distance*cos(treeData[rows, ]$angle)
treeData[rows, ]$y2 <- treeData[rows, ]$distance*sin(treeData[rows, ]$angle)
}
if(nrow(treeData[treeData$x2 > 0,]) == 0){
maxy <- which.max(treeData$y2)
angle <- treeData$angle[maxy] - pi/18
treeData$angle[maxy] <- angle
treeData$y2[maxy] <- treeData$distance[maxy]*sin(angle) + treeData$y1[maxy]
treeData$x2[maxy] <- treeData$distance[maxy]*cos(angle) + treeData$x1[maxy]
}else if(nrow(treeData[treeData$x2 < 0,]) == 0){
maxy <- which.max(treeData$y2)
angle <- treeData$angle[maxy] + pi/18
treeData$angle[maxy] <- angle
treeData$y2[maxy] <- treeData$distance[maxy]*sin(angle) + treeData$y1[maxy]
treeData$x2[maxy] <- treeData$distance[maxy]*cos(angle) + treeData$x1[maxy]
}
## label represents the common evolution path of samples
treeData$label <- ""
label_list <- vapply(seq_len(nrow(treeData)), function(i){
if(treeData$sample[i] == "NORMAL"){
return(branchLabel[[rootNode]])
}else{
if(treeData$end_num[i] > length(tree$tip.label)){
return(branchLabel[[treeData$end_num[i]]])
}
else{
return(as.character(treeData$sample[i]))
}
}
}, FUN.VALUE = character(1))
treeData$label <- label_list
if(!is.null(branchCol) & !compare){
## add signature
if(branchCol == "mutSig"){
tri_matrix <- triMatrix(phyloTree,level = "4")
fit_out <- fitSignatures(tri_matrix,
signaturesRef = signaturesRef,
min.mut.count = min.mut.count)
if(is.na(fit_out)){
message("The number of mutations in all branches are less than min.mut.count in patient ",
phyloTree@patientID)
return(NA)
}
cos_sim_matrix <- fit_out[[1]]$cosine.similarity
sig_level <- colnames(cos_sim_matrix)
signatures <- apply(cos_sim_matrix,1,function(x)names(which.max(x)))
if(any(grepl("Signature", signatures))){
signatures <- gsub('Signature ', '', signatures)
sig_level <- gsub('Signature ', '', sig_level)
}else if(any(grepl("SBS", signatures))){
signatures <- gsub('SBS', '', signatures)
sig_level <- gsub('SBS', '', sig_level)
}
treeData$Signature <- signatures[treeData$label]
treeData[treeData$Signature == ''| is.na(treeData$Signature)]$Signature <- "Unknown"
if("Unknown" %in% treeData$Signature){
sig_level <- c(sig_level,"Unknown")
}
treeData$Signature <- factor(treeData$Signature, levels = sig_level)
treeData <- treeData[order(treeData$Signature), ]
}else{
branch_type <- getBranchType(phyloTree)
types <- as.character(branch_type$Mutation_Type)
names(types) <- as.character(branch_type$Branch_ID)
treeData$Mutation_Type <- types[treeData$label]
# dplyr::mutate(Mutation_Type = types[label]) %>%
# as.data.table()
}
}
treeData <- treeData[treeData$distance!= 0|treeData$sample == rootLabel,]
return(treeData)
}
setPhyloTree <- function(tree, treeEdge, treeData, rootNode, mainTrunk, adjacentPoints, adjacentWs, adjacentAngles,
horizon = pi/2 , W = pi){
trunkPath <- rev(c(mainTrunk,rootNode))
if(length(trunkPath) > 0){
# subdat_list <- lapply(2:length(trunkPath), function(i){
# x1 <- treeData[treeData$end_num == trunkPath[i-1],]$x2
# y1 <- treeData[treeData$end_num == trunkPath[i-1],]$y2
# W <- W
# angle <- horizon
# distance <- treeEdge[treeEdge$endNum == trunkPath[i],]$length
# if(horizon == pi/2){
# x2 <- x1
# y2 <- y1 + distance
# }else{
# x2 <- x1 + distance*cos(angle)
# y2 <- y1 + distance*sin(angle)
# }
# subdat <- data.table::data.table('x1'= x1, 'y1' = y1,
# 'x2' = x2, 'y2' = y2,
# 'horizon' = horizon,'W' = W,
# 'distance' = distance, 'angle' = horizon,
# 'node' = trunkPath[i-1],'end_num' = trunkPath[i])
# return(subdat)
# })
#
# subdat <- dplyr::bind_rows(subdat_list)
# treeData <- rbind(treeData, subdat)
for(i in 2:length(trunkPath)){
x1 <- treeData[treeData$end_num == trunkPath[i-1],]$x2
y1 <- treeData[treeData$end_num == trunkPath[i-1],]$y2
W <- W
angle <- horizon
distance <- treeEdge[treeEdge$endNum == trunkPath[i],]$length
if(horizon == pi/2){
x2 <- x1
y2 <- y1 + distance
}else{
x2 <- x1 + distance*cos(angle)
y2 <- y1 + distance*sin(angle)
}
subdat <- data.table::data.table('x1'= x1, 'y1' = y1,
'x2' = x2, 'y2' = y2,
'horizon' = horizon,'W' = W,
'distance' = distance, 'angle' = horizon,
'node' = trunkPath[i-1],'end_num' = trunkPath[i])
treeData <- rbind(treeData, subdat)
}
}
if(length(adjacentPoints) > 0){
# subdat_list <- lapply(seq_len(length(adjacentPoints)), function(i){
# point <- adjacentPoints[i]
# startnode <- treeEdge[treeEdge$endNum == point, ]$node
# x1 <- treeData[treeData$end_num == startnode,]$x2
# y1 <- treeData[treeData$end_num == startnode,]$y2
# W <- adjacentWs[point]
# angle <- adjacentAngles[point]
# distance <- treeEdge[treeEdge$endNum == point, ]$length
# x2 <- x1 + distance*cos(angle)
# y2 <- y1 + distance*sin(angle)
# subdat <- data.table::data.table('x1'= x1, 'y1' = y1,
# 'x2' = x2, 'y2' = y2,
# 'horizon' = horizon,'W' = W,
# 'distance' = distance, 'angle' = angle,
# 'node' = startnode,'end_num' = point)
# return(subdat)
# })
# subdat <- dplyr::bind_rows(subdat_list)
# treeData <- rbind(treeData, subdat)
for(i in seq_len(length(adjacentPoints))){
point <- adjacentPoints[i]
startnode <- treeEdge[treeEdge$endNum == point, ]$node
x1 <- treeData[treeData$end_num == startnode,]$x2
y1 <- treeData[treeData$end_num == startnode,]$y2
W <- adjacentWs[point]
angle <- adjacentAngles[point]
distance <- treeEdge[treeEdge$endNum == point, ]$length
x2 <- x1 + distance*cos(angle)
y2 <- y1 + distance*sin(angle)
subdat <- data.table::data.table('x1'= x1, 'y1' = y1,
'x2' = x2, 'y2' = y2,
'horizon' = horizon,'W' = W,
'distance' = distance, 'angle' = angle,
'node' = startnode,'end_num' = point)
treeData <- rbind(treeData, subdat)
}
}
return(treeData)
}
calMainTrunk <- function(tree, treeEdge, rootNode, rootLabel = "NORMAL"){
Ntips <- Ntip(tree)
subtips <- treeEdge[!treeEdge$endNum > Ntips,]$endNum
distanceTable <- c(0)
# for(i in seq_len(length(subtips))){
# distanceTable <- cbind(distanceTable, -1)
# }
s <- vapply(seq_len(length(subtips)), function(i)-1, FUN.VALUE = numeric(1))
distanceTable <- c(distanceTable, s)
names(distanceTable) <- c("x",s)
distanceTable <- distanceTable[(-1)]
names(distanceTable) <- tree$tip.label[subtips]
t <- 1
path <- list()
while(t <= length(distanceTable)){
subPath <- c()
end <- which(tree$tip.label == names(distanceTable)[t])
name <- names(distanceTable)[t]
dis <- 0
while (TRUE) {
p <- which(treeEdge$endNum == end)
dis <- dis + treeEdge$length[p]
start <- treeEdge$node[p]
subPath <- append(subPath,end)
end <- start
if(start == rootNode){
break
}
}
path[[name]] <- subPath
distanceTable[name] <- dis
t <- t + 1
}
if(rootLabel == "NORMAL"){
distanceTable <- distanceTable[names(distanceTable) != rootLabel]
}
result <- path[[names(which.max(distanceTable))]]
}
getNodeAngle <- function(tree, treeEdge, mainTrunk,
numList, pointsList, rootNode,
horizon = pi/2, W = pi ,numNORMAL = NULL){
# print(pointsList)
if(!is.null(numNORMAL)){
treeEdge <- treeEdge[treeEdge$endNum != numNORMAL,]
}
nodes <- factor(rev(c(mainTrunk,rootNode)), levels = rev(c(mainTrunk,rootNode)))
adjacentPoints <- c()
adjacentPoints <- lapply(nodes, function(n){
p <- treeEdge[treeEdge$node == n,]$endNum %>% setdiff(nodes)
return(p)
}) %>% unlist()
# for(n in nodes){
# p <- treeEdge[treeEdge$node == n,]$endNum %>% setdiff(nodes)
# adjacentPoints <- append(adjacentPoints,p)
# }
left <- c()
right <- c()
leftList <- c()
rightList <- c()
i <- 1
# adresult <- lapply(seq_len(length(adjacentPoints)),function(i){
# point <- adjacentPoints[i]
# pos <- which(pointsList == point)
# if(i == 1){
# if(horizon < pi/2){
# return(list(l = numList[pos], ll = point, r = NULL, rl = NULL))
# }else{
# return(list(l = NULL, ll = NULL, r = numList[pos], rl = point))
# }
# }else{
# if(sum(right) <= sum(left)){
# return(list(l = NULL, ll = NULL, r = numList[pos], rl = point))
# }else{
# return(list(l = numList[pos], ll = point, r = NULL, rl = NULL))
# }
# }
# })
# left <- lapply(adresult, function(x)x$l) %>% unlist()
# leftList <- lapply(adresult, function(x)x$ll) %>% unlist()
# right <- lapply(adresult, function(x)x$r) %>% unlist()
# rightList <- lapply(adresult, function(x)x$rl) %>% unlist()
for(point in adjacentPoints){
if(i == 1){
if(horizon < pi/2){
leftList <- append(leftList, point)
pos <- which(pointsList == point)
left <- append(left,numList[pos])
}else{
rightList <- append(rightList, point)
pos <- which(pointsList == point)
right <- append(right,numList[pos])
}
}else{
if(sum(right) <= sum(left)){
rightList <- append(rightList, point)
pos <- which(pointsList == point)
right <- append(right,numList[pos])
}else{
leftList <- append(leftList, point)
pos <- which(pointsList == point)
left <- append(left,numList[pos])
}
}
i <- i+1
}
adjacentWs <- rep(0,length(tree$edge.length)+1)
adjacentAngles <- rep(0,length(tree$edge.length)+1)
if(length(rightList) > 0){
startr <- horizon - W/2
wrt <- W/2
totalR <- sum(right)
for(i in seq_len(length(rightList))){
if(rightList[i] <= length(tree$tip.label)){
n <- 1
}else{
n <- right[i]
}
wr <- wrt*(n/totalR)
adjacentWs[rightList[i]] <- wr
if(i == 1){
angler <- wr/2 + startr
c1 <- horizon != pi/2
c2 <- any(startr != 0,n == 1)
c3 <- any(startr != pi/2,n == 1)
if(all(c1,c2,c3)){
if(length(c(leftList,rightList))!=1){
angler <- startr
}
}
adjacentAngles[rightList[i]] <- angler
startr <- angler
}else{
if(n == 1 & horizon != pi/2){
startr <- startr - wr/2
}
angler <- wr/2 + startr + adjacentWs[rightList[i-1]]/2
adjacentAngles[rightList[i]] <- angler
startr <- angler
}
}
}
if(length(leftList) > 0){
wlt <- W/2
startl <- horizon + wlt
totalL <- sum(left)
for(i in seq_len(length(leftList))){
if(leftList[i] <= length(tree$tip.label)){
n <- 1
}else{
n <- left[i]
}
wl <- wlt*(n/totalL)
adjacentWs[leftList[i]] <- wl
if(i == 1){
anglel <- startl - wl/2
c1 <- horizon != pi/2
c2 <- any(startl != pi, n == 1)
c3 <- any(startl != pi/2, n == 1)
# angle extension
if(all(c1,c2,c3)){
if(length(c(leftList,rightList))!=1){
anglel <- startl
}
}
adjacentAngles[leftList[i]] <- anglel
startl <- anglel
}else{
if(n == 1& horizon != pi/2){
startl <- startl + wl/2
}
anglel <- startl- wl/2 - adjacentWs[leftList[i-1]]/2
adjacentAngles[leftList[i]] <- anglel
startl <- anglel
}
}
}
return(list(adjacentWs,adjacentAngles, adjacentPoints))
}
calChildNodeNum <- function(tree, treeEdge, mainTrunk, rootNode, ft = FALSE){
pointsList <- sort(unique(c(treeEdge$node,treeEdge$endNum)))
numList <- rep(1, length(pointsList))
lrnumList <- rep(0, length(pointsList))
# for(i in seq_len(length(pointsList))){
# numList[i] <- 1
# lrnumList[i] <- 0
# }
for(i in seq_len(length(pointsList))){
if(pointsList[i] == rootNode){
next
}
start <- pointsList[i]
end <- as.numeric(treeEdge[treeEdge$endNum == start,]$node)
while(TRUE){
pos <- which(pointsList == end)
numList[pos] <- numList[pos] + 1
if(end == rootNode){
break
}
start <- treeEdge[treeEdge$endNum == end, ]$endNum
end <- treeEdge[treeEdge$endNum == end, ]$node
}
# if(count > 1){
# numList[i] <- 0
# }
}
## samples are not in mainTrunk
nodeRange <- c(Ntip(tree)+1):(length(tree$edge.length)+1)
nodeOnTree <- c(mainTrunk,rootNode) %>% intersect(nodeRange)
nodeNoOnTree <- treeEdge[treeEdge$node %in% c(mainTrunk,rootNode), ]$endNum %>% setdiff(nodeOnTree) %>% intersect(nodeRange)
return(list(numList,pointsList, nodeNoOnTree, nodeOnTree))
}
getSigColors <- function(signatures){
signature_colors <- c("1" = "#E41A1C", "2" = "#377EB8","3" = "#7F0000",
"4" = "#35978f","5" = "#A6D854", "6" = "#FFD92F",
"7" = "#E5C494", "8" = "#8DD3C7","9" = "#6E016B",
"10" = "#BEBADA", "11" = "#e08214", "12" = "#80B1D3",
"13" = "#d6604d", "14" = "#ffff99","15" = "#FCCDE5",
"16" = "#FF6A5A","17" = "#BC80BD","18" = "#CCEBC5",
"19" = "#fb9a99", "20" = "#B6646A", "21" = "#9F994E",
"22" = "#7570B3" , "23" = "#c51b7d" , "24" = "#66A61E" ,
"25" = "#E6AB02" ,"26" = "#003c30", "27" = "#666666",
"28" = "#524762","29" = "#8926AC","30" = "#C42669",
"1A" = "#47B39F","1B" = "#B84C60", "R1" = "#B93183",
"R2" = "#B93183", "R3" = "#845AD1", "U1" = "#7BA52E",
"U2" = "#7A4287","7a" = "#663399", "7b" = "#AD8B5E",
"7c" = "#194682", "7d" = "#FF8931", "10a" = "#AD3D07",
"10b" = "#0D5DFF", "17a" = "#F2A200", "17b" = "#8C3B45",
"31" = "#009F9A", "32" = "#FF2B32", "33" = "#153E71",
"34" = "#FF8C2F", "35" ="#2F8CFF", "36" = "#62BA1F",
"37" = "#9745D6", "38" = "#3826CA", "39" = "#96C4BC",
"40" = "#982B9F", "41" = "#9C3C33", "42" = "#3EB15A",
"43" = "#EF1297", "44" = "#B1098C", "45" = "#883E93",
"46" = "#6ED66B", "47" ="#9F679F" , "48" = "#5B309F",
"49" = "#4E7A9F", "50" = "#9F6878", "51" = "#154E0F",
"52" = "#9D887D", "53" = "#143045", "54" = "#8F1B1D",
"55" = "#544C9F", "56" = "#7F4E10", "57" = "#429F74",
"58" = "#6E6A5F", "59" = "#235F9F", "60" = "#8A1D1B",
"84" = "#4B9F34", "85" = "#598F4D", "Unknown" = "black")
color_scale <- signature_colors[signatures]
return(color_scale)
}
labelBranch <- function(tree){
Root <- which(tree$tip.label == "NORMAL")
internalNodes <- sort(unique(tree$edge[,1]))
result <- list()
subnumList <- list()
subnumList <- lapply(seq_len(length(tree$edge.length)+1), function(i){
if(i > Ntip(tree)){
return(i)
}
else{
return(NA)
}
})
result <- rep(NA, (length(tree$edge.length)+1)) %>% as.list()
# for(i in seq_len(length(tree$edge.length)+1)){
# result[[i]] <- NA
# if(i > Ntip(tree)){
# subnumList[[i]] <- i
# }
# else{
# subnumList[[i]] <- NA
# }
# }
# print(result)
# print(subnumList)
end <- tree$edge[which(tree$edge[,2] == Root),1]
for(i in seq_len(length(tree$tip.label))){
row <- tree$edge[which(tree$edge[,2] == i), ]
if(i == Root){
next
}
while (TRUE) {
node <- row[1]
if(node == end){
result[[node]] <- append(result[[node]], i)
break
}
else{
result[[node]] <- append(result[[node]], i)
row <- tree$edge[which(tree$edge[,2] == node),]
}
}
}
rootNode <- tree$edge[which(tree$edge[,2] == Root),1]
for(i in (length(tree$tip.label)+1):(length(tree$edge.length)+1)){
if(i == rootNode){
next
}
row <- tree$edge[which(tree$edge[,2] == i), ]
while (TRUE) {
node <- row[1]
if(node == rootNode){
subnumList[[node]] <- append(subnumList[[node]], i)
break
}
else{
subnumList[[node]] <- append(subnumList[[node]], i)
row <- tree$edge[which(tree$edge[,2] == node),]
}
}
}
g <- lapply(result, function(x){return(paste(sort(tree$tip.label[x[-1]], decreasing = TRUE),collapse = "&"))})
return(list(g,subnumList))
}
Try the MesKit package in your browser
Any scripts or data that you put into this service are public.
MesKit documentation built on March 28, 2021, 6 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions labelBranch getSigColors calChildNodeNum getNodeAngle calMainTrunk setPhyloTree getTreeData plotTree
plotTree <- function(phyloTree,
treeData = NULL,
branchCol = "mutType",
show.bootstrap = TRUE,
signaturesRef = "cosmic_v2",
min.mut.count = 15,
min.ratio = 1/20,
uncommon.col = "red",
compare.tree.name = "",
use.tumorSampleLabel = FALSE,
show.scale.bar = FALSE,
scale.bar.x = NULL,
scale.bar.y = NULL){
if(min.ratio <= 0|min.ratio > 1){
stop("min.ratio should be within (0,1]")
}
## adjust length of branches by min.ratio
tree <- getTree(phyloTree)
min.len <- max(tree$edge.length)*min.ratio
tree$edge.length[tree$edge.length < min.len] <- min.len
phyloTree <- new('phyloTree',
patientID = getPhyloTreePatient(phyloTree),
tree = tree,
binary.matrix = getBinaryMatrix(phyloTree),
ccf.matrix = getCCFMatrix(phyloTree),
mut.branches = getMutBranches(phyloTree),
branch.type = getBranchType(phyloTree),
ref.build = getPhyloTreeRef(phyloTree),
bootstrap.value = getBootstrapValue(phyloTree),
method = getTreeMethod(phyloTree),
tsb.label = getPhyloTreeTsbLabel(phyloTree))
patient <- getPhyloTreePatient(phyloTree)
if(is.null(treeData)){
treeData <- getTreeData(phyloTree = phyloTree,
branchCol = branchCol,
signaturesRef = signaturesRef,
min.mut.count = min.mut.count)
if(identical(treeData, NA)){
return(NA)
}
compare <- FALSE
}else{
compare <- TRUE
}
# set.seed(1234)
## get bootstrap value
boot_value <- getBootstrapValue(phyloTree)
rootLabel <- "NORMAL"
## plot phylotree
samplePointsSize <- 3
comparePointsSize <- 3.5
sampleTextSize <- 3.5
nodePointsSize <- 1.7
segmentSize <- 1.5
nodeStrokeSize <- 0.25
sampleStrokeSize <- 1.5
bootLabelSize <- 3
bootPaddingSize <- 0.35
bootLabelPaddingSize <- 0.15
legend.text.size <- 10
legend.title.size <- legend.text.size + 0.5
samplesLength <- nrow(treeData[treeData$sample != "internal node",])
if(samplesLength > 7){
samplePointsSize <- 1.5
comparePointsSize <- 2
sampleTextSize <- 3
segmentSize <- 0.8
nodePointsSize <- 0.8
nodeStrokeSize <- 0.15
sampleStrokeSize <- 0.8
bootLabelSize <- 2.5
bootPaddingSize <- 0.1
bootLabelPaddingSize <- 0.1
legend.text.size <- 9
legend.title.size <- legend.text.size + 0.5
}
rootNode <- treeData[treeData$sample == rootLabel,]$node
if(length(boot_value) == 1){
bootsData <- data.frame(x2 = 0, y2 = 0, node = rootNode, end_num = rootNode, boots = boot_value)
}else{
sub <- data.table::data.table(x2 = 0, y2 = 0,node = rootNode, end_num = rootNode)
bootsData <- rbind(treeData[treeData$sample == 'internal node',] %>%
dplyr::select("x2","y2","node","end_num"),sub)
boots <- c()
LN <- min(bootsData$node)-1
boots_list <- lapply(seq_len(nrow(bootsData)), function(i){
if(i == nrow(bootsData)){
boots <- append(boots, boot_value[rootNode - LN])
}else{
boots <- append(boots, boot_value[bootsData$end_num[i] - LN])
}
return(boots)
}) %>% unlist()
bootsData <- cbind(bootsData, boots = boots_list)
}
## initialize
x2 <- NULL
y2 <- NULL
x1 <- NULL
y1 <- NULL
## get the max value of X axis
x_max <- max(abs(treeData$x2))
# print(treeData)
p <- ggplot(data = treeData) +
## balance the space on the left and right sides
geom_segment(aes(x = 0, y = 0, xend = x_max, yend = 0),color = "white",size = 0.01)+
geom_segment(aes(x = 0, y = 0, xend = -x_max, yend = 0),color = "white",size = 0.01)
textAdjust <- mean(as.numeric(treeData$distance))
if(use.tumorSampleLabel){
tsb.label <- getPhyloTreeTsbLabel(phyloTree)
if(nrow(tsb.label) == 0){
stop("Tumor_Sample_Label was not found. Please check clinical data or let use.tumorSampleLabel be FALSE.")
}
ctsb <- treeData$sample[!treeData$sample %in% c("NORMAL", "internal node")]
treeData$sample[!treeData$sample %in% c("NORMAL", "internal node")] <- lapply(as.list(ctsb),
function(x){
tsb.label[which(tsb.label$Tumor_Sample_Barcode == x),]$Tumor_Sample_Label
}) %>%
unlist()
}
if(!is.null(branchCol)){
if(branchCol == "mutSig"){
color_scale <- getSigColors(as.character(unique(treeData$Signature)) )
sig_level <- levels(treeData$Signature)
## initialize
Signature <- NULL
p <- p + geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2, color = Signature), size=segmentSize)
p <- p + scale_color_manual(breaks = sig_level ,values = color_scale) +
theme(legend.title = element_text(size = legend.title.size))
}else{
## sort branch tumor type
all.types <- unique(treeData$Mutation_Type)
public <- all.types[grep("Public", all.types)]
shared <- all.types[grep("Shared", all.types)]
private <- all.types[grep("Private", all.types)]
type.level <- c(public, shared, private)
## get colors
type_all_colors <- c("#7fc97f","#fdc086", "#E64B35FF", "#82166E",
"#B77B42","#6349B7","#D5017D","#B77562",
"#88A4FF", "#439F18", "#971D37","#8C9F3C")
if(length(type.level) > length(type_all_colors)){
left_colors <- sample(colors(),
length(type.level)-length(type_all_colors),
replace = FALSE)
type_all_colors <- append(type_all_colors,left_colors)
}else{
type_all_colors <- type_all_colors[seq_len(length(type.level)) ]
}
names(type_all_colors) <- type.level
treeData$Mutation_Type <- factor(treeData$Mutation_Type, levels = type.level)
## initialize
Mutation_Type <- NULL
p <- p + geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2, color = factor(Mutation_Type)), size=segmentSize)
## remove legend title
p <- p + theme(legend.title = element_blank()) +
scale_color_manual(breaks = type.level, values = type_all_colors)
}
}else{
if(compare){
# p <- p + geom_point(aes(x = x1, y = y1),
# color = uncommon.col,
# data = treeData[treeData$is.match == "NO"],
# size = samplePointsSize)
p <- p + geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2),
color = "black",
data = treeData ,size = segmentSize)
}else{
p <- p + geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2), color = "black",
data = treeData,
size=segmentSize, show.legend = FALSE)
}
}
p <- p + theme(axis.title.x = element_blank(),
axis.text.x = element_blank(),
axis.ticks.x = element_blank(),
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.line = element_blank(),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
panel.border = element_blank(),
legend.text = element_text(size = legend.text.size),
legend.position = 'right') +
scale_x_discrete(expand = expansion(add = mean(treeData$distance)))+
coord_fixed(ratio= 1)
if(compare & is.null(branchCol)){
unmatch_dat <- treeData[treeData$is.match == "NO"&treeData$sample == 'internal node',]
if(nrow(unmatch_dat) > 1){
unmatch_label <- paste0("Clades absent in ", compare.tree.name)
}else if(nrow(unmatch_dat) == 1){
unmatch_label <- paste0("Clade absent in ", compare.tree.name)
}
if(nrow(unmatch_dat) > 0){
rootLabel <- "NORMAL"
numNORMAL <- which(tree$tip.label == rootLabel)
tree <- ape::root(tree, tree$tip.label[numNORMAL])
treeEdge <- data.table::data.table(node = tree$edge[,1], endNum = tree$edge[,2])
rootRow <- which(treeEdge$endNum == numNORMAL)
rootNode <- treeEdge$node[rootRow]
node_point_list <- c()
for(i in seq_len(nrow(unmatch_dat))){
point <- unmatch_dat$end_num[i]
for(i1 in seq_len(length(tree$tip.label))){
tip <- tree$tip.label[i1]
if(tip == rootLabel){
next
}
start <- i1
tip_point_list <- c(i1)
while(TRUE){
end <- as.numeric(treeEdge[treeEdge$endNum == start,]$node)
tip_point_list <- c(tip_point_list, end)
if(end == point){
node_point_list <- unique(c(node_point_list, tip_point_list))
node_point_list <- node_point_list[node_point_list!=point]
break
}
if(end == rootNode){
break
}
start <- end
}
}
}
p <- p + geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2),
color = uncommon.col,
data = treeData[treeData$end_num %in% node_point_list,],
size = segmentSize)
}
}
p <- p + geom_point(aes(x = x2,y = y2),
data = treeData[treeData$sample == 'internal node',],
size = nodePointsSize, color = "#8c510a", fill = "white", shape = 21,
stroke = nodeStrokeSize)
p <- p + geom_point(aes(x = x2, y = y2),
data = treeData[treeData$sample != 'internal node',],
size = samplePointsSize,color = "#67001F", fill = 'white', shape = 21, stroke = sampleStrokeSize)
if(compare & is.null(branchCol)){
if(nrow(unmatch_dat) > 0){
p <- p + geom_point(aes(x = x2, y = y2, color = is.match),
# color = uncommon.col,
data = unmatch_dat,
size = comparePointsSize) +
scale_color_manual(label = unmatch_label,values = uncommon.col) +
theme(legend.title = element_blank())
}
}
Nd <- treeData[sample == rootLabel,]$distance
if(length(Nd)!=0){
if(Nd != 0){
p <- p + geom_point(aes(x =0 , y = 0), size = nodePointsSize, color = "#8c510a",
fill = 'white', shape = 21, stroke = nodeStrokeSize)
}
}
p <- p + geom_text_repel(aes(x = x2, y = y2, label = sample),
nudge_y = textAdjust/10,
nudge_x = textAdjust/10,
segment.color = "grey",
segment.size = 0.25,
max.overlaps = 200,
data = treeData[(!treeData$sample %in% c("internal node",rootLabel)) &
treeData$x2 >= 0, ],
size = sampleTextSize ,force = 10)
p <- p + geom_text_repel(aes(x = x2, y = y2, label = sample),
nudge_y = textAdjust/10,
nudge_x = -textAdjust/10,
segment.color = "grey",
segment.size = 0.25,
max.overlaps = 200,
data = treeData[(!treeData$sample %in% c("internal node",rootLabel)) &
treeData$x2 < 0, ],
size = sampleTextSize ,force = 10)
## label NORMAL
p <- p + geom_text(aes(x = x2,y = y2-textAdjust/5),
label = rootLabel,
data = treeData[treeData$sample == rootLabel,],
size = sampleTextSize)
if(show.bootstrap){
p <- p + geom_label_repel(aes(x = x2, y = y2,label = boots),
data = bootsData,
# nudge_y = textAdjust/6,
# fontface = 'bold',
size = bootLabelSize,
box.padding = unit(bootPaddingSize, "lines"),
label.padding = bootLabelPaddingSize,
segment.colour = "grey50",
segment.size = 0.25,
max.overlaps = 200,
force = 5)
}
# if(compare){
# is.match <- NULL
# p <- p + geom_label_repel(aes(x = x1 + (x2-x1)/2 , y = y1 + (y2 - y1)/2, label = is.match),
# data = treeData[treeData$is.match != "NO",],
# fontface = 'bold',
# size = bootLabelSize, box.padding = unit(bootPaddingSize, "lines"), point.padding = unit(0.5, "lines"),
# segment.colour = "grey50", segment.size = 0.5, force = 5)
# }
#
tree.title <- patient
if(show.scale.bar){
mean_len <- ceiling(mean(treeData$distance))
if(!is.null(scale.bar.x)){
x_bar <- as.numeric(scale.bar.x)
}else{
x_bar <- min(treeData$x2) + min(treeData$x2)/2
}
if(!is.null(scale.bar.y)){
y_bar <- as.numeric(scale.bar.y)
}else{
y_bar <- (max(treeData$y2) - min(treeData$y2))/2
}
df_bar <- data.frame(
x = x_bar,
xend = x_bar,
y = y_bar - mean_len/2,
yend = y_bar + mean_len/2
)
p <- p +
geom_segment(
aes(x = x, y = y, xend = xend, yend = yend),
data = df_bar,
size = 0.6,
color = "black"
) +
geom_text(
aes(x = x, y = yend + mean_len/10),
data = df_bar,
label = as.character(mean_len),
size = 4,
color = "black"
)
}
if(!compare){
p <- p +
ggtitle(tree.title)+
theme(plot.title = element_text(face = "bold",colour = "black", hjust = 0.5,size = 13.5))
}
return(p)
}
getTreeData <- function(phyloTree = NULL,
branchCol = "mutSig",
signaturesRef = "cosmic_v2",
min.mut.count = 15,
compare = FALSE){
tree <- getTree(phyloTree)
rootLabel <- "NORMAL"
tree <- ape::root(tree, tree$tip.label[which(tree$tip.label == rootLabel)])
treeEdge <- data.table::data.table(node = tree$edge[,1], endNum = tree$edge[,2], length = tree$edge.length)
numNORMAL <- which(tree$tip.label == rootLabel)
#the position of NORMAL in edge
rootRow <- which(treeEdge$endNum == numNORMAL)
# the Node connected to NORMAL
rootNode <- treeEdge$node[rootRow]
rootEdge <- treeEdge$length[rootRow]
lb <- labelBranch(tree)
branchLabel <- lb[[1]]
subnumlist <- lb[[2]]
mainTrunk <- calMainTrunk(tree, treeEdge, rootNode = rootNode, rootLabel = rootLabel)
ccn <- calChildNodeNum(tree, treeEdge = treeEdge, mainTrunk, rootNode = rootNode)
numList <- ccn[[1]]
pointsList <- ccn[[2]]
nodeNoOnTree <- ccn[[3]]
nodeOnTree <- ccn[[4]]
if(length(mainTrunk) != 0){
gta <- getNodeAngle(tree = tree, treeEdge = treeEdge,
mainTrunk = mainTrunk, numList = numList,
pointsList = pointsList, rootNode = rootNode,
horizon = pi/2, W = pi, numNORMAL = numNORMAL)
adjacentWs <- gta[[1]]
adjacentAngles <- gta[[2]]
adjacentPoints <- gta[[3]]
}
#X1, y1 are the starting point, x2 and y2 are the end point.Horizon stands for the Angle between the line and the positive x axis;W stands for the Angle occupied by branches (explained in the paper)
#Distance means distance;Node is the internal node with the start;End_num is the number of each node in edge;
treeData <- data.table::data.table('x1'= 0, 'y1' = 0,
'x2' = 0, 'y2' = 0,
'horizon' = pi/2, 'W' = pi,
'angle' = 0,'distance' = 0,
'node' = rootNode, 'end_num' = rootNode)
treeData <- setPhyloTree(tree = tree, treeEdge = treeEdge, treeData = treeData, rootNode = rootNode,
mainTrunk = mainTrunk, adjacentPoints = adjacentPoints,
adjacentWs = adjacentWs, adjacentAngles = adjacentAngles,
horizon = pi/2, W = pi)
# subTrees <- ape::subtrees(tree)
# subroot <- unlist(lapply(subTrees, function(x){return(x$name)}))
if(length(nodeNoOnTree)!=0){
t <- 1
while(TRUE){
node <- nodeNoOnTree[t]
subnodes <- subnumlist[[node]]
subEdge <- treeEdge[treeEdge$node %in% subnodes, ]
rootNode <- node
numNORMAL <- NULL
mainTrunk <- calMainTrunk(tree, subEdge, rootNode = rootNode, rootLabel = "")
ccn <- calChildNodeNum(tree, subEdge, mainTrunk, rootNode = rootNode, ft = TRUE)
numList <- ccn[[1]]
pointsList <- ccn[[2]]
nodeNoOnTree <- append(nodeNoOnTree,ccn[[3]])
nodeOnTree <- append(nodeOnTree,ccn[[4]])
horizon <- treeData[treeData$end_num == rootNode,]$angle
W <- treeData[treeData$end_num == rootNode,]$W
gta <- getNodeAngle(tree = tree, treeEdge = subEdge,
mainTrunk = mainTrunk, numList = numList,
pointsList = pointsList, rootNode = rootNode,
horizon = horizon, W = W, numNORMAL = NULL)
adjacentWs <- gta[[1]]
adjacentAngles <- gta[[2]]
adjacentPoints <- gta[[3]]
treeData <- setPhyloTree(tree = tree, treeEdge = subEdge, treeData = treeData, rootNode = rootNode,
mainTrunk = mainTrunk, adjacentPoints = adjacentPoints,
adjacentWs = adjacentWs, adjacentAngles = adjacentAngles,
horizon = horizon, W = W)
if(length(nodeOnTree) == tree$Nnode){
break
}
t <- t + 1
}
}
## bind NORMAL
rootNode <- treeData[1,]$node
treeData[1,]$end_num <- which(tree$tip.label == "NORMAL")
treeData[1,]$y2 <- -rootEdge
treeData[1,]$distance <- rootEdge
treeData$sample <- ""
sample_list <- vapply(treeData$end_num, function(e){
if(e > length(tree$tip.label)){
return("internal node")
}else{
return(tree$tip.label[e])
}
}, FUN.VALUE = character(1))
treeData$sample <- sample_list
# for(i in seq_len(nrow(treeData))){
# if(treeData$end_num[i] > length(tree$tip.label)){
# treeData$sample[i] <- "internal node"
# }
# else{
# pos <- treeData$end_num[i]
# treeData$sample[i] <- tree$tip.label[pos]
# }
# }
if(nrow(treeData) == 3){
angleList <- c(pi/3, 2*pi/3)
rows <- which(treeData$sample != "NORMAL")
treeData[rows, ]$angle <- angleList
treeData[rows, ]$x2 <- treeData[rows, ]$distance*cos(treeData[rows, ]$angle)
treeData[rows, ]$y2 <- treeData[rows, ]$distance*sin(treeData[rows, ]$angle)
}
if(nrow(treeData[treeData$x2 > 0,]) == 0){
maxy <- which.max(treeData$y2)
angle <- treeData$angle[maxy] - pi/18
treeData$angle[maxy] <- angle
treeData$y2[maxy] <- treeData$distance[maxy]*sin(angle) + treeData$y1[maxy]
treeData$x2[maxy] <- treeData$distance[maxy]*cos(angle) + treeData$x1[maxy]
}else if(nrow(treeData[treeData$x2 < 0,]) == 0){
maxy <- which.max(treeData$y2)
angle <- treeData$angle[maxy] + pi/18
treeData$angle[maxy] <- angle
treeData$y2[maxy] <- treeData$distance[maxy]*sin(angle) + treeData$y1[maxy]
treeData$x2[maxy] <- treeData$distance[maxy]*cos(angle) + treeData$x1[maxy]
}
## label represents the common evolution path of samples
treeData$label <- ""
label_list <- vapply(seq_len(nrow(treeData)), function(i){
if(treeData$sample[i] == "NORMAL"){
return(branchLabel[[rootNode]])
}else{
if(treeData$end_num[i] > length(tree$tip.label)){
return(branchLabel[[treeData$end_num[i]]])
}
else{
return(as.character(treeData$sample[i]))
}
}
}, FUN.VALUE = character(1))
treeData$label <- label_list
if(!is.null(branchCol) & !compare){
## add signature
if(branchCol == "mutSig"){
tri_matrix <- triMatrix(phyloTree,level = "4")
fit_out <- fitSignatures(tri_matrix,
signaturesRef = signaturesRef,
min.mut.count = min.mut.count)
if(is.na(fit_out)){
message("The number of mutations in all branches are less than min.mut.count in patient ",
phyloTree@patientID)
return(NA)
}
cos_sim_matrix <- fit_out[[1]]$cosine.similarity
sig_level <- colnames(cos_sim_matrix)
signatures <- apply(cos_sim_matrix,1,function(x)names(which.max(x)))
if(any(grepl("Signature", signatures))){
signatures <- gsub('Signature ', '', signatures)
sig_level <- gsub('Signature ', '', sig_level)
}else if(any(grepl("SBS", signatures))){
signatures <- gsub('SBS', '', signatures)
sig_level <- gsub('SBS', '', sig_level)
}
treeData$Signature <- signatures[treeData$label]
treeData[treeData$Signature == ''| is.na(treeData$Signature)]$Signature <- "Unknown"
if("Unknown" %in% treeData$Signature){
sig_level <- c(sig_level,"Unknown")
}
treeData$Signature <- factor(treeData$Signature, levels = sig_level)
treeData <- treeData[order(treeData$Signature), ]
}else{
branch_type <- getBranchType(phyloTree)
types <- as.character(branch_type$Mutation_Type)
names(types) <- as.character(branch_type$Branch_ID)
treeData$Mutation_Type <- types[treeData$label]
# dplyr::mutate(Mutation_Type = types[label]) %>%
# as.data.table()
}
}
treeData <- treeData[treeData$distance!= 0|treeData$sample == rootLabel,]
return(treeData)
}
setPhyloTree <- function(tree, treeEdge, treeData, rootNode, mainTrunk, adjacentPoints, adjacentWs, adjacentAngles,
horizon = pi/2 , W = pi){
trunkPath <- rev(c(mainTrunk,rootNode))
if(length(trunkPath) > 0){
# subdat_list <- lapply(2:length(trunkPath), function(i){
# x1 <- treeData[treeData$end_num == trunkPath[i-1],]$x2
# y1 <- treeData[treeData$end_num == trunkPath[i-1],]$y2
# W <- W
# angle <- horizon
# distance <- treeEdge[treeEdge$endNum == trunkPath[i],]$length
# if(horizon == pi/2){
# x2 <- x1
# y2 <- y1 + distance
# }else{
# x2 <- x1 + distance*cos(angle)
# y2 <- y1 + distance*sin(angle)
# }
# subdat <- data.table::data.table('x1'= x1, 'y1' = y1,
# 'x2' = x2, 'y2' = y2,
# 'horizon' = horizon,'W' = W,
# 'distance' = distance, 'angle' = horizon,
# 'node' = trunkPath[i-1],'end_num' = trunkPath[i])
# return(subdat)
# })
#
# subdat <- dplyr::bind_rows(subdat_list)
# treeData <- rbind(treeData, subdat)
for(i in 2:length(trunkPath)){
x1 <- treeData[treeData$end_num == trunkPath[i-1],]$x2
y1 <- treeData[treeData$end_num == trunkPath[i-1],]$y2
W <- W
angle <- horizon
distance <- treeEdge[treeEdge$endNum == trunkPath[i],]$length
if(horizon == pi/2){
x2 <- x1
y2 <- y1 + distance
}else{
x2 <- x1 + distance*cos(angle)
y2 <- y1 + distance*sin(angle)
}
subdat <- data.table::data.table('x1'= x1, 'y1' = y1,
'x2' = x2, 'y2' = y2,
'horizon' = horizon,'W' = W,
'distance' = distance, 'angle' = horizon,
'node' = trunkPath[i-1],'end_num' = trunkPath[i])
treeData <- rbind(treeData, subdat)
}
}
if(length(adjacentPoints) > 0){
# subdat_list <- lapply(seq_len(length(adjacentPoints)), function(i){
# point <- adjacentPoints[i]
# startnode <- treeEdge[treeEdge$endNum == point, ]$node
# x1 <- treeData[treeData$end_num == startnode,]$x2
# y1 <- treeData[treeData$end_num == startnode,]$y2
# W <- adjacentWs[point]
# angle <- adjacentAngles[point]
# distance <- treeEdge[treeEdge$endNum == point, ]$length
# x2 <- x1 + distance*cos(angle)
# y2 <- y1 + distance*sin(angle)
# subdat <- data.table::data.table('x1'= x1, 'y1' = y1,
# 'x2' = x2, 'y2' = y2,
# 'horizon' = horizon,'W' = W,
# 'distance' = distance, 'angle' = angle,
# 'node' = startnode,'end_num' = point)
# return(subdat)
# })
# subdat <- dplyr::bind_rows(subdat_list)
# treeData <- rbind(treeData, subdat)
for(i in seq_len(length(adjacentPoints))){
point <- adjacentPoints[i]
startnode <- treeEdge[treeEdge$endNum == point, ]$node
x1 <- treeData[treeData$end_num == startnode,]$x2
y1 <- treeData[treeData$end_num == startnode,]$y2
W <- adjacentWs[point]
angle <- adjacentAngles[point]
distance <- treeEdge[treeEdge$endNum == point, ]$length
x2 <- x1 + distance*cos(angle)
y2 <- y1 + distance*sin(angle)
subdat <- data.table::data.table('x1'= x1, 'y1' = y1,
'x2' = x2, 'y2' = y2,
'horizon' = horizon,'W' = W,
'distance' = distance, 'angle' = angle,
'node' = startnode,'end_num' = point)
treeData <- rbind(treeData, subdat)
}
}
return(treeData)
}
calMainTrunk <- function(tree, treeEdge, rootNode, rootLabel = "NORMAL"){
Ntips <- Ntip(tree)
subtips <- treeEdge[!treeEdge$endNum > Ntips,]$endNum
distanceTable <- c(0)
# for(i in seq_len(length(subtips))){
# distanceTable <- cbind(distanceTable, -1)
# }
s <- vapply(seq_len(length(subtips)), function(i)-1, FUN.VALUE = numeric(1))
distanceTable <- c(distanceTable, s)
names(distanceTable) <- c("x",s)
distanceTable <- distanceTable[(-1)]
names(distanceTable) <- tree$tip.label[subtips]
t <- 1
path <- list()
while(t <= length(distanceTable)){
subPath <- c()
end <- which(tree$tip.label == names(distanceTable)[t])
name <- names(distanceTable)[t]
dis <- 0
while (TRUE) {
p <- which(treeEdge$endNum == end)
dis <- dis + treeEdge$length[p]
start <- treeEdge$node[p]
subPath <- append(subPath,end)
end <- start
if(start == rootNode){
break
}
}
path[[name]] <- subPath
distanceTable[name] <- dis
t <- t + 1
}
if(rootLabel == "NORMAL"){
distanceTable <- distanceTable[names(distanceTable) != rootLabel]
}
result <- path[[names(which.max(distanceTable))]]
}
getNodeAngle <- function(tree, treeEdge, mainTrunk,
numList, pointsList, rootNode,
horizon = pi/2, W = pi ,numNORMAL = NULL){
# print(pointsList)
if(!is.null(numNORMAL)){
treeEdge <- treeEdge[treeEdge$endNum != numNORMAL,]
}
nodes <- factor(rev(c(mainTrunk,rootNode)), levels = rev(c(mainTrunk,rootNode)))
adjacentPoints <- c()
adjacentPoints <- lapply(nodes, function(n){
p <- treeEdge[treeEdge$node == n,]$endNum %>% setdiff(nodes)
return(p)
}) %>% unlist()
# for(n in nodes){
# p <- treeEdge[treeEdge$node == n,]$endNum %>% setdiff(nodes)
# adjacentPoints <- append(adjacentPoints,p)
# }
left <- c()
right <- c()
leftList <- c()
rightList <- c()
i <- 1
# adresult <- lapply(seq_len(length(adjacentPoints)),function(i){
# point <- adjacentPoints[i]
# pos <- which(pointsList == point)
# if(i == 1){
# if(horizon < pi/2){
# return(list(l = numList[pos], ll = point, r = NULL, rl = NULL))
# }else{
# return(list(l = NULL, ll = NULL, r = numList[pos], rl = point))
# }
# }else{
# if(sum(right) <= sum(left)){
# return(list(l = NULL, ll = NULL, r = numList[pos], rl = point))
# }else{
# return(list(l = numList[pos], ll = point, r = NULL, rl = NULL))
# }
# }
# })
# left <- lapply(adresult, function(x)x$l) %>% unlist()
# leftList <- lapply(adresult, function(x)x$ll) %>% unlist()
# right <- lapply(adresult, function(x)x$r) %>% unlist()
# rightList <- lapply(adresult, function(x)x$rl) %>% unlist()
for(point in adjacentPoints){
if(i == 1){
if(horizon < pi/2){
leftList <- append(leftList, point)
pos <- which(pointsList == point)
left <- append(left,numList[pos])
}else{
rightList <- append(rightList, point)
pos <- which(pointsList == point)
right <- append(right,numList[pos])
}
}else{
if(sum(right) <= sum(left)){
rightList <- append(rightList, point)
pos <- which(pointsList == point)
right <- append(right,numList[pos])
}else{
leftList <- append(leftList, point)
pos <- which(pointsList == point)
left <- append(left,numList[pos])
}
}
i <- i+1
}
adjacentWs <- rep(0,length(tree$edge.length)+1)
adjacentAngles <- rep(0,length(tree$edge.length)+1)
if(length(rightList) > 0){
startr <- horizon - W/2
wrt <- W/2
totalR <- sum(right)
for(i in seq_len(length(rightList))){
if(rightList[i] <= length(tree$tip.label)){
n <- 1
}else{
n <- right[i]
}
wr <- wrt*(n/totalR)
adjacentWs[rightList[i]] <- wr
if(i == 1){
angler <- wr/2 + startr
c1 <- horizon != pi/2
c2 <- any(startr != 0,n == 1)
c3 <- any(startr != pi/2,n == 1)
if(all(c1,c2,c3)){
if(length(c(leftList,rightList))!=1){
angler <- startr
}
}
adjacentAngles[rightList[i]] <- angler
startr <- angler
}else{
if(n == 1 & horizon != pi/2){
startr <- startr - wr/2
}
angler <- wr/2 + startr + adjacentWs[rightList[i-1]]/2
adjacentAngles[rightList[i]] <- angler
startr <- angler
}
}
}
if(length(leftList) > 0){
wlt <- W/2
startl <- horizon + wlt
totalL <- sum(left)
for(i in seq_len(length(leftList))){
if(leftList[i] <= length(tree$tip.label)){
n <- 1
}else{
n <- left[i]
}
wl <- wlt*(n/totalL)
adjacentWs[leftList[i]] <- wl
if(i == 1){
anglel <- startl - wl/2
c1 <- horizon != pi/2
c2 <- any(startl != pi, n == 1)
c3 <- any(startl != pi/2, n == 1)
# angle extension
if(all(c1,c2,c3)){
if(length(c(leftList,rightList))!=1){
anglel <- startl
}
}
adjacentAngles[leftList[i]] <- anglel
startl <- anglel
}else{
if(n == 1& horizon != pi/2){
startl <- startl + wl/2
}
anglel <- startl- wl/2 - adjacentWs[leftList[i-1]]/2
adjacentAngles[leftList[i]] <- anglel
startl <- anglel
}
}
}
return(list(adjacentWs,adjacentAngles, adjacentPoints))
}
calChildNodeNum <- function(tree, treeEdge, mainTrunk, rootNode, ft = FALSE){
pointsList <- sort(unique(c(treeEdge$node,treeEdge$endNum)))
numList <- rep(1, length(pointsList))
lrnumList <- rep(0, length(pointsList))
# for(i in seq_len(length(pointsList))){
# numList[i] <- 1
# lrnumList[i] <- 0
# }
for(i in seq_len(length(pointsList))){
if(pointsList[i] == rootNode){
next
}
start <- pointsList[i]
end <- as.numeric(treeEdge[treeEdge$endNum == start,]$node)
while(TRUE){
pos <- which(pointsList == end)
numList[pos] <- numList[pos] + 1
if(end == rootNode){
break
}
start <- treeEdge[treeEdge$endNum == end, ]$endNum
end <- treeEdge[treeEdge$endNum == end, ]$node
}
# if(count > 1){
# numList[i] <- 0
# }
}
## samples are not in mainTrunk
nodeRange <- c(Ntip(tree)+1):(length(tree$edge.length)+1)
nodeOnTree <- c(mainTrunk,rootNode) %>% intersect(nodeRange)
nodeNoOnTree <- treeEdge[treeEdge$node %in% c(mainTrunk,rootNode), ]$endNum %>% setdiff(nodeOnTree) %>% intersect(nodeRange)
return(list(numList,pointsList, nodeNoOnTree, nodeOnTree))
}
getSigColors <- function(signatures){
signature_colors <- c("1" = "#E41A1C", "2" = "#377EB8","3" = "#7F0000",
"4" = "#35978f","5" = "#A6D854", "6" = "#FFD92F",
"7" = "#E5C494", "8" = "#8DD3C7","9" = "#6E016B",
"10" = "#BEBADA", "11" = "#e08214", "12" = "#80B1D3",
"13" = "#d6604d", "14" = "#ffff99","15" = "#FCCDE5",
"16" = "#FF6A5A","17" = "#BC80BD","18" = "#CCEBC5",
"19" = "#fb9a99", "20" = "#B6646A", "21" = "#9F994E",
"22" = "#7570B3" , "23" = "#c51b7d" , "24" = "#66A61E" ,
"25" = "#E6AB02" ,"26" = "#003c30", "27" = "#666666",
"28" = "#524762","29" = "#8926AC","30" = "#C42669",
"1A" = "#47B39F","1B" = "#B84C60", "R1" = "#B93183",
"R2" = "#B93183", "R3" = "#845AD1", "U1" = "#7BA52E",
"U2" = "#7A4287","7a" = "#663399", "7b" = "#AD8B5E",
"7c" = "#194682", "7d" = "#FF8931", "10a" = "#AD3D07",
"10b" = "#0D5DFF", "17a" = "#F2A200", "17b" = "#8C3B45",
"31" = "#009F9A", "32" = "#FF2B32", "33" = "#153E71",
"34" = "#FF8C2F", "35" ="#2F8CFF", "36" = "#62BA1F",
"37" = "#9745D6", "38" = "#3826CA", "39" = "#96C4BC",
"40" = "#982B9F", "41" = "#9C3C33", "42" = "#3EB15A",
"43" = "#EF1297", "44" = "#B1098C", "45" = "#883E93",
"46" = "#6ED66B", "47" ="#9F679F" , "48" = "#5B309F",
"49" = "#4E7A9F", "50" = "#9F6878", "51" = "#154E0F",
"52" = "#9D887D", "53" = "#143045", "54" = "#8F1B1D",
"55" = "#544C9F", "56" = "#7F4E10", "57" = "#429F74",
"58" = "#6E6A5F", "59" = "#235F9F", "60" = "#8A1D1B",
"84" = "#4B9F34", "85" = "#598F4D", "Unknown" = "black")
color_scale <- signature_colors[signatures]
return(color_scale)
}
labelBranch <- function(tree){
Root <- which(tree$tip.label == "NORMAL")
internalNodes <- sort(unique(tree$edge[,1]))
result <- list()
subnumList <- list()
subnumList <- lapply(seq_len(length(tree$edge.length)+1), function(i){
if(i > Ntip(tree)){
return(i)
}
else{
return(NA)
}
})
result <- rep(NA, (length(tree$edge.length)+1)) %>% as.list()
# for(i in seq_len(length(tree$edge.length)+1)){
# result[[i]] <- NA
# if(i > Ntip(tree)){
# subnumList[[i]] <- i
# }
# else{
# subnumList[[i]] <- NA
# }
# }
# print(result)
# print(subnumList)
end <- tree$edge[which(tree$edge[,2] == Root),1]
for(i in seq_len(length(tree$tip.label))){
row <- tree$edge[which(tree$edge[,2] == i), ]
if(i == Root){
next
}
while (TRUE) {
node <- row[1]
if(node == end){
result[[node]] <- append(result[[node]], i)
break
}
else{
result[[node]] <- append(result[[node]], i)
row <- tree$edge[which(tree$edge[,2] == node),]
}
}
}
rootNode <- tree$edge[which(tree$edge[,2] == Root),1]
for(i in (length(tree$tip.label)+1):(length(tree$edge.length)+1)){
if(i == rootNode){
next
}
row <- tree$edge[which(tree$edge[,2] == i), ]
while (TRUE) {
node <- row[1]
if(node == rootNode){
subnumList[[node]] <- append(subnumList[[node]], i)
break
}
else{
subnumList[[node]] <- append(subnumList[[node]], i)
row <- tree$edge[which(tree$edge[,2] == node),]
}
}
}
g <- lapply(result, function(x){return(paste(sort(tree$tip.label[x[-1]], decreasing = TRUE),collapse = "&"))})
return(list(g,subnumList))
}
Try the MesKit package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.