【StableTraits后续R实现可视化-单个性状】
library(ape)
library(phytools)
library(dplyr)
library(tidyr)
library(ggplot2)
library(ggtree)
data <- read.csv("summary.csv", sep = "\t", header = TRUE)
data_clean <- data %>%
separate(
col = "Trait.Model.Avg_DIC.Avg_PBIC",
into = c("Trait", "Model", "Avg_DIC", "Avg_PBIC"),
sep = ",",
convert = TRUE
)
head(data_clean)
best_dic <- data_clean %>%
group_by(Trait) %>%
arrange(Avg_DIC) %>%
slice(1) %>%
ungroup()
data_long <- data_clean %>%
pivot_longer(
cols = c("Avg_DIC", "Avg_PBIC"),
names_to = "Criterion",
values_to = "Value"
)
ggplot(data_long, aes(x = Trait, y = Value, color = Model, group = Model)) +
geom_line(linewidth = 1) +
geom_point(size = 3) +
facet_wrap(~ Criterion, scales = "free_y", ncol = 1) +
labs(
title = "DIC and PBIC Comparison Across Traits",
x = "Trait",
y = "Value",
color = "Model"
) +
theme_minimal() +
theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
scale_color_manual(values = c("red", "green", "blue"))
ggplot(data_long, aes(
x = reorder(Trait, Value),
y = Value,
fill = Model,
label = round(Value, 1)
)) +
geom_col(position = position_dodge(width = 0.8), width = 0.7) +
geom_text(
aes(y = ifelse(Value > 0, Value + 5, Value - 5)),
position = position_dodge(width = 0.8),
size = 3,
color = "black"
) +
facet_wrap(~ Criterion, scales = "free_x", nrow = 1) +
coord_flip() +
labs(
title = "DIC and PBIC Comparison by Trait and Model",
x = "Trait",
y = "Criterion Value (Lower is Better)",
fill = "Model"
) +
theme_minimal() +
theme(
axis.text.y = element_text(size = 10),
legend.position = "bottom",
panel.grid.major.y = element_blank()
) +
scale_fill_manual(
values = c("red", "green", "blue"),
labels = c("Brownian", "Brownian Strict", "Hetero")
)
tree <- read.tree("bodylength_brownian.tree")
anc_states <- read.table("bodylength_brownian.ancstates", header=TRUE, sep="\t")
p <- ggtree(tree) +
geom_tiplab(size=3) +
theme_tree2()
p <- p %<+% anc_states
p <- p +
geom_nodepoint(aes(color=Median), size=3) +
scale_color_gradient2(low="blue", mid="yellow", high="red",
midpoint=median(anc_states$Median)) +
labs(color="Ancestral Body Length")
print(p)
ggplot(rates_data, aes(x=Branch_midpoint_height, y=Rate)) +
geom_point(aes(color=Orig_len), size=3) +
geom_smooth(method="loess", color="blue") +
scale_color_gradient(low="yellow", high="red", name="Branch Length") +
labs(x="Branch Midpoint Height (Time)", y="Evolutionary Rate",
title="Evolutionary Rate Through Time") +
theme_minimal()
rates_data$Node <- sub("->.*", "", rates_data$Branch)
morphospace <- merge(anc_data, rates_data, by.x="Parameter", by.y="Node")
write.csv(morphospace,"bodylength_morphospace.csv")
morphospace <- read.csv("bodylength_morphospace.csv")
target_path <- c("n0", "n1", "n2", "n3", "n10", "Enas")
path_data <- morphospace %>%
filter(Parameter %in% target_path |
Branch %in% paste(target_path[-length(target_path)],
target_path[-1], sep="->"))
ggplot(path_data, aes(x = Branch_midpoint_height, y = Median)) +
geom_line(aes(group = 1), color = "darkgray", size = 1.5) +
geom_ribbon(aes(ymin = X95.CI_Low, ymax = X95.CI_High),
fill = "lightgray", alpha = 0.3) +
geom_point(aes(color = Rate, size = Rate)) +
geom_label(data = path_data[-1, ],
aes(label = sprintf("%.2f", Rate)),
vjust = 1.8, hjust = 0.5, size = 3.5,
fill = "white", alpha = 0.8) +
geom_text(aes(label = Parameter),
vjust = -1.5, hjust = 0.5, size = 4) +
scale_color_gradientn(colors = c("blue", "green", "yellow", "red"),
name = "Evolutionary Rate") +
scale_size_continuous(range = c(3, 8), guide = "none") +
labs(x = "Time (Branch Midpoint Height)",
y = "Median Body Length",
title = "Detailed View of Evolutionary Changes Along n0-n1-n2-n3-n10-Enas Path",
subtitle = "Branch-specific evolutionary rates shown below each segment") +
theme_minimal() +
theme(plot.title = element_text(hjust = 0.5, face = "bold"),
plot.subtitle = element_text(hjust = 0.5),
legend.position = "bottom")
tree <- read.tree("bodylength_brownian.tree")
morphospace <- read.csv("bodylength_morphospace.csv")
node_ages <- node.depth.edgelength(tree)
max_age <- max(node_ages)
x_position <- node_ages - max_age
tree_data <- ggtree(tree)$data %>%
left_join(morphospace, by = c("label" = "Parameter")) %>%
mutate(
x_position = x_position[node],
y_bodylength = ifelse(isTip, Brownian, NA)
)
internal_nodes <- tree_data %>%
filter(!isTip) %>%
select(node, label) %>%
left_join(
morphospace %>%
filter(grepl("^n", Parameter)) %>%
select(Parameter, Brownian),
by = c("label" = "Parameter")
)
tree_data <- tree_data %>%
left_join(internal_nodes %>% select(node, Brownian_int = Brownian), by = "node") %>%
mutate(
y_bodylength = coalesce(y_bodylength, Brownian_int)
)
p_right_aligned <- ggplot(tree_data, aes(x = x_position, y = y_bodylength)) +
geom_segment(
aes(x = x_position[parent], xend = x_position,
y = y_bodylength[parent], yend = y_bodylength,
color = Rate),
linewidth = 1.2, alpha = 0.8
) +
geom_point(
aes(size = coalesce(Brownian, Brownian_int), color = Rate),
alpha = 0.8
) +
geom_text(
aes(label = label),
hjust = -0.2, vjust = 0.5, size = 3
) +
labs(
x = "Time Before Present (Myr)",
y = "Body Length (mm)",
title = "Phylogeny in Time-Morphospace (Bodylength)"
) +
scale_color_gradientn(
colors = c("blue", "green", "yellow", "red"),
name = "Evolutionary Rate",
na.value = "gray50"
) +
scale_size_continuous(name = "Body Length") +
scale_x_continuous(labels = function(x) abs(x)) +
theme_minimal() +
theme(legend.position = "right")
print(p_right_aligned)