From bac2cb5c3add09c456a7cdb09dfd59ecc6959710 Mon Sep 17 00:00:00 2001 From: sclayton33 Date: Thu, 19 Mar 2026 16:40:52 -0400 Subject: [PATCH 1/2] add gastroc dge comparison code --- .../DGE_comparison_gastroc_clayton.py | 72 +++++++++++++++++++ 1 file changed, 72 insertions(+) create mode 100644 Manuscript_Code/NBISC100_project/gastrocnemius/DGE_assessment/DGE_comparison_gastroc_clayton.py diff --git a/Manuscript_Code/NBISC100_project/gastrocnemius/DGE_assessment/DGE_comparison_gastroc_clayton.py b/Manuscript_Code/NBISC100_project/gastrocnemius/DGE_assessment/DGE_comparison_gastroc_clayton.py new file mode 100644 index 0000000..e97ebb6 --- /dev/null +++ b/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) \ No newline at end of file From 2295b54aef6337297e37dcbbceddf3691176c9cc Mon Sep 17 00:00:00 2001 From: sclayton33 Date: Sat, 11 Jul 2026 17:48:29 -0400 Subject: [PATCH 2/2] hard-code osd and kit number orders --- .../Gastrocnemius_Genome_mapping_clayton.R | 93 +++++++----- ...trocnemius_Transcriptome_mapping_clayton.R | 75 ++++++---- ...Gastrocnemius_rRNA_contamination_clayton.R | 45 +++--- .../Gastrocnemius_read_depth_clayton.R | 71 +++++---- .../Gastrocnemius_RseQC_genebody_clayton.R | 43 ++++-- ...trocnemius_RseQC_inferExperiment_clayton.R | 45 +++--- ...ocnemius_RseQC_read_distribution_clayton.R | 137 ++++++++++-------- 7 files changed, 306 insertions(+), 203 deletions(-) diff --git a/Manuscript_Code/NBISC100_project/gastrocnemius/alignment_metrics/Gastrocnemius_Genome_mapping_clayton.R b/Manuscript_Code/NBISC100_project/gastrocnemius/alignment_metrics/Gastrocnemius_Genome_mapping_clayton.R index 8861b7e..186f730 100644 --- a/Manuscript_Code/NBISC100_project/gastrocnemius/alignment_metrics/Gastrocnemius_Genome_mapping_clayton.R +++ b/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") diff --git a/Manuscript_Code/NBISC100_project/gastrocnemius/alignment_metrics/Gastrocnemius_Transcriptome_mapping_clayton.R b/Manuscript_Code/NBISC100_project/gastrocnemius/alignment_metrics/Gastrocnemius_Transcriptome_mapping_clayton.R index afb9464..81a7aa1 100644 --- a/Manuscript_Code/NBISC100_project/gastrocnemius/alignment_metrics/Gastrocnemius_Transcriptome_mapping_clayton.R +++ b/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") diff --git a/Manuscript_Code/NBISC100_project/gastrocnemius/rRNA_contamination/Gastrocnemius_rRNA_contamination_clayton.R b/Manuscript_Code/NBISC100_project/gastrocnemius/rRNA_contamination/Gastrocnemius_rRNA_contamination_clayton.R index cd9286a..5f5e7b8 100644 --- a/Manuscript_Code/NBISC100_project/gastrocnemius/rRNA_contamination/Gastrocnemius_rRNA_contamination_clayton.R +++ b/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") diff --git a/Manuscript_Code/NBISC100_project/gastrocnemius/read_depth/Gastrocnemius_read_depth_clayton.R b/Manuscript_Code/NBISC100_project/gastrocnemius/read_depth/Gastrocnemius_read_depth_clayton.R index fa89658..e25b764 100644 --- a/Manuscript_Code/NBISC100_project/gastrocnemius/read_depth/Gastrocnemius_read_depth_clayton.R +++ b/Manuscript_Code/NBISC100_project/gastrocnemius/read_depth/Gastrocnemius_read_depth_clayton.R @@ -1,48 +1,65 @@ -#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 specific read_depth plots -read_depth.gastroc.metrics.df <- gastroc.metrics.df %>% select(osd_num, read_depth, sequencing_instrument) +read_depth.gastroc.metrics.df <- gastroc.metrics.df %>% + select(osd_num, read_depth, sequencing_instrument) + +# 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") +sequencer_order <- c("Illumina HiSeq 3000", "Illumina NovaSeq 6000", "Illumina HiSeq 4000") #Add new columns for library_kit read_depth.gastroc.metrics.df <- read_depth.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), + sequencing_instrument = factor(sequencing_instrument, levels = sequencer_order) +) # Box_plot1: gastroc_read_depth by library kit ggplot(read_depth.gastroc.metrics.df, aes(x = osd_num, y = read_depth/1e6, 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 = 12))+ - scale_fill_brewer(palette = "Set2")+ - labs(title = "Read Depth of Gastroc Datasets by Library Kits", x = "OSD-number", y = "Read depth (millions)") + - 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 = 12))+ + scale_fill_brewer(palette = "Set2")+ + labs(title = "Read Depth of Gastroc Datasets by Library Kits", x = "OSD-number", y = "Read depth (millions)") + + theme_classic() + + theme(legend.position = "none")+ + theme( + plot.title = element_text(hjust = 0.5, size = 11), + plot.title.position = "plot" + ) ggsave("read_depth_gastroc_libraryKits_SAC.png", dpi = 300, - width = 6.7, height = 4, units = "in") + width = 6.7, height = 4, units = "in") # Box_plot2: gastroc_read_depth by sequencing instrument ggplot(read_depth.gastroc.metrics.df, aes(x = osd_num, y = read_depth/1e6, fill= osd_num)) + - geom_boxplot(size = 0.1, varwidth = TRUE) + - stat_boxplot(geom = "errorbar", width = 0.2, size= 0.1)+ - facet_wrap(~sequencing_instrument, scales = "free_x", drop = TRUE)+ - scale_y_continuous(breaks = pretty_breaks(n = 12))+ - scale_fill_brewer(palette = "Set2")+ - labs(title = "Read Depth of Gastroc Datasets by Sequencer", x = "OSD-number", y = "Read depth (millions)") + - 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(~sequencing_instrument, scales = "free_x", drop = TRUE) + + scale_y_continuous(breaks = pretty_breaks(n = 12))+ + scale_fill_brewer(palette = "Set2")+ + labs(title = "Read Depth of Gastroc Datasets by Sequencer", x = "OSD-number", y = "Read depth (millions)") + + theme_classic() + + theme(legend.position = "none")+ + theme( + plot.title = element_text(hjust = 0.5, size = 11), + plot.title.position = "plot" + ) ggsave("read_depth_gastroc_sequencer_SAC.png", dpi = 600, - width = 6.7, height = 6, units = "in") + width = 6.7, height = 6, units = "in") diff --git a/Manuscript_Code/NBISC100_project/gastrocnemius/rseqc_metrics/Gastrocnemius_RseQC_genebody_clayton.R b/Manuscript_Code/NBISC100_project/gastrocnemius/rseqc_metrics/Gastrocnemius_RseQC_genebody_clayton.R index d6b4f34..977aa6d 100644 --- a/Manuscript_Code/NBISC100_project/gastrocnemius/rseqc_metrics/Gastrocnemius_RseQC_genebody_clayton.R +++ b/Manuscript_Code/NBISC100_project/gastrocnemius/rseqc_metrics/Gastrocnemius_RseQC_genebody_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 RseQC genebody coverage RseQC_genebody.gastroc.metrics.df <-gastroc.metrics.df %>% select(osd_num, ratio_genebody_cov_3_to_5) +# 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 RseQC_genebody.gastroc.metrics.df <- RseQC_genebody.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: gatroc_RseQC_geneBody_coverage by library kit ggplot(RseQC_genebody.gastroc.metrics.df, aes(x = osd_num, y = ratio_genebody_cov_3_to_5, 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 = "Gene Body Coverage of Gastroc Datasets by Library Kits", x = "OSD-number", y = "3' to 5' gene body coverage ratio") + - 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 = "Gene Body Coverage of Gastroc Datasets by Library Kits", x = "OSD-number", y = "3' to 5' gene body coverage ratio") + + theme_classic() + + theme(legend.position = "none")+ + theme( + plot.title = element_text(hjust = 0.5, size = 11), + plot.title.position = "plot" + ) ggsave("RseQC_geneBodyCoverage_gastroc_libraryKits_SAC.png", dpi = 300, - width = 6.7, height = 4, units = "in") + width = 6.7, height = 4, units = "in") diff --git a/Manuscript_Code/NBISC100_project/gastrocnemius/rseqc_metrics/Gastrocnemius_RseQC_inferExperiment_clayton.R b/Manuscript_Code/NBISC100_project/gastrocnemius/rseqc_metrics/Gastrocnemius_RseQC_inferExperiment_clayton.R index 9e377e3..9c486d5 100644 --- a/Manuscript_Code/NBISC100_project/gastrocnemius/rseqc_metrics/Gastrocnemius_RseQC_inferExperiment_clayton.R +++ b/Manuscript_Code/NBISC100_project/gastrocnemius/rseqc_metrics/Gastrocnemius_RseQC_inferExperiment_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 RseQC infer experiment RseQC_inferExp.gastroc.metrics.df <-gastroc.metrics.df %>% - select(osd_num, pct_sense) + select(osd_num, pct_sense) + +# 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 RseQC_inferExp.gastroc.metrics.df <- RseQC_inferExp.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_RseQC_strandedness by library kit ggplot(RseQC_inferExp.gastroc.metrics.df, aes(x = osd_num, y = pct_sense, 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 = "Strandedness of Gastroc Datasets by Library Kits", x = "OSD-number", y = "sense(%)") + - 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 = "Strandedness of Gastroc Datasets by Library Kits", x = "OSD-number", y = "sense(%)") + + theme_classic() + + theme(legend.position = "none")+ + theme( + plot.title = element_text(hjust = 0.5, size = 11), + plot.title.position = "plot" + ) ggsave("RseQC_strandedned_gastroc_libraryKits_SAC.png", dpi = 300, - width = 6.7, height = 4, units = "in") + width = 6.7, height = 4, units = "in") diff --git a/Manuscript_Code/NBISC100_project/gastrocnemius/rseqc_metrics/Gastrocnemius_RseQC_read_distribution_clayton.R b/Manuscript_Code/NBISC100_project/gastrocnemius/rseqc_metrics/Gastrocnemius_RseQC_read_distribution_clayton.R index f07be70..a329710 100644 --- a/Manuscript_Code/NBISC100_project/gastrocnemius/rseqc_metrics/Gastrocnemius_RseQC_read_distribution_clayton.R +++ b/Manuscript_Code/NBISC100_project/gastrocnemius/rseqc_metrics/Gastrocnemius_RseQC_read_distribution_clayton.R @@ -1,99 +1,122 @@ -#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 RseQC read distribution RseQC_read_distribution.gastroc.metrics.df <-gastroc.metrics.df %>% - select(osd_num, cds_exons_pct,`5_utr_exons_pct`,`3_utr_exons_pct`,introns_pct, tss_up_5kb_10kb_pct) + select(osd_num, cds_exons_pct,`5_utr_exons_pct`,`3_utr_exons_pct`,introns_pct, tss_up_5kb_10kb_pct) + +# 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 RseQC_read_distribution.gastroc.metrics.df <- RseQC_read_distribution.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) +) ## cds exons---------------------- # Box_plot1: gastroc_RseQC_cds_exons by library kit ggplot(RseQC_read_distribution.gastroc.metrics.df, aes(x = osd_num, y = cds_exons_pct, 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 = "Read Distribution (CDS Exons) of Gastroc Datasets by Library Kits", x = "OSD-number", y = "cds exons(%)") + - 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 = "Read Distribution (CDS Exons) of Gastroc Datasets by Library Kits", x = "OSD-number", y = "cds exons(%)") + + theme_classic() + + theme(legend.position = "none")+ + theme( + plot.title = element_text(hjust = 0.5, size = 11), + plot.title.position = "plot" + ) ggsave("RseQC_cds_exons_gastroc_libraryKits_SAC.png", dpi = 300, - width = 7.7, height = 4, units = "in") + width = 7.7, height = 4, units = "in") ## 5'UTRs--------------- # Box_plot1: gastroc_RseQC_5'UTR by library kit ggplot(RseQC_read_distribution.gastroc.metrics.df, aes(x = osd_num, y = `5_utr_exons_pct`, 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 = "Read Distribution (5' UTRs) of Gastroc Datasets by Library Kits", x = "OSD-number", y = "5' UTRs (%)") + - 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 = "Read Distribution (5' UTRs) of Gastroc Datasets by Library Kits", x = "OSD-number", y = "5' UTRs (%)") + + theme_classic() + + theme(legend.position = "none")+ + theme( + plot.title = element_text(hjust = 0.5, size = 11), + plot.title.position = "plot" + ) ggsave("RseQC_5UTRs_gastroc_libraryKits_SAC.png", dpi = 300, - width = 7.7, height = 4, units = "in") + width = 7.7, height = 4, units = "in") ## 3' UTRs -------------- # Box_plot1: gastroc_RseQC_3'UTR by library kit ggplot(RseQC_read_distribution.gastroc.metrics.df, aes(x = osd_num, y = `3_utr_exons_pct`, 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 = "Read Distribution (3' UTRs) of Gastroc Datasets by Library Kits", x = "OSD-number", y = "3' UTRs (%)") + - 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 = "Read Distribution (3' UTRs) of Gastroc Datasets by Library Kits", x = "OSD-number", y = "3' UTRs (%)") + + theme_classic() + + theme(legend.position = "none")+ + theme( + plot.title = element_text(hjust = 0.5, size = 11), + plot.title.position = "plot" + ) ggsave("RseQC_3UTRs_gastroc_libraryKits_SAC.png", dpi = 300, - width = 7.7, height = 4, units = "in") + width = 7.7, height = 4, units = "in") ## introns ------------ # Box_plot1: gastroc_RseQC_3'UTR by library kit ggplot(RseQC_read_distribution.gastroc.metrics.df, aes(x = osd_num, y = introns_pct, 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 = "Read Distribution (Introns) of Gastroc Datasets by Library Kits", x = "OSD-number", y = "introns (%)") + - 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 = "Read Distribution (Introns) of Gastroc Datasets by Library Kits", x = "OSD-number", y = "introns (%)") + + theme_classic() + + theme(legend.position = "none")+ + theme( + plot.title = element_text(hjust = 0.5, size = 11), + plot.title.position = "plot" + ) ggsave("RseQC_introns_gastroc_libraryKits_SAC.png", dpi = 300, - width = 7.7, height = 4, units = "in") + width = 7.7, height = 4, units = "in") ## tss_up 5 to 10 kb ------- # Box_plot1: gastroc_RseQC_introns by library kit ggplot(RseQC_read_distribution.gastroc.metrics.df, aes(x = osd_num, y = tss_up_5kb_10kb_pct, 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 = "Read Distribution (TSS up 5-10 kb) of Gastroc Datasets by Library Kits", x = "OSD-number", y = "tss up 5-10 kb (%)") + - 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 = "Read Distribution (TSS up 5-10 kb) of Gastroc Datasets by Library Kits", x = "OSD-number", y = "tss up 5-10 kb (%)") + + theme_classic() + + theme(legend.position = "none")+ + theme( + plot.title = element_text(hjust = 0.5, size = 11), + plot.title.position = "plot" + ) ggsave("RseQC_tssUp5_gastroc_libraryKits_SAC.png", dpi = 300, - width = 8.3, height = 4, units = "in") + width = 8.3, height = 4, units = "in")