Ensure same OSD number and library kit order for each gastroc QC plot

Manuscript_Code/NBISC100_project/gastrocnemius/DGE_assessment/DGE_comparison_gastroc_clayton.py

@@ -0,0 +1,72 @@
+import pandas as pd
+from matplotlib_venn import venn3
+import matplotlib.pyplot as plt
+
+from supervenn import supervenn  # optional, but useful for comparing > 4 datasets
+
+# skip first row since it just lists dataset name
+# read in each sheet as separate df
+# some gene symbols are nan, so drop them
+OSD_101_gastroc = pd.read_excel("DEGs OSD101_419_401.xlsx", sheet_name="OSD-101 DGE", skiprows=1).dropna()
+OSD_401_gastroc = pd.read_excel("DEGs OSD101_419_401.xlsx", sheet_name="OSD-401 DGE", skiprows=1).dropna()
+OSD_419_gastroc = pd.read_excel("DEGs OSD101_419_401.xlsx", sheet_name="OSD-419 DGE", skiprows=1).dropna()
+
+# create a set for each and only include genes with an adjusted p-value < 0.05
+OSD_101_set = set(OSD_101_gastroc[OSD_101_gastroc["ADJP"] < 0.05]["Symbol"])
+OSD_401_set = set(OSD_401_gastroc[OSD_401_gastroc["ADJP"] < 0.05]["Symbol"])
+OSD_419_set = set(OSD_419_gastroc[OSD_419_gastroc["ADJP"] < 0.05]["Symbol"])
+
+# compare intersections (shared genes)
+shared_genes_in_all_3 = OSD_101_set & OSD_401_set & OSD_419_set
+shared_genes_101_419 = OSD_101_set & OSD_419_set
+shared_genes_401_419 = OSD_401_set & OSD_419_set
+shared_genes_101_401 = OSD_101_set & OSD_401_set
+
+# unique genes in each dataset
+unique_OSD_101 = OSD_101_set - (OSD_401_set | OSD_419_set)
+unique_OSD_401 = OSD_401_set - (OSD_101_set | OSD_419_set)
+unique_OSD_419 = OSD_419_set - (OSD_101_set | OSD_401_set)
+
+# count and list shared and unique genes
+print(f"Shared across all 3: {len(shared_genes_in_all_3)}")
+print(f"Shared across all 3: {shared_genes_in_all_3}")
+
+print(f"Shared between 101 & 401: {len(shared_genes_101_401)}")
+print(f"Shared between 101 & 401: {shared_genes_101_401}")
+
+print(f"Shared between 101 & 419: {len(shared_genes_101_419)}")
+print(f"Shared between 101 & 419: {shared_genes_101_419}")
+
+print(f"Shared between 401 & 419: {len(shared_genes_401_419)}")
+print(f"Shared between 401 & 419: {shared_genes_401_419}")
+
+print(f"Unique to OSD 101: {len(unique_OSD_101)}")
+print(f"Unique to OSD 101: {unique_OSD_101}")
+
+print(f"Unique to OSD 401: {len(unique_OSD_401)}")
+print(f"Unique to OSD 401: {unique_OSD_401}")
+
+print(f"Unique to OSD 419: {len(unique_OSD_419)}")
+print(f"Unique to OSD 419: {unique_OSD_419}")
+
+# classic venn diagram
+plt.figure(figsize=(6,6))
+venn3(
+    [OSD_101_set, OSD_401_set, OSD_419_set],
+    ("OSD-101", "OSD-401", "OSD-419"),
+    set_colors=['steelblue', 'orange', 'green'],
+    alpha=0.8
+)
+plt.title("Differentially Expressed Gastroc Genes (padj < 0.05)")
+# umcomment below line to save
+# plt.savefig("venn3_all_3_gastroc_OSD_SAC.png", dpi=300)
+
+# supervenn version
+venn_sets = [OSD_101_set, OSD_401_set, OSD_419_set]
+OSD_nums = ["OSD-101", "OSD-401", "OSD-419"]
+plt.figure(figsize=(8,3), dpi=300)
+supervenn(venn_sets, OSD_nums, rotate_col_annotations=True,
+          col_annotations_area_height=1.2, sets_ordering='minimize gaps',side_plots=False)
+plt.title("Differentially Expressed Gastroc Genes (padj < 0.05)")
+# umcomment below line to save
+# plt.savefig("supervenn_all_3_gastroc_OSD_DEGs_SAC.png", dpi=300)

Manuscript_Code/NBISC100_project/gastrocnemius/alignment_metrics/Gastrocnemius_Genome_mapping_clayton.R

@@ -1,49 +1,64 @@
-#load library
+# load libraries
 library(ggplot2)
 library(tidyverse)
 library(scales)
 library(RColorBrewer)
 
-#load data
+# load QC dataset
 gastroc.metrics.df <- read_csv("gastrocnemius_qc_metrics.csv")
 
-#prepare data for genome mapping
- genome.alignment.gastroc.metrics.df <- gastroc.metrics.df %>% select(osd_num,uniquely_mapped_percent,multimapped_percent,multimapped_toomany_percent,unmapped_tooshort_percent, unmapped_other_percent)
+# prepare data for genome mapping
+genome.alignment.gastroc.metrics.df <- gastroc.metrics.df %>%
+    select(osd_num,uniquely_mapped_percent,multimapped_percent,multimapped_toomany_percent,unmapped_tooshort_percent, unmapped_other_percent)
 
- #Add a new column for library_kit
- genome.alignment.gastroc.metrics.df <-  genome.alignment.gastroc.metrics.df %>%
-   mutate(library_kit= recode(osd_num,
-                              "OSD-401"="polyA-nonUPX kit",
-                              "OSD-101"="ribo-deplete kit",
-                              "OSD-419"="polyA-UPX kit"))
+# hard-coded OSD and kit orders
+osd_order <- c("OSD-401", "OSD-419", "OSD-101")
+kit_order <- c("polyA-nonUPX kit", "polyA-UPX kit", "ribo-deplete kit")
+
+# add a new column for library_kit
+genome.alignment.gastroc.metrics.df <- genome.alignment.gastroc.metrics.df %>%
+mutate(
+    osd_num = factor(osd_num, levels = osd_order),
+    library_kit= recode(osd_num,
+                        "OSD-401"="polyA-nonUPX kit",
+                        "OSD-101"="ribo-deplete kit",
+                        "OSD-419"="polyA-UPX kit"),
+    library_kit = factor(library_kit, levels = kit_order)
+)
 
 # Box_plot1: percentage of uniquely mapped data by library kit
- ggplot(genome.alignment.gastroc.metrics.df, aes(x = osd_num, y= uniquely_mapped_percent, fill = osd_num)) +
-   geom_boxplot(linewidth = 0.1, varwidth = TRUE) +
-   stat_boxplot(geom = "errorbar", width = 0.2, size= 0.1)+
-   facet_wrap(~library_kit, scales = "free_x", drop = TRUE)+
-   scale_y_continuous(breaks = pretty_breaks(n = 8))+
-   scale_fill_brewer(palette = "Set2")+
-   labs(title = "Uniquely Mapped Percentage of Gastroc Datasets by Library Kits", x = "OSD-number", y = "Uniquely mapped (%)") +
-   theme_classic() +
-   theme(legend.position = "none")+
-   theme(plot.title = element_text(hjust = 0.5))
-
- ggsave("genome_uniquely_mapped_gastroc_LibraryKits_SAC.png", dpi = 300,
-        width = 7, height = 4, units = "in")
-
-
- # Box_plot2: percentage of total mapped data by library kit
- ggplot(genome.alignment.gastroc.metrics.df, aes(x = osd_num, y= 100 - (unmapped_tooshort_percent + unmapped_other_percent), fill = osd_num)) +
-   geom_boxplot(linewidth = 0.1, varwidth = TRUE) +
-   stat_boxplot(geom = "errorbar", width = 0.2, size= 0.1)+
-   facet_wrap(~library_kit, scales = "free_x", drop = TRUE)+
-   scale_y_continuous(breaks = pretty_breaks(n = 10))+
-   scale_fill_brewer(palette = "Set2")+
-   labs(title = "Total Mapped Percentage of Gastroc Datasets by Library Kits", x = "OSD-number", y = "Total mapped (%)") +
-   theme_classic() +
-   theme(legend.position = "none")+
-   theme(plot.title = element_text(hjust = 0.5))
-
- ggsave("genome_total_mapped_gastroc_LibraryKits_SAC.png", dpi = 300,
-        width = 7, height = 4, units = "in")
+ggplot(genome.alignment.gastroc.metrics.df, aes(x = osd_num, y= uniquely_mapped_percent, fill = osd_num)) +
+    geom_boxplot(linewidth = 0.1, varwidth = TRUE) +
+    stat_boxplot(geom = "errorbar", width = 0.2, size= 0.1)+
+    facet_wrap(~library_kit, scales = "free_x", drop = TRUE)+
+    scale_y_continuous(breaks = pretty_breaks(n = 8))+
+    scale_fill_brewer(palette = "Set2")+
+    labs(title = "Uniquely Mapped Percentage of Gastroc Datasets by Library Kits", x = "OSD-number", y = "Uniquely mapped (%)") +
+    theme_classic() +
+    theme(legend.position = "none")+
+    theme(
+        plot.title = element_text(hjust = 0.5, size = 11),
+        plot.title.position = "plot"
+    )
+
+ggsave("genome_uniquely_mapped_gastroc_LibraryKits_SAC.png", dpi = 300,
+    width = 7, height = 4, units = "in")
+
+
+# Box_plot2: percentage of total mapped data by library kit
+ggplot(genome.alignment.gastroc.metrics.df, aes(x = osd_num, y= 100 - (unmapped_tooshort_percent + unmapped_other_percent), fill = osd_num)) +
+    geom_boxplot(linewidth = 0.1, varwidth = TRUE) +
+    stat_boxplot(geom = "errorbar", width = 0.2, size= 0.1)+
+    facet_wrap(~library_kit, scales = "free_x", drop = TRUE)+
+    scale_y_continuous(breaks = pretty_breaks(n = 10))+
+    scale_fill_brewer(palette = "Set2")+
+    labs(title = "Total Mapped Percentage of Gastroc Datasets by Library Kits", x = "OSD-number", y = "Total mapped (%)") +
+    theme_classic() +
+    theme(legend.position = "none")+
+    theme(
+        plot.title = element_text(hjust = 0.5, size = 11),
+        plot.title.position = "plot"
+    )
+
+ggsave("genome_total_mapped_gastroc_LibraryKits_SAC.png", dpi = 300,
+    width = 7, height = 4, units = "in")

Manuscript_Code/NBISC100_project/gastrocnemius/alignment_metrics/Gastrocnemius_Transcriptome_mapping_clayton.R

@@ -1,48 +1,63 @@
-#load library
+# load library
 library(ggplot2)
 library(tidyverse)
 library(scales)
 library(RColorBrewer)
 
-#load data
+# load QC dataset
 gastroc.metrics.df <- read_csv("gastrocnemius_qc_metrics.csv")
 
-#prepare data for transcriptome alignment
-transcriptome.alignment.gastroc.metrics.df <- gastroc.metrics.df %>% select(osd_num, pct_uniquely_aligned, pct_multi_aligned, pct_unalignable)
+# prepare data for transcriptome alignment
+transcriptome.alignment.gastroc.metrics.df <- gastroc.metrics.df %>%
+    select(osd_num, pct_uniquely_aligned, pct_multi_aligned, pct_unalignable)
 
-#Add new a column for library_kit
-transcriptome.alignment.gastroc.metrics.df <-  transcriptome.alignment.gastroc.metrics.df %>%
-  mutate(library_kit= recode(osd_num,
-                            "OSD-401"="polyA-nonUPX kit",
-                            "OSD-101"="ribo-deplete kit",
-                            "OSD-419"="polyA-UPX kit"))
+# hard-coded OSD and kit orders
+osd_order <- c("OSD-401", "OSD-419", "OSD-101")
+kit_order <- c("polyA-nonUPX kit", "polyA-UPX kit", "ribo-deplete kit")
+
+# add a new column for library_kit
+transcriptome.alignment.gastroc.metrics.df <- transcriptome.alignment.gastroc.metrics.df %>%
+mutate(
+    osd_num = factor(osd_num, levels = osd_order),
+    library_kit= recode(osd_num,
+                        "OSD-401"="polyA-nonUPX kit",
+                        "OSD-101"="ribo-deplete kit",
+                        "OSD-419"="polyA-UPX kit"),
+    library_kit = factor(library_kit, levels = kit_order)
+)
 
 # Box_plot1: percentage of uniquely aligned data to transcriptome by library kit
 ggplot(transcriptome.alignment.gastroc.metrics.df, aes(x = osd_num, y= pct_uniquely_aligned, fill = osd_num)) +
-  geom_boxplot(linewidth = 0.1, varwidth = TRUE) +
-  stat_boxplot(geom = "errorbar", width = 0.2, size= 0.1)+
-  facet_wrap(~library_kit, scales = "free_x", drop = TRUE)+
-  scale_y_continuous(breaks = pretty_breaks(n = 8))+
-  scale_fill_brewer(palette = "Set2")+
-  labs(title = "Uniquely Aligned Percentage of Gastroc Datasets by Library Kits", x = "OSD-number", y = "Uniquely aligned (%)") +
-  theme_classic() +
-  theme(legend.position = "none")+
-  theme(plot.title = element_text(hjust = 0.5))
+    geom_boxplot(linewidth = 0.1, varwidth = TRUE) +
+    stat_boxplot(geom = "errorbar", width = 0.2, size= 0.1)+
+    facet_wrap(~library_kit, scales = "free_x", drop = TRUE)+
+    scale_y_continuous(breaks = pretty_breaks(n = 8))+
+    scale_fill_brewer(palette = "Set2")+
+    labs(title = "Uniquely Aligned Percentage of Gastroc Datasets by Library Kits", x = "OSD-number", y = "Uniquely aligned (%)") +
+    theme_classic() +
+    theme(legend.position = "none")+
+    theme(
+        plot.title = element_text(hjust = 0.5, size = 11),
+        plot.title.position = "plot"
+    )
 
 ggsave("transcriptome_unique_alignment_gastroc_libraryKits_SAC.png", dpi = 300,
-       width = 7.3, height = 4, units = "in")
+    width = 7.3, height = 4, units = "in")
 
 # Box_plot2: percentage of total aligned data to transcriptome by library kit
 ggplot(transcriptome.alignment.gastroc.metrics.df, aes(x = osd_num, y= (pct_uniquely_aligned+pct_multi_aligned), fill = osd_num)) +
-  geom_boxplot(linewidth = 0.1, varwidth = TRUE) +
-  stat_boxplot(geom = "errorbar", width = 0.2, size= 0.1)+
-  facet_wrap(~library_kit, scales = "free_x", drop = TRUE)+
-  scale_y_continuous(breaks = pretty_breaks(n = 10))+
-  scale_fill_brewer(palette = "Set2")+
-  labs(title = "Total Aligned Percentage of Gastroc Datasets by Library Kits", x = "OSD-number", y = "Total aligned (%)") +
-  theme_classic() +
-  theme(legend.position = "none")+
-  theme(plot.title = element_text(hjust = 0.5))
+    geom_boxplot(linewidth = 0.1, varwidth = TRUE) +
+    stat_boxplot(geom = "errorbar", width = 0.2, size= 0.1)+
+    facet_wrap(~library_kit, scales = "free_x", drop = TRUE)+
+    scale_y_continuous(breaks = pretty_breaks(n = 10))+
+    scale_fill_brewer(palette = "Set2")+
+    labs(title = "Total Aligned Percentage of Gastroc Datasets by Library Kits", x = "OSD-number", y = "Total aligned (%)") +
+    theme_classic() +
+    theme(legend.position = "none")+
+    theme(
+        plot.title = element_text(hjust = 0.5, size = 11),
+        plot.title.position = "plot"
+    )
 
 ggsave("transcriptome_total_alignment_gastroc_libraryKits_SAC.png", dpi = 300,
-       width = 7.3, height = 4, units = "in")
+    width = 7.3, height = 4, units = "in")

Manuscript_Code/NBISC100_project/gastrocnemius/rRNA_contamination/Gastrocnemius_rRNA_contamination_clayton.R

@@ -1,34 +1,45 @@
-#load library
+# load library
 library(ggplot2)
 library(tidyverse)
 library(scales)
 library(RColorBrewer)
 
-#load data
+# load QC dataset
 gastroc.metrics.df <- read_csv("gastrocnemius_qc_metrics.csv")
 
 # prepare data for rRNA contamination
 rRNA_contamination.gastroc.metrics.df <-gastroc.metrics.df %>%
-  select(osd_num, rrna_contamination)
+    select(osd_num, rrna_contamination)
+
+# hard-coded OSD and kit orders
+osd_order <- c("OSD-401", "OSD-419", "OSD-101")
+kit_order <- c("polyA-nonUPX kit", "polyA-UPX kit", "ribo-deplete kit")
 
 # Add a new column for library kit
 rRNA_contamination.gastroc.metrics.df <- rRNA_contamination.gastroc.metrics.df %>%
-  mutate(library_kit= recode(osd_num,
-                            "OSD-401"="polyA-nonUPX kit",
-                            "OSD-101"="ribo-deplete kit",
-                            "OSD-419"="polyA-UPX kit"))
+mutate(
+    osd_num = factor(osd_num, levels = osd_order),
+    library_kit= recode(osd_num,
+                        "OSD-401"="polyA-nonUPX kit",
+                        "OSD-101"="ribo-deplete kit",
+                        "OSD-419"="polyA-UPX kit"),
+    library_kit = factor(library_kit, levels = kit_order)
+)
 
 # Box_plot1: gastroc_rRNA_contamination by library kit
 ggplot(rRNA_contamination.gastroc.metrics.df, aes(x = osd_num, y = rrna_contamination, fill= osd_num)) +
-  geom_boxplot(size = 0.1, varwidth = TRUE) +
-  stat_boxplot(geom = "errorbar", width = 0.2, size= 0.1)+
-  facet_wrap(~library_kit, scales = "free_x", drop = TRUE)+
-  scale_y_continuous(breaks = pretty_breaks(n = 10))+
-  scale_fill_brewer(palette = "Set2") +
-  labs(title = "rRNA Contamination of Gastroc Datasets by Library Kits", x = "OSD-number", y = "rRNA contamination %") +
-  theme_classic() +
-  theme(legend.position = "none")+
-  theme(plot.title = element_text(hjust = 0.5))
+    geom_boxplot(size = 0.1, varwidth = TRUE) +
+    stat_boxplot(geom = "errorbar", width = 0.2, size= 0.1)+
+    facet_wrap(~library_kit, scales = "free_x", drop = TRUE)+
+    scale_y_continuous(breaks = pretty_breaks(n = 10))+
+    scale_fill_brewer(palette = "Set2") +
+    labs(title = "rRNA Contamination of Gastroc Datasets by Library Kits", x = "OSD-number", y = "rRNA contamination %") +
+    theme_classic() +
+    theme(legend.position = "none")+
+    theme(
+        plot.title = element_text(hjust = 0.5, size = 11),
+        plot.title.position = "plot"
+    )
 
 ggsave("rRNA_contamination_gastroc_libraryKits_SAC.png", dpi = 300,
-       width = 7, height = 4, units = "in")
+    width = 7, height = 4, units = "in")