Code
suppressPackageStartupMessages({
library(purrr)
library(dplyr)
library(tibble)
library(tidyr)
library(ggplot2)
library(colorblindr)
library(patchwork)
library(ggbeeswarm)
library(ggnewscale)
})
set.seed(42)Visualisation: Common DNA sequence variation influences epigenetic aging in African Populations
Visualisation
Richard J. Acton 
manuscript doi: 10.1101/2024.08.26.608843
Overview
This paper (Meeks et al. 2024) has quite a lot of data which lends itself well to ‘small-multiples’ plots where different facets of the same dataset are represented in the same or similar ways. Small multiples plots are a good reason to have multiple panels in a single figure unlike multi-panel figures where unrelated sub-panels are placed in the same figure, as is commonly done when the ‘natural’ number of figures is more than the limit of figures allowed by a journal.
As is common the data underlying the figures ‘as plotted’ is not provided as a supplement, though this is a pre-print and it might accompany the final paper. Consequently we have not done much re-visualisation of these figures
We have some minor critiques of the data visualisation to share, focused on improving ease of readability the graphs.
I will quickly restate my perennial complaint about bioRxiv’s layout and note that figures should be adjacent to the text referencing them and not placed all together at the end of the manuscript.
Figure 1
Given that the sample size here appears to range from 35-130 I would consider a beeswarm plot preferable to a violin plot for these figures. Plotting the raw points is helpful for illustrating the relative sample sizes of your groups which when comparing the amount of error can be quite relevant to bear in mind. Violin plots are best suited to applications where it is impractical to plot the individual points due their shear number. Much like a histogram their shape can end up being a little misleading especially for sparse data, and if the bin size/bandwidth is not optimal. Violin plotting tools tend to obfuscate the parameters used in the kernel density estimation function used to smooth the curve so it is not simple to check their suitability for your data.
There is quite a lot of axis label duplication in this plot matrix that may not be necessary. “Age-Adjusted prediction error” only really needs to appear once, at most twice to indicate the Y axis values. The presence of multiple redundant copies of this label make the plot noisier than necessary. The long axis labels mean the reader is inclined to read each or them and scan for differences as in similar plot layouts these may be different values, removing them removes the possibility that they might be different in a way the reader needs to pay attention to. Similarly the group names on the x axis tick labels likely only need to appear on the bottom row of plots and not the top, especially as the group variable is already encoded by the colour redundantly with its position on the x axis.
The data ranges on the plot differ slightly in the y axis but the primary point of comparison is different in error between groups not clocks so this is reasonable, though given the relatively minor difference in the data range it may have been possible to use the same range across plots to further minimise differences between the plot that need to be conveyed to the reader.
Note the order in which the groups appear: 1) Himba, 2) ‡Khomani San, 3) Baka, 4) Hispanic/Latino 5) European. This is consistent with some subsequent figures, but not with others, it would be preferable to preserve the order throughout. The colours used for these group are consistent throughout, though it may have been possible to make more use of these colours to indicate the groups and make it easier to identify at a glance which group’s data was being plotted, we’ll come back to this.
Code
comparison_names <- c(
"Baka vs ‡Khomani San", "Himba vs ‡Khomani San", "Himba vs Baka"
)
properties <- c("Effect Size", "Heritability", "Weight")
orig_group_colours <- c(
"Himba" = "#CD7C90", "‡Khomani San" = "#B583BC", "Baka" = "#CFA373",
"Hispanic / Latino" = "#93A1DE", "European" = "#99C48C"
)
group_colours <- c(
"Himba" = "#FBE523", "‡Khomani San" = "#1F8F89", "Baka" = "#75207F",
"Hispanic / Latino" = "#FB9064", "European" = "#271057"
)
estimators <- c(
"Horvath2013", "Hannum et al. 2013", "Skin & Blood Clock", "Zhang et al. 2019",
"PhenoAge", "FitAge", "GrimAge2 (true age)", "GrimAge2 (predicted age)"
)
f1_dummy_data <- tibble(
estimator = factor(estimators, levels = estimators, ordered = TRUE),
group = map(1:8, ~factor(
names(group_colours), levels = names(group_colours), ordered = TRUE
)),
mean_error = list(
c( 0, -1, -1, -1, 0),
c( -3, -1, 0, 1, 0),
c( -7, 0, 1, 0, 1),
c( -4, -2, 0, 1, 1),
c(-11, 0, 3, -1, 1),
c( -2, 4, 3, -1, -1),
c( 1, 7, 4, -1, -3),
c( -3, 6, 5, -1,- 2)
),
n = map(1:8, ~c(49, 36, 45, 69, 130))
) %>%
unnest(cols = c(group, mean_error, n)) %>%
mutate(error = map2(mean_error, n, \(m, n) { rnorm(n = n, mean = m) })) %>%
unnest(cols = error)
f1_plotf <- \(data, group_colours = orig_group_colours) {
ggplot(data, aes(group, error)) +
geom_quasirandom(aes(colour = group), show.legend = FALSE, size = 0.5) +
geom_boxplot(
aes(group = group), outlier.shape = NA, width = 0.2, fill = NA, colour = "darkgrey"
) +
scale_colour_manual(values = group_colours) +
geom_hline(yintercept = 0) +
facet_wrap(~estimator, nrow = 2) +
labs(y = "Age-Adjusted Preediction Error", x = "") +
theme_bw() +
theme(axis.text.x = element_text(angle = 30, vjust = 1, hjust = 1))
}
f1_dummy_data %>% f1_plotf(group_colours = group_colours) 
The colour pallette used appears fairly accessible with the three main groups being distinguishable in all but the completely de-saturated case.
Code
f1_dummy_data %>%
filter(estimator == "Horvath2013") %>%
f1_plotf(group_colours = orig_group_colours) %>%
cvd_grid()
However we can probably do better using some colours spaced out over a couple of the virisis palettes.
Code
f1_dummy_data %>%
filter(estimator == "Horvath2013") %>%
f1_plotf(group_colours = group_colours) %>%
cvd_grid()
Now even the desaturated case has quite distinguishable shades of grey.
Figure 2
In A-C the background might be better represented by a density colour scale due to the high number of overplotted points. Given that this is not the main point of the plot, which is primarily conveyed by the points highlighted in red. It might make sense to do the density colour scale in a simple greyscale with a relatively narrow colour range to keep the background points relatively flat to provide a backdrop to the highlighted foreground points. I’ve tried out a version with a greyscale density plot for the non-significant sites and I think it still does a decent job of conveying the distribution of these points. This also has the advantage of having many fewer objects in the plot so it is smaller and easier to render when saved in a vector format like svg or pdf. I have also coloured the axes to correspond to the group the value for which are being plotted on that axis. This lets us continue to use the colour key for the groups that we established in figure 1 to provide a convenient perceptual que for which groups are being compared.
Note that the order of the comparisons between the groups was: 1) Baka vs ‡Khomani San, 2) Himba vs ‡Khomani San, 3) Himba vs Baka of the form (y-axis vs x-axis) This is not consistent with subsequent orders in which group comparisions where represented.
Whilst fine to do so it is not strictly necessary to have a different colour for every chromosome in a Manhattan plot, they can just alternate. In a context where you might want to use some similar colours to mean something other than chromosome number elsewhere, this approach may be better so as not to overwhelm the reader to many colours at the same time.
The p-value of the horizontal red dotted line plotted presumably to indicate the significance threshold is not indicated on the figure or its legend only mentioned in the text.
Code
# Dummy data generation functions
gen_random_correlated <- function(n = 100, sd = 0.3) {
r_norm <- rnorm(n, sd = sd)
r_unif <- c(runif(n/2, max = sd), runif(n/2, max = sd) * -1)
c(r_norm, r_norm + r_unif)
}
gen_random_norm_uncor <- function(n = 100, sd = 5) {
c(rnorm(n, sd = sd), rnorm(n, sd = sd))
}
f2_dummy_data <- tibble(
group = c(rep("sig", 2e3), rep("nsig", 2e5)),
axis = c(rep(c("x", "y"), each = 1e3), rep(c("x", "y"), each = 1e5)),
value = c(gen_random_correlated(1e3, 300), gen_random_norm_uncor(1e5, 2e3))
) %>%
pivot_wider(names_from = axis, values_from = value, values_fn = list) %>%
unnest(cols = everything())
f2_dummy_plot <- f2_dummy_data %>% {
ggplot(., aes(x, y)) +
geom_density2d_filled(
data = . %>% filter(group == "nsig"), show.legend = FALSE, alpha = 0.7
) +
scale_fill_grey(start = 1, end = 0.2) +
geom_point(data = . %>% filter(group == "sig"), color = "red", size = 0.5) +
geom_abline(colour = "red", slope = 1) +
labs(
title = "A", y = "Effect Size in Baka", x = "Effect Size in ‡Khomani San"
) +
coord_equal() +
theme_bw() +
theme(
axis.line.y.left = element_line(
colour = group_colours["Baka"], size = 1
),
axis.line.x.bottom = element_line(
colour = group_colours["‡Khomani San"], size = 1
)
)
}
f2_dummy_plot 
Figure 3
There is potential here for a clearer and better aligned layout of this faceted grid. Using the row titles: Effect size, Heritability, Weight means that these do not need to be repeated in the x axis labels. Using column titles consisting of the comparisons listed in the section above on figure 2 would make it clear that the groups being compared are consistent in each column. Also keeping the axes on which pairwise comparisons are made consistent between these two figures instead of changing the order of comparisons, and which axis values are on, would improve readability as this plot would conform to the expectation established in the previous one. We can re-enforce this with the use of the coloured axes as we did in figure 2. Using row titles would permit you to drop that variable from the axis names and simply include the group name. This is especially helpful for improving the ability to quickly scan the plot and know which groups are being compared given that order of the group and the statistic are in inconsistent places in the names. Specifically: “Effect size in <Group>”, and “Weight in <group>” start with the statistic whereas “<Group> Heritability” ends with it. This make the labels more difficult to process at a glance as you at not looking at a consistent place in the labal.
Note Comparison Order: 1) Baka vs Himba, 2) ‡Khomani San vs Baka, 3) ‡Khomani San vs Himba
Code
# Dummy data generation functions
gen_random_unif_uncor <- function(n = 100) {
c(runif(n, min = 0.3), runif(n, min = 0.3))
}
# Correlation plot function
cor_plot <- function(x) {
group1 <- sub(".* vs (.*)", "\\1", x$comparison[1])
group2 <- sub("(.*) vs .*", "\\1", x$comparison[1])
x %>%
ggplot2::ggplot(aes(x, y)) +
ggplot2::geom_point() +
ggplot2::geom_abline(
slope = 1, linewidth = 0.5, colour = "red", linetype = "dashed"
) +
# ggplot2::geom_smooth(
# method = "lm", se = FALSE, formula = "y ~ x",
# linewidth = 0.5, colour = "red", linetype = "dashed"
# ) +
labs(x = group1, y = group2) +
switch (as.character(x$property[1]),
"Effect Size" = lims(x = c(-1, 1), y = c(-1, 1)),
"Heritability" = lims(x = c(0, 1), y = c(0, 1)),
"Weight" = lims(x = c(-20, 20), y = c(-20, 20))
) +
theme_bw() +
theme(
axis.line.x.bottom = element_line(
colour = group_colours[group1], size = 1
),
axis.line.y.left = element_line(
colour = group_colours[group2], size = 1
)
) +
facet_grid(property ~ comparison)
}
# Generating Dummy Data and individual plots
f3_dummy_data <- tibble(
value = c(
gen_random_correlated(),
gen_random_correlated(),
gen_random_correlated(),
gen_random_unif_uncor(),
gen_random_unif_uncor(),
gen_random_unif_uncor(),
gen_random_norm_uncor(),
gen_random_norm_uncor(),
gen_random_norm_uncor()
),
group = factor(
x = rep(c(rep("x", 100), rep("y", 100)), 9),
levels = c("x", "y")
),
property = factor(c(rep(properties, each = 600)), levels = properties),
comparison = factor(
rep(rep(comparison_names, each = 200), 3), levels = comparison_names
)
) %>%
pivot_wider(names_from = group, values_from = value, values_fn = list) %>%
unnest(cols = everything()) %>%
group_by(comparison, property) %>%
nest() %>%
mutate(
data = map(data, ~{
.x %>%
mutate(property = property, comparison = comparison) # %>%
# rename_with(~map_chr(.x, ~{
# if(.x == "x") {
# sub("(.*) vs .*", "\\1", comparison[1])
# } else if (.x == "y"){
# sub(".* vs (.*)", "\\1", comparison[1])
# } else {.x}
# }))
}),
plot = map(data, cor_plot)
)
# Composing plots
no_strip_theme <- theme(
strip.background = element_blank(), strip.text = element_blank()
)
x_strip_theme <- theme(
strip.background.x = element_blank(), strip.text.x = element_blank()
)
y_strip_theme <- theme(
strip.background.y = element_blank(), strip.text.y = element_blank()
)
patchwork <- (
(f3_dummy_data$plot[[1]] + y_strip_theme) +
(f3_dummy_data$plot[[2]] + y_strip_theme) +
(f3_dummy_data$plot[[3]])
) / (
(f3_dummy_data$plot[[4]] + no_strip_theme) +
(f3_dummy_data$plot[[5]] + no_strip_theme) +
(f3_dummy_data$plot[[6]] + x_strip_theme)
) / (
(f3_dummy_data$plot[[7]] + no_strip_theme) +
(f3_dummy_data$plot[[8]] + no_strip_theme) +
(f3_dummy_data$plot[[9]] + x_strip_theme)
) + plot_annotation(tag_levels = "A")
patchwork 
Figure 4
Given that panels A-C on the first row and panels E-G on the second correspond to each other they should ideally line up in columns A/E etc. This would leave a gap under the meta analysis panel D but this could readily be filled with the colour keys from both of these sets of plots which are currently missing from the figure. It might not be relevant to include the specifics of which genotypes are denoted by the colours for each CpG site but it would be good to have a visual key showing that the different colours correspond to different genotypes, not just a note in the caption.
It is also worth noting that orange points are used in both plots types here but they have different meanings in each plot so ideally a different colour should have been chosen to represent either the permutation results or genotype.
The emphasis on the red point in the first row indicating the CpG for which the genotype / age / methylation data is plotted below could perhaps be a little clearer, possibly with the use of a different shape. See the re-works of figure 5 A-C below for an attempt to add this kind of emphasis for those plots.
Figure 5
Similarly to figure 5 the first row of figures and the second are not aligned despite the second depicting genotype / age / methlation data for the points highlighted in the first. This alignment can likewise be achieved by adjusting the key / legend.
The significance of the red points, the highlighted points, and the genotype are missing from legend, I have added the genotype in the re-work but not the other two. For the highlighted points I’ve used a slighly transleucent group colour background of an oversized diamond point with crossed diamond lines overlayed in black to ensure both that the point is more visible and that its position is clear. Unfortunately it’s not very straightforward to programatically generate a legend for the highlighted meQTLs though it would be fairly simple to export this plot to an .svg and do a little manual editing in Inkscape to make one.
I’ve also attempted a version of panels A-C with a density plot for the 1000 Genomes data instead of individual points.
I’ve also added coloured backgrounds to the column title strips that correspond to the colours I’ve been using to denote the different populations elsewhere to clearly convey the mapping of the column to the group.
Code
linscale <- \(x){ (x - min(x)) / (max(x) - min(x)) }
r_beta <- linscale(rbeta(3000, 1, 5))
r_comb <- linscale(r_beta + abs(rnorm(3000)))
f5abc_dummy_data <- tibble(
x = r_comb, y = r_beta, group = factor(
rep(names(group_colours)[1:3], 1000), levels = names(group_colours)[1:3],
ordered = TRUE
),
meQTL_CpGs = c(
rep("highlighted", 3), rep("age associated", 900),
rep("not age associated", 2097)
),
row = "meQTL"
)
f5abc_plotf <- \(x, fill_guide = TRUE) {
x %>% {
ggplot(., aes(x, y)) +
geom_density2d_filled(
data = . %>% filter(meQTL_CpGs == "not age associated"),
alpha = 0.6, show.legend = fill_guide
) +
scale_fill_viridis_d(option = "E") +
guides(colour = "none") +
geom_abline(
slope = 1, linewidth = 0.6, colour = "white", linetype = "dashed"
) +
geom_point(
data = . %>% filter(meQTL_CpGs == "age associated") ,
size = 0.6, colour = "red"
) +
labs(y = "European Frequency", x = "Frequency", fill = "Density") +
ggnewscale::new_scale("fill") +
guides(fill = "none") +
scale_fill_manual(values = group_colours) +
scale_colour_manual(values = group_colours) +
geom_point(
data = . %>% filter(meQTL_CpGs == "highlighted"),
aes(fill = group, colour = group), size = 4, shape = 23, alpha = 0.8
) +
geom_point(
data = . %>% filter(meQTL_CpGs == "highlighted"),
size = 4, shape = 9, colour = "black"
) +
lims(x = 0:1, y = 0:1) +
theme_bw() +
theme(
strip.text = element_text(face = "bold", colour = "darkgrey", size = 12)
) +
facet_grid(row ~ group)
}
}
f5a_plot <- f5abc_dummy_data %>%
filter(group == "Himba") %>%
f5abc_plotf(fill_guide = FALSE) +
labs(x = "") +
y_strip_theme +
theme(strip.background.x = element_rect(fill = group_colours["Himba"]))
f5b_plot <- f5abc_dummy_data %>%
filter(group == "‡Khomani San") %>%
f5abc_plotf(fill_guide = FALSE) +
labs(y = "") +
y_strip_theme +
theme(strip.background.x = element_rect(
fill = group_colours["‡Khomani San"]
))
f5c_plot <- f5abc_dummy_data %>%
filter(group == "Baka") %>%
f5abc_plotf() +
labs(y = "", x = "") +
theme(
strip.background.x = element_rect(fill = group_colours["Baka"]),
strip.background.y = element_rect(fill = "white")
)
f5efg_dummy_data <- tibble(
Methylation = c(runif(300, 0.6, 1), runif(30, 0.4, 0.6), runif(3, 0.1, 0.3)),
Age = sample(5:95, 333, replace = TRUE),
Genotype = factor(
c(rep("A", 300), rep("B", 30), rep("C", 3)), levels = LETTERS[1:3],
ordered = TRUE
),
CpG = c(
rep("cg25256723", 100), rep("cg05163071", 100), rep("cg21860825", 100),
rep("cg25256723", 10), rep("cg05163071", 10), rep("cg21860825", 10),
rep("cg25256723", 1), rep("cg05163071", 1), rep("cg21860825", 1)
),
row = "CpG"
)
f5efg_plotf <- \(x) {
x %>%
ggplot(aes(Age, Methylation)) +
geom_point(aes(colour = Genotype)) +
scale_colour_manual(values = c(
"A" = "#E49D00", "B" = "#56B4E8", "C" = "#004848"
)) +
theme_bw() +
theme(
strip.text = element_text(face = "bold", colour = "darkgrey", size = 12)
) +
lims(x = c(5, 95), y = 0:1) +
facet_grid(row ~ CpG)
}
f5e_plot <- f5efg_dummy_data %>%
filter(CpG == "cg25256723") %>%
f5efg_plotf() +
labs(x = "") +
guides(colour = "none") +
y_strip_theme +
theme(strip.background.x = element_rect(fill = group_colours["Himba"]))
f5f_plot <- f5efg_dummy_data %>%
filter(CpG == "cg05163071") %>%
f5efg_plotf() +
labs(y = "") +
guides(colour = "none") +
y_strip_theme +
theme(strip.background.x = element_rect(
fill = group_colours["‡Khomani San"]
))
f5g_plot <- f5efg_dummy_data %>%
filter(CpG == "cg21860825") %>%
f5efg_plotf() +
labs(y = "", x = "") +
theme(
strip.background.x = element_rect(fill = group_colours["Baka"]),
strip.background.y = element_rect(fill = "white")
)
f5plot <-
(f5a_plot + f5b_plot + f5c_plot) /
(f5e_plot + f5f_plot + f5g_plot) +
plot_annotation(tag_levels = "A")
f5plot
Code
f5plot %>% cvd_grid()
Figure 6
I like the use of dark / light shades to convey heritable vs non-heritable predictors it allows the dual encoding of group with colour and the x axis, although a larger difference in lightness may have been desirable here for greater clarity
The order of African, European, Hispanic / Latino is consistent within figure 6 but inconsistent with figure 1 though the colours are consistent.
Figure 7
There is a lot of axis label duplication in this figure, EAS and the metric name need only be present once per row.
The group colours are consistent with consistent between figures 1 and 7 as is the order, though the order is not consistent with other figures such as figure 5.
There is no need to tell the reader it is a scatter plot in the figure legend, they can see that. This might be useful information in alt text but the legend is not alt text.
Session Info
Code
sessionInfo()R version 4.3.1 (2023-06-16)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 22.04.4 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.20.so; LAPACK version 3.10.0
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
time zone: Etc/UTC
tzcode source: system (glibc)
attached base packages:
[1] stats graphics grDevices datasets utils methods base
other attached packages:
[1] ggnewscale_0.4.7 ggbeeswarm_0.6.0 patchwork_1.1.1 colorblindr_0.1.0
[5] colorspace_2.0-3 ggplot2_3.3.6 tidyr_1.2.0 tibble_3.2.1
[9] dplyr_1.1.4 purrr_0.3.4
loaded via a namespace (and not attached):
[1] gtable_0.3.0 jsonlite_1.8.0 compiler_4.3.1 renv_0.15.5
[5] tidyselect_1.2.0 stringr_1.4.0 scales_1.3.0 yaml_2.3.5
[9] fastmap_1.1.0 R6_2.5.1 labeling_0.4.2 generics_0.1.2
[13] isoband_0.2.5 knitr_1.39 MASS_7.3-60 htmlwidgets_1.6.4
[17] munsell_0.5.0 pillar_1.9.0 rlang_1.1.0 utf8_1.2.2
[21] stringi_1.7.6 xfun_0.38 viridisLite_0.4.0 cli_3.6.2
[25] withr_2.5.0 magrittr_2.0.3 digest_0.6.29 grid_4.3.1
[29] cowplot_1.1.1 beeswarm_0.4.0 lifecycle_1.0.3 vipor_0.4.5
[33] vctrs_0.6.5 evaluate_0.15 glue_1.6.2 farver_2.1.0
[37] fansi_1.0.3 rmarkdown_2.14 ellipsis_0.3.2 tools_4.3.1
[41] pkgconfig_2.0.3 htmltools_0.5.7 See also the companion piece to this post on openness and reproducibility in this manuscript.
References
Meeks, Gillian L., Brooke Scelza, Hana M. Asnake, Sean Prall, Etienne Patin, Alain Froment, Maud Fagny, Lluis Quintana-Murci, Brenna M. Henn, and Shyamalika Gopalan. 2024. “Common DNA Sequence Variation Influences Epigenetic Aging in African Populations.” bioRxiv: The Preprint Server for Biology, August, 2024.08.26.608843. https://doi.org/10.1101/2024.08.26.608843.