# function to plot a single [simple] PCA plot
single_pca_plot <- function(region = c('GLOBAL', 'AFR', 'EUR', 'AMR', 'EAS', 'CSA', 'OCE', 'MID'),
pc_scores,
which_pc1,
which_pc2){
# PCA plot specifics
pca_plot <- ggplot(aes=(text='s')) +
geom_point(size=2) +
theme_classic() +
labs(shape="Projects") +
theme(text = element_text(size=14, color='black'),
legend.title = element_text(size = 20),
legend.text = element_text(size = 18),
legend.key.size = unit(0.7, "lines")) +
guides(shape = guide_legend(override.aes = list(size = 0.85)),
color = guide_legend(override.aes = list(size = 5), nrow=3)) # nrow=1 for GLOBAL, MID and OCE; nrow=2 for the rest
# GLOBAL PCA
if (region == 'GLOBAL'){
region_color <- unique(pc_scores[c("genetic.region", "continent.colors")]) %>% pull(continent.colors, genetic.region)
pca_plot <- pca_plot +
geom_point(pc_scores, mapping = aes_string(x=which_pc1, y=which_pc2, color='genetic.region', shape='project'), alpha = 0.4, size=3) +
scale_color_manual(values=region_color, name='Regions') +
guides(shape="none")
x_lim = ggplot_build(pca_plot)$layout$panel_scales_x[[1]]$range$range
y_lim = ggplot_build(pca_plot)$layout$panel_scales_y[[1]]$range$range
plot_pca_density <- ggplot(pc_scores, aes_string(x=which_pc1, y=which_pc2)) +
geom_hex(bins=50) +
lims(x=x_lim, y=y_lim) +
scale_fill_gradientn(colours = rev(brewer.pal(5,'Spectral')), name = "Count") +
theme_classic() +
theme(legend.text = element_text(size = 12),
axis.text=element_text(size=10),
axis.title=element_text(size=12))
} else {
pop_count <- length(unique(pc_scores$population))
pop_color <- colorRampPalette(brewer.pal(7, "Set1"))(pop_count)
pca_plot <- pca_plot +
geom_point(pc_scores, mapping = aes_string(x=which_pc1, y=which_pc2, color='population', shape='project'), alpha = 0.4, size=3) +
scale_color_manual(values=pop_color, name='Populations') +
guides(shape="none")
x_lim = ggplot_build(pca_plot)$layout$panel_scales_x[[1]]$range$range
y_lim = ggplot_build(pca_plot)$layout$panel_scales_y[[1]]$range$range
plot_pca_density <- ggplot(pc_scores, aes_string(x=which_pc1, y=which_pc2)) +
geom_hex(bins=50) +
lims(x=x_lim, y=y_lim) +
scale_fill_gradientn(colours = rev(brewer.pal(5,'Spectral')), name = "Count") +
theme_classic() +
theme(legend.text = element_text(size = 12),
axis.text=element_text(size=11),
axis.title=element_text(size=13))
}
return(list(pca_plot, plot_pca_density))
}
plots_in_grid <- function(region, outlier_status, longitude_lim, latitude_lim){
region <- toupper(region)
outlier_status <- tolower(outlier_status)
region_pc_scores <- read.table(gzfile('./EAS_scores_without_outliers.txt.bgz'), header=T)
samples_to_remove <- c('LP6005441-DNA_A09', 'HGDP01371')
region_pc_scores <- region_pc_scores[!region_pc_scores$s %in% samples_to_remove, ]
plot_metadata <- read.delim('./gnomad_meta_for_plots.tsv', header=T, sep='\t')
pc_scores_meta <- region_pc_scores %>% left_join(plot_metadata) %>% arrange(population)
pc_scores_meta <- na.omit(pc_scores_meta)
data(wrld_simpl)
world <- fortify(wrld_simpl)
if (region == 'GLOBAL'){
color_vec <- unique(pc_scores_meta[c("genetic.region", "continent.colors")]) %>% pull(continent.colors, genetic.region)
p_map <- ggplot() +
geom_polygon(data = world, aes(long, lat, group=group), fill='lightyellow', color='lightgrey') +
geom_point(data = pc_scores_meta, aes(longitude, latitude, color=genetic.region, shape=project), size=3) +
coord_fixed(xlim = longitude_lim, ylim = latitude_lim) +
labs(x='Longitude', y='Latitude') +
theme_classic() +
scale_color_manual(values = color_vec) +
theme(panel.background = element_rect(fill = "lightblue"),
plot.background = element_rect(fill = "transparent", color = NA),
legend.position='none',
text = element_text(size=14),
axis.text = element_text(color='black'))
} else {
cont_metadata <- subset(pc_scores_meta, genetic.region == region)
color_vec <- colorRampPalette(brewer.pal(7, "Set1"))(length(unique(cont_metadata$population)))
p_map <- ggplot() +
geom_polygon(data = world, aes(long, lat, group=group), fill='lightyellow', color='lightgrey') +
geom_point(data = cont_metadata, aes(longitude, latitude, color=population, shape=project), size=3.5) +
coord_fixed(xlim = longitude_lim, ylim = latitude_lim) +
labs(x='Longitude', y='Latitude') +
theme_classic() +
scale_color_manual(values = color_vec) +
theme(panel.background = element_rect(fill = "lightblue"),
plot.background = element_rect(fill = "transparent", color = NA),
legend.position='none',
text = element_text(size=14),
axis.text = element_text(color='black'))
}
# plot PCs 1-4 with their respective density plots
pca_1_2 <- single_pca_plot(region, pc_scores_meta, 'PC1', 'PC2')
pca_3_4 <- single_pca_plot(region, pc_scores_meta, 'PC3', 'PC4')
### junk code ###
all_pca <- plot_grid(pca_1_2[[1]] + theme(legend.position="none"), pca_3_4[[1]] + theme(legend.position="none"), nrow=1) # put 2 plots together
legend_b <- get_legend(pca_1_2)
all_pca_add_legend <- plot_grid(all_pca, legend_b, rel_widths = c(4.8, 1)) # add legend
# arrange PCA plots together without a legend
pca_no_legend <- plot_grid(pca_1_2[[1]]+theme(legend.position="none"), pca_3_4[[1]]+theme(legend.position="none"))
# arrange density plots together
density_pca <- plot_grid(pca_1_2[[2]], pca_3_4[[2]])
# arrange PCA and density plots together
combined_plots <- plot_grid(pca_no_legend, density_pca, nrow=2)
# get the legend of the PCA plots and add it to the combined plots - the legends correspond to the colors found in the continental maps
legend = gtable_filter(ggplot_gtable(ggplot_build(pca_1_2[[1]] + theme(legend.position="bottom", legend.direction = "horizontal"))), "guide-box")
# create a ggplot grob from the combined plots (with legend) and the continental map so they can be arranged together using ggdraw() + draw_grob()
g_combined_plots <- ggplotGrob(combined_plots)
g_map <- ggplotGrob(p_map)
# arrange all plots together into one - adjust the coordinates and size as needed
combined_with_map <- ggdraw() +
draw_grob(g_map, x = 0, y = 0, width = 0.36, height = 1) +
draw_grob(g_combined_plots, x = 0.38, y = 0, width = 0.58, height = 1) +
draw_plot_label(label = c("A", "B", "C"), size = 15, x = c(0, 0.38, 0.38), y = c(1, 1, 0.5)) # add labels for plot reference
final_plot <- arrangeGrob(combined_with_map, top=textGrob(region, x = 0, hjust = 0, gp=gpar(fontsize=20))) # add region name as title to final plot
final_plot_w_legend <- grid.arrange(final_plot, legend, heights=c(1.1, .1), nrow = 2)
# save plots as png - when you run this function, this automatically writes out the plots
return(final_plot_w_legend)
}
EAS <- plots_in_grid('EAS', 'without_outliers', c(60,148), c(0,70))
ggsave('EAS_without_outliers_NEW.png', EAS, width=15, height=8) # 20 and 10
ggsave('EAS_without_outliers_NEW.png', EAS, width=15, height=9) # 20 and 10
