diff --git a/05-data_visualization/amel_bterall_samples_heatmap.Rmd b/05-data_visualization/amel_bterall_samples_heatmap.Rmd
deleted file mode 100644
index 249ccb4..0000000
--- a/05-data_visualization/amel_bterall_samples_heatmap.Rmd
+++ /dev/null
@@ -1,152 +0,0 @@
----
-title: "Differential expression analysis of nAChR subunits in _Apis mellifera_ and _Bombus terrestris_"
-author: "Federico Lopez, Alicja Witwicka"
-date: '`r Sys.Date()`'
-output:
-  github_document:
-    toc: yes
-  pdf_document: 
-    fig_caption: yes
-    toc: yes
-  html_document:
-    toc: yes
-number_sections: true
-
----
-
-```{r setup, echo = FALSE}
-knitr::opts_chunk$set(
-  echo = TRUE,
-  message = FALSE,
-  warning = FALSE,
-  cache.lazy = FALSE,
-  include = FALSE,
-  out.height = "\textheight",
-  out.width = "\textwidth"
-)
-```
-
-# Load libraries
-```{r}
-load_pkgs <- function(pkg) {
-  sapply(pkg, require, character.only = TRUE)
-}
-
-cran_pkgs <- c(
-  "tidyverse", "pheatmap", "styler", "here", "patchwork", "ggplot2"
-)
-load_pkgs(cran_pkgs)
-```
-
-# Plot heatmap for Apis mellifera samples
-```{r}
-here::here()
-
-subunit_levels <- c(
-  "alpha9", "alpha8", "alpha7", "alpha6", "alpha5", "alpha4", "alpha3", "alpha2", "alpha1", "beta2", "beta1"
-  )
-
-amel_heatmap <- read.csv(
-  here::here(
-    "results",
-    "amel_vst_counts.csv"
-    ),
-  header = TRUE,
-  stringsAsFactors = TRUE
-)
-
-levels(as.factor(amel_heatmap$group))
-
-sample_levels_amel <- c("Queen brain Liberti et al. 2019", "Queen brain Manfredini et al. 2015",   
-                   "Worker brain Christen et al. 2018" , "Worker brain Jasper et al. 2015", 
-                   "Thoracic ganglia Jasper et al. 2015", "Antenna Jasper et al. 2015", 
-                   "Hypopharyngeal gland Jasper et al. 2015", "Malpighian tubule Jasper et al. 2015", 
-                   "Mandibular gland Jasper et al. 2015", "Midgut Jasper et al. 2015", 
-                   "Muscle Jasper et al. 2015", "Nasonov gland Jasper et al. 2015")
-
-gg_amel_heatmap <- ggplot2::ggplot(amel_heatmap, aes(
-  x = sample, y = factor(subunit, levels = subunit_levels), fill = vst_count
-)) +
-  geom_tile(color = "white") +
-  theme_minimal() +
-  theme(
-    plot.title = element_text(hjust = 0.5, size = 14, face = "plain"),
-    plot.subtitle = element_text(hjust = 0.5, size = 14, face = "plain"),
-    legend.title = element_blank(),
-    legend.text = element_text(size = 10),
-    axis.title = element_text(size = 12),
-    axis.text = element_text(size = 8),
-    #axis.text.x = element_text(angle = 70, hjust = 1, vjust = 0.5),
-    axis.text.x = element_blank(),
-    panel.spacing = unit(2.5, "mm"),
-    strip.text = element_text(size = 8, angle = 0, hjust = 0, vjust = 0.5)
-  ) +
-  labs(
-    title = "Heatmap of transformed counts from honey bee castes and tissues",
-    x = "Sample",
-    y = expression(paste("Log"[2], " (VST) transformed count"))
-  ) +
-  scale_fill_gradientn(colours = c("#004d99", "#3B9AB2", "#78B7C5", "#EBCC2A")) +
-  facet_grid(~ factor(group, , levels = sample_levels_amel), scales = "free", space = "free") 
-
-```
-
-# Plot heatmap for Bombus terrestris samples
-```{r}
-bter_heatmap <- read.csv(
-  here::here(
-    "results",
-    "bter_vst_counts.csv"
-    ),
-  header = TRUE,
-  stringsAsFactors = TRUE
-)
-
-levels(as.factor(bter_heatmap$group))
-
-sample_levels_bter <- c("Queen head Colgan et al. 2019", "Worker head Colgan et al. 2019",
-                   "Queen whole body Harrison et al. 2015", "Worker whole body Harrison et al. 2015",
-                   "Worker larva Harrison et al. 2015", "Worker pupa Harrison et al. 2015")
-
-gg_bter_heatmap <- ggplot2::ggplot(bter_heatmap, aes(
-  x = sample, y = factor(subunit, levels = subunit_levels), fill = vst_count
-)) +
-  geom_tile(color = "white") +
-  theme_minimal() +
-  theme(
-    plot.title = element_text(hjust = 0.5, size = 14, face = "plain"),
-    plot.subtitle = element_text(hjust = 0.5, size = 14, face = "plain"),
-    legend.title = element_blank(),
-    legend.text = element_text(size = 10),
-    axis.title = element_text(size = 12),
-    axis.text = element_text(size = 8),
-    #axis.text.x = element_text(angle = 70, hjust = 1, vjust = 0.5),
-    axis.text.x = element_blank(),
-    panel.spacing = unit(2.5, "mm"),
-    strip.text = element_text(size = 8, angle = 0, hjust = 0, vjust = 0.5)
-  ) +
-  labs(
-    title = "Heatmap of transformed counts from bumble bee castes and developmental stages",
-    x = "Sample",
-    y = expression(paste("Log"[2], " (VST) transformed count"))
-  ) +
-  scale_fill_gradientn(colours = c("#004d99", "#3B9AB2", "#78B7C5", "#EBCC2A")) +
-  facet_grid(~ factor(group, , levels = sample_levels_bter), scales = "free", space = "free") 
-
-```
-
-# Combine plots
-```{r}
-heatmap <- gg_amel_heatmap / gg_bter_heatmap
-
-# Save
-ggsave(heatmap,
-       filename = here::here(
-         "results", "2020-07-28-gene_expression",
-         "figure_three.pdf"
-       ),
-       width = 9,
-       height = 10,
-       units = "in"
-)
-```