DXR DE Analysis

workflowr

Ignored files:
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    data/Anat HM/
    Ignored:    data/CDKN1A_geneplot_Dox.RDS
    Ignored:    data/Cormotif_dox.RDS
    Ignored:    data/Cormotif_prob_gene_list.RDS
    Ignored:    data/Cormotif_prob_gene_list_doxonly.RDS
    Ignored:    data/DMSO_TNN13_plot.RDS
    Ignored:    data/DOX_TNN13_plot.RDS
    Ignored:    data/DOXgeneplots.RDS
    Ignored:    data/ExpressionMatrix_EMP.csv
    Ignored:    data/Ind1_DOX_Spearman_plot.RDS
    Ignored:    data/Ind6REP_Spearman_list.csv
    Ignored:    data/Ind6REP_Spearman_plot.RDS
    Ignored:    data/Ind6REP_Spearman_set.csv
    Ignored:    data/Ind6_Spearman_plot.RDS
    Ignored:    data/MDM2_geneplot_Dox.RDS
    Ignored:    data/Metadata.csv
    Ignored:    data/SIRT1_geneplot_Dox.RDS
    Ignored:    data/Sample_annotated.csv
    Ignored:    data/SpearmanHeatmapMatrix_EMP
    Ignored:    data/VennDiagram_DEGs.csv.2025-05-19_13-04-29.204911.log
    Ignored:    data/VolcanoPlot_D144R_original_EMP_250623.png
    Ignored:    data/VolcanoPlot_D24R_original_EMP_250623.png
    Ignored:    data/VolcanoPlot_D24T_original_EMP_250623.png
    Ignored:    data/annot_dox.RDS
    Ignored:    data/annot_list_hm.csv
    Ignored:    data/cormotifARclust_pp.RDS
    Ignored:    data/cormotif_dxr_1.RDS
    Ignored:    data/cormotif_dxr_2.RDS
    Ignored:    data/counts/
    Ignored:    data/counts_DE_df_dox.RDS
    Ignored:    data/counts_DE_raw_data.RDS
    Ignored:    data/counts_raw_filt.RDS
    Ignored:    data/counts_raw_matrix.RDS
    Ignored:    data/counts_raw_matrix_EMP_250514.csv
    Ignored:    data/d24_Spearman_plot.RDS
    Ignored:    data/dge_calc_dxr.RDS
    Ignored:    data/dge_calc_matrix.RDS
    Ignored:    data/ensembl_backup_dox.RDS
    Ignored:    data/fC_AllCounts.RDS
    Ignored:    data/fC_DOXCounts.RDS
    Ignored:    data/featureCounts_Concat_Matrix_DOXSamples_EMP_250430.csv
    Ignored:    data/filcpm_colnames_matrix.csv
    Ignored:    data/filcpm_matrix.csv
    Ignored:    data/filt_gene_list_dox.RDS
    Ignored:    data/filter_gene_list_final.RDS
    Ignored:    data/final_data/
    Ignored:    data/gene_clustlike_motif.RDS
    Ignored:    data/gene_postprob_motif.RDS
    Ignored:    data/genedf_dxr.RDS
    Ignored:    data/genematrix_dox.RDS
    Ignored:    data/genematrix_dxr.RDS
    Ignored:    data/heartgenes.csv
    Ignored:    data/heartgenes_dox.csv
    Ignored:    data/heatmap_group_Anat.RDS
    Ignored:    data/ind_num_dox.RDS
    Ignored:    data/ind_num_dxr.RDS
    Ignored:    data/initial_cormotif.RDS
    Ignored:    data/initial_cormotif_dox.RDS
    Ignored:    data/new/
    Ignored:    data/new_cormotif_dox.RDS
    Ignored:    data/plot_leg_d.RDS
    Ignored:    data/plot_leg_d_horizontal.RDS
    Ignored:    data/plot_leg_d_vertical.RDS
    Ignored:    data/process_gene_data_funct.RDS
    Ignored:    data/ruv/
    Ignored:    data/tableED_GOBP.RDS
    Ignored:    data/tableESR_GOBP_postprob.RDS
    Ignored:    data/tableLD_GOBP.RDS
    Ignored:    data/tableLR_GOBP_postprob.RDS
    Ignored:    data/tableNR_GOBP.RDS
    Ignored:    data/tableNR_GOBP_postprob.RDS
    Ignored:    data/table_motif1_GOBP_d.RDS
    Ignored:    data/table_motif1_genes_GOBP.RDS
    Ignored:    data/table_motif2_GOBP_d.RDS
    Ignored:    data/table_motif3_genes_GOBP.RDS
    Ignored:    data/top.table_V.D144r_dox.RDS
    Ignored:    data/top.table_V.D24_dox.RDS
    Ignored:    data/top.table_V.D24r_dox.RDS

editor_options: chunk_output_type: inline

I did this part with Sayan according to his analysis to ensure that my matrix was consistent with his - allowing me to work on downstream analysis #Read in Libraries

#Read in Counts Files from FeatureCounts

#load in libraries needed
#these counts files are from featureCounts, all saved as RDS objects

# ####Individual 1 - 84-1####
# Counts_84_DOX_24 <- readRDS("data/counts/Counts_84_DOX_24.RDS")
# Counts_84_DMSO_24 <- readRDS("data/counts/Counts_84_DMSO_24.RDS")
# Counts_84_DOX_24rec <- readRDS("data/counts/Counts_84_DOX_24rec.RDS")
# Counts_84_DMSO_24rec <- readRDS("data/counts/Counts_84_DMSO_24rec.RDS")
# Counts_84_DOX_144rec <- readRDS("data/counts/Counts_84_DOX_144rec.RDS")
# Counts_84_DMSO_144rec <- readRDS("data/counts/Counts_84_DMSO_144rec.RDS")
# 
# ####Individual 2 - 87-1####
# Counts_87_DOX_24 <- readRDS("data/counts/Counts_87_DOX_24.RDS")
# Counts_87_DMSO_24 <- readRDS("data/counts/Counts_87_DMSO_24.RDS")
# Counts_87_DOX_24rec <- readRDS("data/counts/Counts_87_DOX_24rec.RDS")
# Counts_87_DMSO_24rec <- readRDS("data/counts/Counts_87_DMSO_24rec.RDS")
# Counts_87_DOX_144rec <- readRDS("data/counts/Counts_87_DOX_144rec.RDS")
# Counts_87_DMSO_144rec <- readRDS("data/counts/Counts_87_DMSO_144rec.RDS")
# 
# ####Individual 3 - 78-1####
# Counts_78_DOX_24 <- readRDS("data/counts/Counts_78_DOX_24.RDS")
# Counts_78_DMSO_24 <- readRDS("data/counts/Counts_78_DMSO_24.RDS")
# Counts_78_DOX_24rec <- readRDS("data/counts/Counts_78_DOX_24rec.RDS")
# Counts_78_DMSO_24rec <- readRDS("data/counts/Counts_78_DMSO_24rec.RDS")
# Counts_78_DOX_144rec <- readRDS("data/counts/Counts_78_DOX_144rec.RDS")
# Counts_78_DMSO_144rec <- readRDS("data/counts/Counts_78_DMSO_144rec.RDS")
# 
# ####Individual 4 - 75-1####
# Counts_75_DOX_24 <- readRDS("data/counts/Counts_75_DOX_24.RDS")
# Counts_75_DMSO_24 <- readRDS("data/counts/Counts_75_DMSO_24.RDS")
# Counts_75_DOX_24rec <- readRDS("data/counts/Counts_75_DOX_24rec.RDS")
# Counts_75_DMSO_24rec <- readRDS("data/counts/Counts_75_DMSO_24rec.RDS")
# Counts_75_DOX_144rec <- readRDS("data/counts/Counts_75_DOX_144rec.RDS")
# Counts_75_DMSO_144rec <- readRDS("data/counts/Counts_75_DMSO_144rec.RDS")
# 
# ####Individual 5 - 17-3####
# Counts_17_DOX_24 <- readRDS("data/counts/Counts_17_DOX_24.RDS")
# Counts_17_DMSO_24 <- readRDS("data/counts/Counts_17_DMSO_24.RDS")
# Counts_17_DOX_24rec <- readRDS("data/counts/Counts_17_DOX_24rec.RDS")
# Counts_17_DMSO_24rec <- readRDS("data/counts/Counts_17_DMSO_24rec.RDS")
# Counts_17_DOX_144rec <- readRDS("data/counts/Counts_17_DOX_144rec.RDS")
# Counts_17_DMSO_144rec <- readRDS("data/counts/Counts_17_DMSO_144rec.RDS")
# 
# ####Individual 6 - 90-1####
# Counts_90_DOX_24 <- readRDS("data/counts/Counts_90_DOX_24.RDS")
# Counts_90_DMSO_24 <- readRDS("data/counts/Counts_90_DMSO_24.RDS")
# Counts_90_DOX_24rec <- readRDS("data/counts/Counts_90_DOX_24rec.RDS")
# Counts_90_DMSO_24rec <- readRDS("data/counts/Counts_90_DMSO_24rec.RDS")
# Counts_90_DOX_144rec <- readRDS("data/counts/Counts_90_DOX_144rec.RDS")
# Counts_90_DMSO_144rec <- readRDS("data/counts/Counts_90_DMSO_144rec.RDS")
# 
# ####Individual 7 - 90-1REP####
# Counts_90REP_DOX_24 <- readRDS("data/counts/Counts_90REP_DOX_24.RDS")
# Counts_90REP_DMSO_24 <- readRDS("data/counts/Counts_90REP_DMSO_24.RDS")
# Counts_90REP_DOX_24rec <- readRDS("data/counts/Counts_90REP_DOX_24rec.RDS")
# Counts_90REP_DMSO_24rec <- readRDS("data/counts/Counts_90REP_DMSO_24rec.RDS")
# Counts_90REP_DOX_144rec <- readRDS("data/counts/Counts_90REP_DOX_144rec.RDS")
# Counts_90REP_DMSO_144rec <- readRDS("data/counts/Counts_90REP_DMSO_144rec.RDS")

##Create my Dataframe from Counts Files

# counts_raw_df <-
#    data.frame(
#      Counts_84_DOX_24,
#      Counts_84_DMSO_24$MCW_EMP_JT_R29_R1.bam,
#      Counts_84_DOX_24rec$MCW_EMP_JT_R30_R1.bam,
#      Counts_84_DMSO_24rec$MCW_EMP_JT_R32_R1.bam,
#      Counts_84_DOX_144rec$MCW_EMP_JT_R33_R1.bam,
#      Counts_84_DMSO_144rec$MCW_EMP_JT_R35_R1.bam,
#      Counts_87_DOX_24$MCW_EMP_JT_R36_R1.bam,
#      Counts_87_DMSO_24$MCW_EMP_JT_R38_R1.bam,
#      Counts_87_DOX_24rec$MCW_EMP_JT_R39_R1.bam,
#      Counts_87_DMSO_24rec$MCW_EMP_JT_R41_R1.bam,
#      Counts_87_DOX_144rec$MCW_EMP_JT_R42_R1.bam,
#      Counts_87_DMSO_144rec$MCW_EMP_JT_R44_R1.bam,
#      Counts_78_DOX_24$MCW_EMP_JT_R45_R1.bam,
#      Counts_78_DMSO_24$MCW_EMP_JT_R47_R1.bam,
#      Counts_78_DOX_24rec$MCW_EMP_JT_R48_R1.bam,
#      Counts_78_DMSO_24rec$MCW_EMP_JT_R50_R1.bam,
#      Counts_78_DOX_144rec$MCW_EMP_JT_R51_R1.bam,
#      Counts_78_DMSO_144rec$MCW_EMP_JT_R53_R1.bam,
#      Counts_75_DOX_24$MCW_EMP_JT_R54_R1.bam,
#      Counts_75_DMSO_24$MCW_EMP_JT_R56_R1.bam,
#      Counts_75_DOX_24rec$MCW_EMP_JT_R57_R1.bam,
#      Counts_75_DMSO_24rec$MCW_EMP_JT_R59_R1.bam,
#      Counts_75_DOX_144rec$MCW_EMP_JT_R60_R1.bam,
#      Counts_75_DMSO_144rec$MCW_EMP_JT_R62_R1.bam,
#      Counts_17_DOX_24$MCW_EMP_JT_R63_R1.bam,
#      Counts_17_DMSO_24$MCW_EMP_JT_R65_R1.bam,
#      Counts_17_DOX_24rec$MCW_EMP_JT_R66_R1.bam,
#      Counts_17_DMSO_24rec$MCW_EMP_JT_R68_R1.bam,
#      Counts_17_DOX_144rec$MCW_EMP_JT_R69_R1.bam,
#      Counts_17_DMSO_144rec$MCW_EMP_JT_R71_R1.bam,
#      Counts_90_DOX_24$MCW_EMP_JT_R72_R1.bam,
#      Counts_90_DMSO_24$MCW_EMP_JT_R74_R1.bam,
#      Counts_90_DOX_24rec$MCW_EMP_JT_R75_R1.bam,
#      Counts_90_DMSO_24rec$MCW_EMP_JT_R77_R1.bam,
#      Counts_90_DOX_144rec$MCW_EMP_JT_R78_R1.bam,
#      Counts_90_DMSO_144rec$MCW_EMP_JT_R80_R1.bam,
#      Counts_90REP_DOX_24$MCW_EMP_JT_R81_R1.bam,
#      Counts_90REP_DMSO_24$MCW_EMP_JT_R83_R1.bam,
#      Counts_90REP_DOX_24rec$MCW_EMP_JT_R84_R1.bam,
#      Counts_90REP_DMSO_24rec$MCW_EMP_JT_R86_R1.bam,
#      Counts_90REP_DOX_144rec$MCW_EMP_JT_R87_R1.bam,
#      Counts_90REP_DMSO_144rec$MCW_EMP_JT_R89_R1.bam
#    )

#now save this as a matrix
#counts_raw_matrix <- counts_raw_df %>% column_to_rownames(var = "X") %>% as.matrix()

counts_raw_matrix <- readRDS("data/new/counts_raw_matrix.RDS")

dim(counts_raw_matrix)

[1] 28395    42

#28395 is my initial amount of genes prior to filtering

#write this to a csv so I can save it for later
#write.csv(counts_raw_matrix, "C:/Users/emmap/RDirectory/Recovery_RNAseq/Recovery_5FU/data/new/counts_raw_matrix_EMP_250514.csv")

#I also want to save this as an R object so I don't have to run the counts every time
#saveRDS(counts_raw_matrix, "data/new/counts_raw_matrix.RDS")

###Include colors and factors for analysis

#I want to include the color schemes I have for my treatment, individuals, and timepoints
####Colors####
tx_col <- c("DOX" = "#499FBD", "DMSO" = "#BBBBBC")
col_tx_large <- rep(c("#499FBD" , "#BBBBBC"), 21)
col_tx_large_2 <- c(rep("#499FBD" , 3), rep("#BBBBBC", 3), 21)

ind_col <- c("#003F5C", "#45AE91",  "#58209D", "#8B3E9B", "#FF6361", "#BC4169", "#FF2362")

ind_col_norep <- c("#003F5C", "#45AE91",  "#58209D", "#8B3E9B", "#FF6361", "#BC4169")

time_col <- c("#238B45", "#74C476", "#C7E9C0")

cond_col <- c("#003F5C", "#45AE91",  "#58209D", "#8B3E9B", "#FF6361", "#BC4169")

#QC Mapping Plots

#this dataframe contains my alignment percentages from featureCounts
##already filtered to only include DOX + DMSO samples
fC_DOXCounts <- readRDS("data/fC_DOXCounts.RDS")

#Now I want to plot these values out

####Reads by Sample####
reads_by_sample <- c("DOX" = "#499FBD", "DMSO" = "#BBBBBC")
fC_DOXCounts %>% 
  ggplot(., aes (x = Conditions, y = Total_Align, fill = Treatment, group_by = Line))+
  geom_col()+
 geom_hline(aes(yintercept=20000000))+
 scale_fill_manual(values=reads_by_sample)+
  ggtitle(expression("Total number of reads by sample"))+
  xlab("")+
  ylab(expression("RNA-sequencing reads"))+
  theme_bw()+
  theme(plot.title = element_text(size = rel(2), hjust = 0.5),
        axis.title = element_text(size = 15, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.text.y = element_text(size =10, color = "black", angle = 0, hjust = 0.8, vjust = 0.5),
        axis.text.x = element_text(size =10, color = "black", angle = 90, hjust = 1, vjust = 0.2),
        #strip.text.x = element_text(size = 15, color = "black", face = "bold"),
        strip.text.y = element_text(color = "white"))

Past versions of "QC Mapping Plots-1.png"

####Read Counts by Treatment####
fC_DOXCounts %>% 
  ggplot(., aes (x =Treatment, y= Total_Align, fill = Treatment))+
  geom_boxplot()+
 scale_fill_manual(values=reads_by_sample)+
  ggtitle(expression("Total number of reads by treatment"))+
  xlab("")+
  ylab(expression("RNA-sequencing reads"))+
  theme_bw()+
  
  theme(plot.title = element_text(size = rel(2), hjust = 0.5),
        axis.title = element_text(size = 15, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.text.y = element_text(size =10, color = "black", angle = 0, hjust = 0.8, vjust = 0.5),
        axis.text.x = element_text(size =10, color = "black", angle = 90, hjust = 1, vjust = 0.2),
        #strip.text.x = element_text(size = 15, color = "black", face = "bold"),
        strip.text.y = element_text(color = "white"))

Past versions of "QC Mapping Plots-2.png"

####Total Reads Per Individual####
fC_DOXCounts %>% 
  ggplot(., aes (x =as.factor(Line), y=Total_Align))+
  geom_boxplot(aes(fill=as.factor(Line)))+
 scale_fill_brewer(palette = "Dark2", name = "Individual")+
  ggtitle(expression("Total number of reads by individual"))+
  xlab("")+
  ylab(expression("RNA-sequencing reads"))+
  theme_bw()+

  theme(plot.title = element_text(size = rel(2), hjust = 0.5),
        axis.title = element_text(size = 15, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.text.y = element_text(size =10, color = "black", angle = 0, hjust = 0.8, vjust = 0.5),
        axis.text.x = element_text(size =10, color = "black", angle = 0, hjust = 1),
        #strip.text.x = element_text(size = 15, color = "black", face = "bold"),
        strip.text.y = element_text(color = "white"))

Past versions of "QC Mapping Plots-3.png"

####Total Mapped Reads Per Drug####

reads_by_sample <- c("DOX" = "#499FBD", "DMSO" = "#BBBBBC")
fC_DOXCounts %>% 
  ggplot(., aes (x = Conditions, y = Assigned_Align, fill = Treatment, group_by = Line))+
  geom_col()+
 geom_hline(aes(yintercept=20000000))+
 scale_fill_manual(values=reads_by_sample)+
  ggtitle(expression("Total number of mapped reads by sample"))+
  xlab("")+
  ylab(expression("RNA-sequencing reads"))+
  theme_bw()+
  theme(plot.title = element_text(size = rel(2), hjust = 0.5),
        axis.title = element_text(size = 15, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.text.y = element_text(size =10, color = "black", angle = 0, hjust = 0.8, vjust = 0.5),
        axis.text.x = element_text(size =10, color = "black", angle = 90, hjust = 1, vjust = 0.2),
        #strip.text.x = element_text(size = 15, color = "black", face = "bold"),
        strip.text.y = element_text(color = "white"))

Past versions of "QC Mapping Plots-4.png"

####Read Counts by Treatment####
fC_DOXCounts %>% 
  ggplot(., aes (x =Treatment, y= Assigned_Align, fill = Treatment))+
  geom_boxplot()+
 scale_fill_manual(values=reads_by_sample)+
  ggtitle(expression("Total number of mapped reads by treatment"))+
  xlab("")+
  ylab(expression("RNA-sequencing reads"))+
  theme_bw()+
  
  theme(plot.title = element_text(size = rel(2), hjust = 0.5),
        axis.title = element_text(size = 15, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.text.y = element_text(size =10, color = "black", angle = 0, hjust = 0.8, vjust = 0.5),
        axis.text.x = element_text(size =10, color = "black", angle = 90, hjust = 1, vjust = 0.2),
        #strip.text.x = element_text(size = 15, color = "black", face = "bold"),
        strip.text.y = element_text(color = "white"))

Past versions of "QC Mapping Plots-5.png"

####Total Reads Per Individual####
fC_DOXCounts %>% 
  ggplot(., aes (x =as.factor(Line), y=Assigned_Align))+
  geom_boxplot(aes(fill=as.factor(Line)))+
 scale_fill_brewer(palette = "Dark2", name = "Individual")+
  ggtitle(expression("Total number of mapped reads by individual"))+
  xlab("")+
  ylab(expression("RNA-sequencing reads"))+
  theme_bw()+

  theme(plot.title = element_text(size = rel(2), hjust = 0.5),
        axis.title = element_text(size = 15, color = "black"),
        axis.ticks = element_line(linewidth = 1.5),
        axis.text.y = element_text(size =10, color = "black", angle = 0, hjust = 0.8, vjust = 0.5),
        axis.text.x = element_text(size =10, color = "black", angle = 0, hjust = 1),
        #strip.text.x = element_text(size = 15, color = "black", face = "bold"),
        strip.text.y = element_text(color = "white"))

Past versions of "QC Mapping Plots-6.png"

Now, I want to filter my dataframe Before I can filter by rowMeans, I must convert to log2cpm #Filter my Dataframe and Convert to log2cpm

#transform counts to cpm as a first step
counts_cpm_unfilt <- cpm(counts_raw_matrix, log = TRUE)
dim(counts_cpm_unfilt)

[1] 28395    42

#I should have 28395 genes here since this is unfiltered

hist(counts_cpm_unfilt,  
     main = "Histogram of Unfiltered Counts", 
     xlab = expression("Log"[2]*" counts-per-million"), 
     col = 4)

Past versions of "Dataframe Filtering-1.png"

###filter my data by rowMeans > 0 to exclude lowly expressed genes

filcpm_matrix <- subset(counts_cpm_unfilt, (rowMeans(counts_cpm_unfilt) > 0))
dim(filcpm_matrix)

[1] 14319    42

#I should have 14319 genes here

#now let's make a histogram of this to check the difference
hist(filcpm_matrix,  
     main = "Histogram of Filtered Counts by rowMeans > 0", 
     xlab = expression("Log"[2]*" counts-per-million"), 
     col = 2)

Past versions of "Dataframe Filtering-2.png"

#change the column names to match my samples - make sure that they are in the right order
#Individual 1 = 84-1 (M)
#Individual 2 = 87-1 (F)
#Individual 3 = 78-1 (F)
#Individual 4 = 75-1 (F)
#Individual 5 = 17-3 (M)
#Individual 6 = 90-1 (M)
#Individual 6REP = 90-1REP (M)

#Treatment/time should follow this order:
#DOX24tx
#DMSO24tx
#DOX24rec
#DMSO24rec
#DOX144rec
#DMSO144rec

colnames(filcpm_matrix) <- c("DOX_24T_Ind1",
                                "DMSO_24T_Ind1",
                                "DOX_24R_Ind1",
                                "DMSO_24R_Ind1",
                                "DOX_144R_Ind1",
                                "DMSO_144R_Ind1",
                                "DOX_24T_Ind2",
                                "DMSO_24T_Ind2",
                                "DOX_24R_Ind2",
                                "DMSO_24R_Ind2",
                                "DOX_144R_Ind2",
                                "DMSO_144R_Ind2",
                                "DOX_24T_Ind3",
                                "DMSO_24T_Ind3",
                                "DOX_24R_Ind3",
                                "DMSO_24R_Ind3",
                                "DOX_144R_Ind3",
                                "DMSO_144R_Ind3",
                                "DOX_24T_Ind4",
                                "DMSO_24T_Ind4",
                                "DOX_24R_Ind4",
                                "DMSO_24R_Ind4",
                                "DOX_144R_Ind4",
                                "DMSO_144R_Ind4",
                                "DOX_24T_Ind5",
                                "DMSO_24T_Ind5",
                                "DOX_24R_Ind5",
                                "DMSO_24R_Ind5",
                                "DOX_144R_Ind5",
                                "DMSO_144R_Ind5",
                                "DOX_24T_Ind6",
                                "DMSO_24T_Ind6",
                                "DOX_24R_Ind6",
                                "DMSO_24R_Ind6",
                                "DOX_144R_Ind6",
                                "DMSO_144R_Ind6",
                                "DOX_24T_Ind6REP",
                                "DMSO_24T_Ind6REP",
                                "DOX_24R_Ind6REP",
                                "DMSO_24R_Ind6REP",
                                "DOX_144R_Ind6REP",
                                "DMSO_144R_Ind6REP")

#export this as a csv
#write.csv(filcpm_matrix, "data/new/filcpm_final_matrix.csv")

##QC Boxplots of Filtered vs Unfiltered Data

#make boxplots of all counts vs log2cpm filtered counts

#set the margins so the x axis isn't cut off
##I don't mind if this one is partially cut off since all you need is the library number and not the whole name
par(mar = c(8,4,2,2))
#boxplot of unfiltered cpm matrix
boxplot(counts_cpm_unfilt, 
        main = "Boxplots of Unfiltered log2cpm", 
        names = colnames(counts_cpm_unfilt), 
        adj=1, las = 2, cex.axis = 0.7)

Past versions of "QC Filtering Boxplots-1.png"

#set the margins so the x axis isn't cut off
par(mar = c(8,4,2,2))
#boxplot of filtered cpm matrix
boxplot(filcpm_matrix, 
        main = "Boxplots of Filtered log2cpm (rowMeans > 0)", 
        names = colnames(filcpm_matrix), 
        adj=1, las = 2, cex.axis = 0.7)

Past versions of "QC Filtering Boxplots-2.png"

After making my final matrix, I pulled the gene symbols from the entrez IDs I had as my rownames I ran this initially and then moved the column into my final matrix My final matrix is called filcpm_final_matrix.csv saved under data #Put together my Data for Filtered Gene List

##I did this earlier so don't run again, I put the list into the filcpm_final_matrix.csv
# # ----------------- Load Required Libraries -----------------
# library(dplyr)
# library(readr)
# library(org.Hs.eg.db)
# library(AnnotationDbi)
# # ----------------- Load Data -----------------
# sample_data <- read_csv("data/filcpm_final_matrix.csv", show_col_types = FALSE)
# # ----------------- Ensure Entrez_ID is Present and in Character Format -----------------
# # Check column names
# print(colnames(sample_data))
# # Rename if needed (adjust if the column name is not exactly 'Entrez_ID')
# # sample_data <- sample_data %>% rename(Entrez_ID = `actual_column_name`)
# # Convert Entrez_ID to character
# sample_data <- sample_data %>%
#   mutate(Entrez_ID = as.character(Entrez_ID))
# # ----------------- Map Entrez_ID to Gene Symbol -----------------
# gene_symbols <- AnnotationDbi::select(
#   org.Hs.eg.db,
#   keys = sample_data$Entrez_ID,
#   columns = c("SYMBOL"),
#   keytype = "ENTREZID"
# )
# # ----------------- Join Back to Main Data -----------------
# sample_annotated <- left_join(sample_data, gene_symbols, by = c("Entrez_ID" = "ENTREZID"))
# # ----------------- Save Annotated Output -----------------
# #write_csv(sample_annotated, "data/Sample_annotated.csv")

#Since I ran this before, Sample_annotated.csv columns of EntrezID and Symbol have been copied into my final matrix - so disregard this file except for record-keeping

Now that I have my final matrix, I would like to check some key genes I want to make sure that these genes are responding as we expect We have triple checked this dataset to ensure that columns are in order #Check my Response Genes in log2cpm format

#Load in my count matrix
boxplot1 <- read.csv("data/new/filcpm_final_matrix.csv") %>% 
  as.data.frame()

#save boxplot1 as an object filcpm_matrix_genes
#saveRDS(boxplot1, "data/new/filcpm_matrix_genes.RDS")

#Define gene list(s)
initial_test_genes <- c("CDKN1A", "MDM2", "BAX", "RARG", "TP53", "TOP2B", "TOP2A")  
#Add more gene symbols as needed or add more categories

#Now put in the function I want to use to generate boxplots of genes
process_gene_data <- function(gene) {
  gene_data <- boxplot1 %>% filter(SYMBOL == gene)
  long_data <- gene_data %>%
    pivot_longer(cols = -c(Entrez_ID, SYMBOL), names_to = "Sample", values_to = "log2CPM") %>%
    mutate(
      Drug = case_when(
        grepl("DOX", Sample) ~ "DOX",
        grepl("DMSO", Sample) ~ "DMSO",
        TRUE ~ NA_character_
      ),
      Timepoint = case_when(
        grepl("_24T_", Sample) ~ "24T",
        grepl("_24R_", Sample) ~ "24R",
        grepl("_144R_", Sample) ~ "144R",
        TRUE ~ NA_character_
      ),
      Indv = case_when(
        grepl("Ind1$", Sample) ~ "1",
        grepl("Ind2$", Sample) ~ "2",
        grepl("Ind3$", Sample) ~ "3",
        grepl("Ind4$", Sample) ~ "4",
        grepl("Ind5$", Sample) ~ "5",
        grepl("Ind6$", Sample) ~ "6",
        grepl("Ind6REP$", Sample) ~ "6R",
        TRUE ~ NA_character_
      ),
      Condition = paste(Drug, Timepoint, sep = "_")
    )
  long_data$Condition <- factor(
    long_data$Condition,
    levels = c(
      "DOX_24T", 
      "DMSO_24T", 
      "DOX_24R", 
      "DMSO_24R", 
      "DOX_144R", 
      "DMSO_144R"
    )
  )
  return(long_data)
}

#this function is saved under process_gene_data so I will save as an R object

#saveRDS(process_gene_data, "data/new/process_gene_data_funct.RDS")

#Generate Boxplots from the above function using our gene list above
for (gene in initial_test_genes) {
  gene_data <- process_gene_data(gene)
  p <- ggplot(gene_data, aes(x = Condition, y = log2CPM, fill = Drug)) +
    geom_boxplot(outlier.shape = NA) +
    geom_point(aes(color = Indv), size = 2, alpha = 0.5, position = position_jitter(width = -1, height = 0)) +
    scale_fill_manual(values = c("DOX" = "#499FBD", "DMSO" = "#BBBBBC")) +
    ggtitle(paste("Log2CPM Expression of", gene)) +
    labs(x = "Treatment", y = "log2CPM") +
    theme_bw() +
    theme(
      plot.title = element_text(size = rel(2), hjust = 0.5),
      axis.title = element_text(size = 15, color = "black"),
      axis.text.x = element_text(size = 10, color = "black", angle = 90, hjust = 1)
    )
  print(p)
}

Past versions of "Check Response Genes log2cpm-1.png"

Past versions of "Check Response Genes log2cpm-2.png"

Past versions of "Check Response Genes log2cpm-3.png"

Past versions of "Check Response Genes log2cpm-4.png"

Past versions of "Check Response Genes log2cpm-5.png"

Past versions of "Check Response Genes log2cpm-6.png"

Past versions of "Check Response Genes log2cpm-7.png"

Now I’ve confirmed with some boxplots that my genes are present and (mostly) behaving as they should - Sayan’s CDKN1A and MDM2 are initially high at 24hr in DOX 0.5 - My CDKN1A and MDM2 are similar to DMSO at 24hr DOX 0.5 - These genes increase at DOX24R - These genes are also high at DOX144R but not as high as 24R However, TP53 and BAX are acting similarly across our data #PCA Analysis

#Now I want to check if my data is as expected on a PCA plot
#perform PCA calculations
prcomp_res_unfilt <- prcomp(t(counts_cpm_unfilt %>% as.matrix()), center =  TRUE)

prcomp_res_filt <- prcomp(t(filcpm_matrix %>% as.matrix()), center =  TRUE)

#read in my metadata annotations
Metadata <- read.csv("data/new/Metadata.csv")

#add in labels for individual numbers
ind_num <- c("1", "1", "1", "1", "1", "1", 
             "2", "2", "2", "2", "2", "2", 
             "3", "3", "3", "3", "3", "3", 
             "4", "4", "4", "4", "4", "4", 
             "5", "5", "5", "5", "5", "5", 
             "6", "6", "6", "6", "6", "6", 
             "6R", "6R", "6R", "6R", "6R", "6R")

# saveRDS(ind_num, "data/new/ind_num.RDS")

#now plot my PCA for unfiltered log2cpm
####PC1/PC2####
ggplot2::autoplot(prcomp_res_unfilt, data = Metadata, colour = "Condition", shape = "Time", size =4, x=1, y=2) +
  ggrepel::geom_text_repel(label=ind_num, vjust = -.5, max.overlaps = 30) +
  scale_color_manual(values=cond_col) +
  ggtitle(expression("PCA of Unfiltered log"[2]*"cpm")) +
  theme_bw()

Past versions of "PCA Analysis-1.png"

####PC2/PC3####
ggplot2::autoplot(prcomp_res_unfilt, data = Metadata, colour = "Condition", shape = "Time", size =4, x=2, y=3) +
  ggrepel::geom_text_repel(label=ind_num, vjust = -.5, max.overlaps = 30) +
  scale_color_manual(values=cond_col) +
  ggtitle(expression("PCA of Unfiltered log"[2]*"cpm")) +
  theme_bw()

Past versions of "PCA Analysis-2.png"

####PC3/PC4####
ggplot2::autoplot(prcomp_res_unfilt, data = Metadata, colour = "Condition", shape = "Time", size =4, x=3, y=4) +
  ggrepel::geom_text_repel(label=ind_num, vjust = -.5, max.overlaps = 30) +
  scale_color_manual(values=cond_col) +
  ggtitle(expression("PCA of Unfiltered log"[2]*"cpm")) +
  theme_bw()

Past versions of "PCA Analysis-3.png"

#Now plot my PCA for filtered log2cpm
####PC1/PC2####
ggplot2::autoplot(prcomp_res_filt, data = Metadata, colour = "Condition", shape = "Time", size =4, x=1, y=2) +
  ggrepel::geom_text_repel(label=ind_num, vjust = -.5, max.overlaps = 30) +
  scale_color_manual(values=cond_col) +
  ggtitle(expression("PCA of Filtered log"[2]*"cpm")) +
  theme_bw()

Past versions of "PCA Analysis-4.png"

ggplot2::autoplot(prcomp_res_filt, data = Metadata, colour = "Condition", shape = "Time", size =4, x=1, y=2) +
  ggrepel::geom_text_repel(label=ind_num, vjust = -.5, max.overlaps = 30) +
  scale_color_manual(values=cond_col) +
  ggtitle(expression("PCA of Filtered log"[2]*"cpm")) +
  theme_bw()

Past versions of "PCA Analysis-5.png"

####PC2/PC3####
ggplot2::autoplot(prcomp_res_filt, data = Metadata, colour = "Condition", shape = "Time", size =4, x=2, y=3) +
  ggrepel::geom_text_repel(label=ind_num, vjust = -.5, max.overlaps = 30) +
  scale_color_manual(values=cond_col) +
  ggtitle(expression("PCA of Filtered log"[2]*"cpm")) +
  theme_bw()

Past versions of "PCA Analysis-6.png"

####PC3/PC4####
ggplot2::autoplot(prcomp_res_filt, data = Metadata, colour = "Condition", shape = "Time", size =4, x=3, y=4) +
  ggrepel::geom_text_repel(label=ind_num, vjust = -.5, max.overlaps = 30) +
  scale_color_manual(values=cond_col) +
  ggtitle(expression("PCA of Filtered log"[2]*"cpm")) +
  theme_bw()

Past versions of "PCA Analysis-7.png"

#Correlation Heatmaps

#check to make sure that the column names are correct
lcpm_2 <- filcpm_matrix
colnames(lcpm_2) <- Metadata$Final_sample_name

#compute the correlation matrices, one pearson and one spearman
cor_matrix_pearson <- cor(lcpm_2, 
                          y = NULL, 
                          use = "everything",  
                          method = "pearson")
cor_matrix_spearman <- cor(lcpm_2,
                           y = NULL,
                           use = "everything",
                           method = "spearman")

# Extract metadata columns
Individual <- as.character(Metadata$Ind)
Time <- as.character(Metadata$Time)
Treatment <- as.character(Metadata$Drug)

# Define color palettes for annotations
annot_col_cor = list(drugs = c("DOX" = "#499FBD", 
                               "DMSO" = "#BBBBBC"),
                     individuals = c("1" = "#003F5C", 
                                     "2" = "#45AE91", 
                                     "3" = "#58209D", 
                                     "4" = "#8B3E9B", 
                                     "5" = "#FF6361", 
                                     "6" = "#BC4169", 
                                     "6R" = "#FF2362"),
                     timepoints = c("24T" = "#238B45",
                                    "24R" = "#74C476", 
                                    "144R" = "#C7E9C0"))

drug_colors <- c("DOX" = "#499FBD", 
                 "DMSO" = "#BBBBBC")
ind_colors <- c("1" = "red",
                "2" = "orange",
                "3" = "yellow",
                "4" = "green",
                "5" = "blue",
                "6" = "violet",
                "6R" = "purple")
time_colors <- c("24T" = "#238B45",
                 "24R" = "#74C476",
                 "144R" = "#C7E9C0")

# Create annotations
top_annotation <- HeatmapAnnotation(
  Individual = Individual, 
  Time = Time,
  Treatment = Treatment,
  col = list(
    Individual = ind_colors, 
    Time = time_colors,
    Treatment = drug_colors
  )
)

####ANNOTATED HEATMAPS####
# pheatmap(cor_matrix_pearson, border_color = "black", legend = TRUE, angle_col = 90, display_numbers = FALSE, number_color = "black", fontsize = 10, fontsize_number = 5, annotation_col = top_annotation, annotation_colors = annot_col)

####Pearson Heatmap####
heatmap_pearson <- Heatmap(cor_matrix_pearson,
                           name = "Pearson",
                           top_annotation = top_annotation,
                           show_row_names = TRUE,
                           show_column_names = TRUE,
                           cluster_rows = TRUE,
                           cluster_columns = TRUE,
                           border = TRUE)

# Draw the heatmap
draw(heatmap_pearson)

Past versions of "Correlation Heatmaps-1.png"

####Spearman Heatmap####
heatmap_spearman <- Heatmap(cor_matrix_spearman,
                           name = "Spearman",
                           top_annotation = top_annotation,
                           show_row_names = TRUE,
                           show_column_names = TRUE,
                           cluster_rows = TRUE,
                           cluster_columns = TRUE,
                           border = TRUE)

# Draw the heatmap
draw(heatmap_spearman)

Past versions of "Correlation Heatmaps-2.png"

#Make my Filtered Gene List

#Now I want to make a filtered gene list (my rownames)
##I will use this to filter my counts for limma + Cormotif

filt_gene_list <- rownames(filcpm_matrix)
#save this filtered gene list as I'll use it to filter my counts
#saveRDS(filt_gene_list, "data/new/filt_gene_list.RDS")

##Filter my Counts For DE Analysis

counts_raw_matrix <- readRDS("data/new/counts_raw_matrix.RDS")
#change column names to match samples for my raw counts matrix
colnames(counts_raw_matrix) <- c("DOX_24T_Ind1",
                                "DMSO_24T_Ind1",
                                "DOX_24R_Ind1",
                                "DMSO_24R_Ind1",
                                "DOX_144R_Ind1",
                                "DMSO_144R_Ind1",
                                "DOX_24T_Ind2",
                                "DMSO_24T_Ind2",
                                "DOX_24R_Ind2",
                                "DMSO_24R_Ind2",
                                "DOX_144R_Ind2",
                                "DMSO_144R_Ind2",
                                "DOX_24T_Ind3",
                                "DMSO_24T_Ind3",
                                "DOX_24R_Ind3",
                                "DMSO_24R_Ind3",
                                "DOX_144R_Ind3",
                                "DMSO_144R_Ind3",
                                "DOX_24T_Ind4",
                                "DMSO_24T_Ind4",
                                "DOX_24R_Ind4",
                                "DMSO_24R_Ind4",
                                "DOX_144R_Ind4",
                                "DMSO_144R_Ind4",
                                "DOX_24T_Ind5",
                                "DMSO_24T_Ind5",
                                "DOX_24R_Ind5",
                                "DMSO_24R_Ind5",
                                "DOX_144R_Ind5",
                                "DMSO_144R_Ind5",
                                "DOX_24T_Ind6",
                                "DMSO_24T_Ind6",
                                "DOX_24R_Ind6",
                                "DMSO_24R_Ind6",
                                "DOX_144R_Ind6",
                                "DMSO_144R_Ind6",
                                "DOX_24T_Ind6REP",
                                "DMSO_24T_Ind6REP",
                                "DOX_24R_Ind6REP",
                                "DMSO_24R_Ind6REP",
                                "DOX_144R_Ind6REP",
                                "DMSO_144R_Ind6REP")


#subset my count matrix based on filtered CPM matrix
x <- counts_raw_matrix[row.names(filcpm_matrix),]
dim(x)

[1] 14319    42

#14319 genes as expected!
#this is still in counts form

#remove my replicate individual at this time
x_norep <- x[,1:36]

#modify my metadata to match
Metadata_2 <- Metadata[1:36,]
rownames(Metadata_2) <- Metadata_2$Sample_bam
colnames(x_norep) <- Metadata_2$Sample_ID
rownames(Metadata_2) <- Metadata_2$Sample_ID

Metadata_2$Condition <- make.names(Metadata_2$Condition)
Metadata_2$Ind <- as.character(Metadata_2$Ind)

#saveRDS(Metadata_2, "data/new/Metadata_2_norep.RDS")

#Perform Differential Expression Analysis

#create DGEList object
dge <- DGEList(counts = x_norep)
dge$samples$group <- factor(Metadata_2$Condition)
dge <- calcNormFactors(dge, method = "TMM")

#saveRDS(dge, "data/new/dge_matrix.RDS")

#check normalization factors from TMM normalization of LIBRARIES
dge$samples

                     group lib.size norm.factors
84-1_DOX_24        DOX_24T 23393931    0.9745263
84-1_DMSO_24      DMSO_24T 22853195    0.9565797
84-1_DOX_24+24     DOX_24R 23846995    1.1659432
84-1_DMSO_24+24   DMSO_24R 21299355    0.9649641
84-1_DOX_24+144   DOX_144R 18222568    0.9913625
84-1_DMSO_24+144 DMSO_144R 28115884    0.9653464
87-1_DOX_24        DOX_24T 19935097    1.0526605
87-1_DMSO_24      DMSO_24T 21302879    0.9773889
87-1_DOX_24+24     DOX_24R 25636959    1.0751043
87-1_DMSO_24+24   DMSO_24R 26319662    0.9940323
87-1_DOX_24+144   DOX_144R 23463426    0.9003102
87-1_DMSO_24+144 DMSO_144R 25840938    0.9888449
78-1_DOX_24        DOX_24T 23085807    0.7676077
78-1_DMSO_24      DMSO_24T 25610495    1.0077383
78-1_DOX_24+24     DOX_24R 18083930    1.1682704
78-1_DMSO_24+24   DMSO_24R 24331177    0.9906872
78-1_DOX_24+144   DOX_144R 19754391    0.9941834
78-1_DMSO_24+144 DMSO_144R 22641509    1.0010734
75-1_DOX_24        DOX_24T 20583626    1.0676786
75-1_DMSO_24      DMSO_24T 28166198    1.0031906
75-1_DOX_24+24     DOX_24R 25831427    1.1530208
75-1_DMSO_24+24   DMSO_24R 26081158    1.0058953
75-1_DOX_24+144   DOX_144R 24659898    0.9261599
75-1_DMSO_24+144 DMSO_144R 25412931    0.9703454
17-3_DOX_24        DOX_24T 22518848    0.9766893
17-3_DMSO_24      DMSO_24T 24589534    0.9612345
17-3_DOX_24+24     DOX_24R 24797547    1.1703079
17-3_DMSO_24+24   DMSO_24R 25977536    0.9509690
17-3_DOX_24+144   DOX_144R 27447106    0.9422729
17-3_DMSO_24+144 DMSO_144R 24893583    0.9356377
90-1_DOX_24        DOX_24T 25187428    1.0311957
90-1_DMSO_24      DMSO_24T 25630519    1.0283437
90-1_DOX_24+24     DOX_24R 26138399    1.1183471
90-1_DMSO_24+24   DMSO_24R 24430396    0.9988688
90-1_DOX_24+144   DOX_144R 23323463    0.9496884
90-1_DMSO_24+144 DMSO_144R 25424152    0.9872926

#create my design matrix for DE
design <- model.matrix(~ 0 + Metadata_2$Condition)
colnames(design) <- gsub("Metadata_2\\$Condition", "", colnames(design))

#take care that the matrix automatically sorts cols alphabetically
##currently DMSO144R, DMSO24R, DMSO24T, DOX144R, DOX24R, DOX24T

#run duplicate correlation for individual effect
corfit <- duplicateCorrelation(object = dge$counts, design = design, block = Metadata_2$Ind)

#voom transformation and plot
v <- voom(dge, design, block = Metadata_2$Ind, correlation = corfit$consensus.correlation, plot = TRUE)

Past versions of "Differential Expression Analysis-1.png"

#fit my linear model
fit <- lmFit(v, design, block = Metadata_2$Ind, correlation = corfit$consensus.correlation)

#make my contrast matrix to compare across tx and veh
contrast_matrix <- makeContrasts(
  V.D24T = DOX_24T - DMSO_24T,
  V.D24R = DOX_24R - DMSO_24R, 
  V.D144R = DOX_144R - DMSO_144R,
  levels = design
)

#apply these contrasts to compare DOX to DMSO VEH
fit2 <- contrasts.fit(fit, contrast_matrix)
fit2 <- eBayes(fit2)

#plot the mean-variance trend
plotSA(fit2, main = "Final model: Mean-Variance trend")

Past versions of "Differential Expression Analysis-2.png"

#look at the summary of your results
##this tells you the number of DEGs in each condition
results_summary <- decideTests(fit2, adjust.method = "BH", p.value = 0.05)
summary(results_summary)

       V.D24T V.D24R V.D144R
Down     4723   3593     359
NotSig   5076   7151   13810
Up       4520   3575     150

#        V.D24 V.D24r V.D144r
# Down    4723   3593     359
# NotSig  5076   7151   13810
# Up      4520   3575     150

vennDiagram(object = results_summary, include = c("up", "down"))

Past versions of "Differential Expression Analysis-3.png"

##Create Toptables of my DEGs

# Generate Top Table for Specific Comparisons

Toptable_V.D24T <- topTable(fit = fit2, coef = "V.D24T", number = nrow(x), adjust.method = "BH", p.value = 1, sort.by = "none")
#write.csv(Toptable_V.D24T, "data/new/DEGs/Toptable_V.D24T.csv")

Toptable_V.D24R <- topTable(fit = fit2, coef = "V.D24R", number = nrow(x), adjust.method = "BH", p.value = 1, sort.by = "none")
#write.csv(Toptable_V.D24R, "data/new/DEGs/Toptable_V.D24R.csv")

Toptable_V.D144R <- topTable(fit = fit2, coef = "V.D144R", number = nrow(x), adjust.method = "BH", p.value = 1, sort.by = "none")
#write.csv(Toptable_V.D144R, "data/new/DEGs/Toptable_V.D144R.csv")

#save all of these toptables as R objects
# saveRDS(list(
#   V.D24T = Toptable_V.D24T,
#   V.D24R = Toptable_V.D24R,
#   V.D144R = Toptable_V.D144R
# ), file = "data/new/Toptable_list.RDS")

Toptable_list <- readRDS("data/new/Toptable_list.RDS")

##Top 5 DEGs after DE Analysis per Condition

#use your three toptables so I can pull out top 5 genes from each based on adj. p val

top5_D24T <- Toptable_V.D24T[order(Toptable_V.D24T$adj.P.Val), ][1:5,] %>% 
  rownames_to_column(., var = "Entrez_ID")
top5_D24R <- Toptable_V.D24R[order(Toptable_V.D24R$adj.P.Val), ][1:5,] %>% 
  rownames_to_column(., var = "Entrez_ID")
top5_D144R <- Toptable_V.D144R[order(Toptable_V.D144R$adj.P.Val), ][1:5,] %>% 
  rownames_to_column(., var = "Entrez_ID")

#now that I've pulled the top 5 DEGs from each, make a list to pull them from my log2cpm data

boxplot1 <- read.csv("data/new/filcpm_final_matrix.csv") %>% 
  as.data.frame()

#Define gene list
#these are the top 5 genes pulled from my toptables
top5_D24T_geneslist <- c(top5_D24T$Entrez_ID)
top5_D24R_geneslist <- c(top5_D24R$Entrez_ID)
top5_D144R_geneslist <- c(top5_D144R$Entrez_ID)
#Add more gene symbols as needed or add more categories

#now pull these from my log2cpm matrix
top5_D24T_genes <- boxplot1[boxplot1$Entrez_ID %in% top5_D24T_geneslist,]
dim(top5_D24T_genes)

[1]  5 44

#5 genes in 44 cols
top5_D24R_genes <- boxplot1[boxplot1$Entrez_ID %in% top5_D24R_geneslist,]
dim(top5_D24R_genes)

[1]  5 44

#5 genes in 44 cols
top5_D144R_genes <- boxplot1[boxplot1$Entrez_ID %in% top5_D144R_geneslist,]
dim(top5_D144R_genes)

[1]  5 44

#5 genes in 44 cols

#Now put in the function I want to use to generate boxplots of genes
#####D24T#####
process_top5_D24T <- function(gene) {
  gene_data <- top5_D24T_genes %>% filter(Entrez_ID == gene)
  long_data <- gene_data %>%
    pivot_longer(cols = -c(Entrez_ID, SYMBOL), names_to = "Sample", values_to = "log2cpm") %>%
    mutate(
      Drug = case_when(
        grepl("DOX", Sample) ~ "DOX",
        grepl("DMSO", Sample) ~ "DMSO",
        TRUE ~ NA_character_
      ),
      Timepoint = case_when(
        grepl("_24T_", Sample) ~ "24T",
        grepl("_24R_", Sample) ~ "24R",
        grepl("_144R_", Sample) ~ "144R",
        TRUE ~ NA_character_
      ),
      Indv = case_when(
        grepl("Ind1$", Sample) ~ "1",
        grepl("Ind2$", Sample) ~ "2",
        grepl("Ind3$", Sample) ~ "3",
        grepl("Ind4$", Sample) ~ "4",
        grepl("Ind5$", Sample) ~ "5",
        grepl("Ind6$", Sample) ~ "6",
        grepl("Ind6R$", Sample) ~ "6R",
        TRUE ~ NA_character_
      ),
      Condition = paste(Drug, Timepoint, sep = "_")
    )
  long_data$Condition <- factor(
    long_data$Condition,
    levels = c(
      "DOX_24T", 
      "DMSO_24T", 
      "DOX_24R", 
      "DMSO_24R", 
      "DOX_144R", 
      "DMSO_144R"
    )
  )
  return(long_data)
}


#Generate Boxplots from the above function using our gene list above
for (gene in top5_D24T_geneslist) {
  gene_data <- process_top5_D24T(gene)
  p <- ggplot(gene_data, aes(x = Condition, y = log2cpm, fill = Drug)) +
    geom_boxplot(outlier.shape = NA) +
    geom_point(aes(color = Indv), size = 2, alpha = 0.5, position = position_jitter(width = -1, height = 0)) +
    scale_fill_manual(values = c("DOX" = "#499FBD", "DMSO" = "#BBBBBC")) +
    ggtitle(paste("Log2cpm Expression of", gene, "top 5 DEGs D24T")) +
    labs(x = "Treatment", y = "log2cpm") +
    theme_bw() +
    theme(
      plot.title = element_text(size = rel(2), hjust = 0.5),
      axis.title = element_text(size = 15, color = "black"),
      axis.text.x = element_text(size = 10, color = "black", angle = 90, hjust = 1)
    )
  print(p)
}

Past versions of "Top 5 DEGs Limma-1.png"

Past versions of "Top 5 DEGs Limma-2.png"

Past versions of "Top 5 DEGs Limma-3.png"

Past versions of "Top 5 DEGs Limma-4.png"

Past versions of "Top 5 DEGs Limma-5.png"

#####D24R#####
process_top5_D24R <- function(gene) {
  gene_data <- top5_D24R_genes %>% filter(Entrez_ID == gene)
  long_data <- gene_data %>%
    pivot_longer(cols = -c(Entrez_ID, SYMBOL), names_to = "Sample", values_to = "log2cpm") %>%
    mutate(
      Drug = case_when(
        grepl("DOX", Sample) ~ "DOX",
        grepl("DMSO", Sample) ~ "DMSO",
        TRUE ~ NA_character_
      ),
      Timepoint = case_when(
        grepl("_24T_", Sample) ~ "24T",
        grepl("_24R_", Sample) ~ "24R",
        grepl("_144R_", Sample) ~ "144R",
        TRUE ~ NA_character_
      ),
      Indv = case_when(
        grepl("Ind1$", Sample) ~ "1",
        grepl("Ind2$", Sample) ~ "2",
        grepl("Ind3$", Sample) ~ "3",
        grepl("Ind4$", Sample) ~ "4",
        grepl("Ind5$", Sample) ~ "5",
        grepl("Ind6$", Sample) ~ "6",
        grepl("Ind6R$", Sample) ~ "6R",
        TRUE ~ NA_character_
      ),
      Condition = paste(Drug, Timepoint, sep = "_")
    )
  long_data$Condition <- factor(
    long_data$Condition,
    levels = c(
      "DOX_24T", 
      "DMSO_24T", 
      "DOX_24R", 
      "DMSO_24R", 
      "DOX_144R", 
      "DMSO_144R"
    )
  )
  return(long_data)
}


#Generate Boxplots from the above function using our gene list above
for (gene in top5_D24R_geneslist) {
  gene_data <- process_top5_D24R(gene)
  p <- ggplot(gene_data, aes(x = Condition, y = log2cpm, fill = Drug)) +
    geom_boxplot(outlier.shape = NA) +
    geom_point(aes(color = Indv), size = 2, alpha = 0.5, position = position_jitter(width = -1, height = 0)) +
    scale_fill_manual(values = c("DOX" = "#499FBD", "DMSO" = "#BBBBBC")) +
    ggtitle(paste("Log2cpm Expression of", gene, "top 5 DEGs D24R")) +
    labs(x = "Treatment", y = "log2cpm") +
    theme_bw() +
    theme(
      plot.title = element_text(size = rel(2), hjust = 0.5),
      axis.title = element_text(size = 15, color = "black"),
      axis.text.x = element_text(size = 10, color = "black", angle = 90, hjust = 1)
    )
  print(p)
}

Past versions of "Top 5 DEGs Limma-6.png"

Past versions of "Top 5 DEGs Limma-7.png"

Past versions of "Top 5 DEGs Limma-8.png"

Past versions of "Top 5 DEGs Limma-9.png"

Past versions of "Top 5 DEGs Limma-10.png"

#####D144R#####
process_top5_D144R <- function(gene) {
  gene_data <- top5_D144R_genes %>% filter(Entrez_ID == gene)
  long_data <- gene_data %>%
    pivot_longer(cols = -c(Entrez_ID, SYMBOL), names_to = "Sample", values_to = "log2cpm") %>%
    mutate(
      Drug = case_when(
        grepl("DOX", Sample) ~ "DOX",
        grepl("DMSO", Sample) ~ "DMSO",
        TRUE ~ NA_character_
      ),
      Timepoint = case_when(
        grepl("_24T_", Sample) ~ "24T",
        grepl("_24R_", Sample) ~ "24R",
        grepl("_144R_", Sample) ~ "144R",
        TRUE ~ NA_character_
      ),
      Indv = case_when(
        grepl("Ind1$", Sample) ~ "1",
        grepl("Ind2$", Sample) ~ "2",
        grepl("Ind3$", Sample) ~ "3",
        grepl("Ind4$", Sample) ~ "4",
        grepl("Ind5$", Sample) ~ "5",
        grepl("Ind6$", Sample) ~ "6",
        grepl("Ind6R$", Sample) ~ "6R",
        TRUE ~ NA_character_
      ),
      Condition = paste(Drug, Timepoint, sep = "_")
    )
  long_data$Condition <- factor(
    long_data$Condition,
    levels = c(
      "DOX_24T", 
      "DMSO_24T", 
      "DOX_24R", 
      "DMSO_24R", 
      "DOX_144R", 
      "DMSO_144R"
    )
  )
  return(long_data)
}


#Generate Boxplots from the above function using our gene list above
for (gene in top5_D144R_geneslist) {
  gene_data <- process_top5_D144R(gene)
  p <- ggplot(gene_data, aes(x = Condition, y = log2cpm, fill = Drug)) +
    geom_boxplot(outlier.shape = NA) +
    geom_point(aes(color = Indv), size = 2, alpha = 0.5, position = position_jitter(width = -1, height = 0)) +
    scale_fill_manual(values = c("DOX" = "#499FBD", "DMSO" = "#BBBBBC")) +
    ggtitle(paste("Log2cpm Expression of", gene, "top 5 DEGs D144R")) +
    labs(x = "Treatment", y = "log2cpm") +
    theme_bw() +
    theme(
      plot.title = element_text(size = rel(2), hjust = 0.5),
      axis.title = element_text(size = 15, color = "black"),
      axis.text.x = element_text(size = 10, color = "black", angle = 90, hjust = 1)
    )
  print(p)
}

Past versions of "Top 5 DEGs Limma-11.png"

Past versions of "Top 5 DEGs Limma-12.png"

Past versions of "Top 5 DEGs Limma-13.png"

Past versions of "Top 5 DEGs Limma-14.png"

Past versions of "Top 5 DEGs Limma-15.png"

##LogFC Boxplots for All Conditions

# Toptable_V.D24T

# Load DEGs Data
DOX_24T <- read.csv("data/new/DEGs/Toptable_V.D24T.csv")
DOX_24R <- read.csv("data/new/DEGs/Toptable_V.D24R.csv")
DOX_144R <- read.csv("data/new/DEGs/Toptable_V.D144R.csv")

#make a list of all of the genes in this set so I can plot the logFC in other sets
D24T_DEGs <- DOX_24T$Entrez_ID[DOX_24T$adj.P.Val < 0.05]
length(D24T_DEGs)

[1] 9243

#9243 genes in length after adj. p value cutoff
#if I did p value only - 9658 genes

#now that I have a list of my DEGs from D24T - pull these genes out from the other DEG lists
# D24R_DEGs <- Toptable_V.D24R[Toptable_V.D24R$adj.P.Val < 0.05,]

# D24R_DEGs_D24T <- D24R_DEGs %>%  
#   rownames_to_column(., var = "Entrez_ID")

# D24R_DEGs_D24T <- D24R_DEGs_D24T[D24R_DEGs_D24T$Entrez_ID %in% D24T_DEGs,]
# dim(D24R_DEGs_D24T)
#4831 genes in common here after adj p value cutoff

# D144R_DEGs <- Toptable_V.D144R[Toptable_V.D144R$adj.P.Val < 0.05,]
# 
# D144R_DEGs_D24T <- D144R_DEGs %>%  
#   rownames_to_column(., var = "Entrez_ID")
# 
# D144R_DEGs_D24T <- D144R_DEGs_D24T[D144R_DEGs_D24T$Entrez_ID %in% D24T_DEGs,]
# dim(D144R_DEGs_D24T)
#322 genes in common after adj p value cutoff

#now I want to plot the logFC of these

#ignore the above for now, just plot those full gene sets in logFC
D24T_DEGs <- DOX_24T$Entrez_ID[DOX_24T$adj.P.Val < 0.05]
length(D24T_DEGs)

[1] 9243

D24R_DEGs <- DOX_24R$Entrez_ID[DOX_24R$adj.P.Val < 0.05]
length(D24R_DEGs)

[1] 7168

D144R_DEGs <- DOX_144R$Entrez_ID[DOX_144R$adj.P.Val < 0.05]
length(D144R_DEGs)

[1] 509

#now that I have the full list of genes, I want to plot the logFC across conditions
#to do this - make a combined toptable

# Toptable_list

#plot the set of genes across conditions
# 
# logFC_long_allsets <- imap_dfr(Toptable_list, function(tbl, condition) {
#   tbl %>% 
#     rownames_to_column(var = "Entrez_ID") %>% 
#     dplyr::select(Entrez_ID, logFC) %>% 
#     mutate(Condition = condition)
# })
# 
# 
# ggplot(logFC_long_allsets, aes(x = Condition, y = logFC, group = Entrez_ID, color = Entrez_ID)) +
#   geom_boxplot() +
#   geom_point(size = 2) +
#   labs(
#     title = "logFC of Genes Across Conditions",
#     x = "Condition",
#     y = "logFC"
#   ) +
#   theme_minimal()


####Try this####

#Combine the toptables I have from pairwise analysis into a single dataframe
d24_toptable_dxr <- Toptable_V.D24T %>% 
  rownames_to_column(var = "Entrez_ID") %>% 
  mutate(Time = "24")

d24r_toptable_dxr <- Toptable_V.D24R %>% 
  rownames_to_column(var = "Entrez_ID") %>% 
  mutate(Time = "24R")

d144r_toptable_dxr <- Toptable_V.D144R %>% 
  rownames_to_column(var = "Entrez_ID") %>% 
  mutate(Time = "144R")

combined_toptables_dxr <- bind_rows(
  d24_toptable_dxr,
  d24r_toptable_dxr,
  d144r_toptable_dxr)

#Filter the data based on each motif
filt_toptable_dxr <- combined_toptables_dxr %>% 
  dplyr::filter(Entrez_ID  %in% D24T_DEGs) %>% 
  mutate(absFC = abs(logFC)) %>% 
  mutate(time = factor(Time, levels = c("24", "24R", "144R"), labels = c("24hr", "24hr Recovery", "144hr Recovery"))) %>% 
  ggplot(., aes(x = time, y = logFC))+
  geom_boxplot(aes(fill = time))+
  scale_fill_manual(values = time_col)+
  guides(fill = guide_legend(title = "Time"))+
  theme_bw()+
  xlab(" ")+
  ylab("logFC")+
  theme_bw()+
  ggtitle("LogFC for D24T DEGs Across Conditions")+
  theme(plot.title = element_text(size = rel(1.5), hjust = 0.5),
        axis.title = element_text(size = 15, colour = "black"),
        strip.background = element_rect(fill = "#CAD899"),
        axis.text.x = element_text(size = 8, colour = "white", angle = 15),
        strip.text.x = element_text(size = 12, colour = "black", face = "bold"))

#now do the same with the other conditions

filt_toptable_dxr_24r <- combined_toptables_dxr %>% 
  dplyr::filter(Entrez_ID  %in% D24R_DEGs) %>% 
  mutate(absFC = abs(logFC)) %>% 
  mutate(time = factor(Time, levels = c("24", "24R", "144R"), labels = c("24hr", "24hr Recovery", "144hr Recovery"))) %>% 
  ggplot(., aes(x = time, y = logFC))+
  geom_boxplot(aes(fill = time))+
  scale_fill_manual(values = time_col)+
  guides(fill = guide_legend(title = "Time"))+
  theme_bw()+
  xlab(" ")+
  ylab("logFC")+
  theme_bw()+
  ggtitle("LogFC for D24R DEGs Across Conditions")+
  theme(plot.title = element_text(size = rel(1.5), hjust = 0.5),
        axis.title = element_text(size = 15, colour = "black"),
        strip.background = element_rect(fill = "#CAD899"),
        axis.text.x = element_text(size = 8, colour = "white", angle = 15),
        strip.text.x = element_text(size = 12, colour = "black", face = "bold"))

#D144R

filt_toptable_dxr_144r <- combined_toptables_dxr %>% 
  dplyr::filter(Entrez_ID  %in% D144R_DEGs) %>% 
  mutate(absFC = abs(logFC)) %>% 
  mutate(time = factor(Time, levels = c("24", "24R", "144R"), labels = c("24hr", "24hr Recovery", "144hr Recovery"))) %>% 
  ggplot(., aes(x = time, y = logFC))+
  geom_boxplot(aes(fill = time))+
  scale_fill_manual(values = time_col)+
  guides(fill = guide_legend(title = "Time"))+
  theme_bw()+
  xlab(" ")+
  ylab("logFC")+
  theme_bw()+
  ggtitle("LogFC for D144R DEGs Across Conditions")+
  theme(plot.title = element_text(size = rel(1.5), hjust = 0.5),
        axis.title = element_text(size = 15, colour = "black"),
        strip.background = element_rect(fill = "#CAD899"),
        axis.text.x = element_text(size = 8, colour = "white", angle = 15),
        strip.text.x = element_text(size = 12, colour = "black", face = "bold"))

#now put the names of these graphs to print them
filt_toptable_dxr

Past versions of "DOX24T DEGs logFC all Conditions-1.png"

filt_toptable_dxr_24r

Past versions of "DOX24T DEGs logFC all Conditions-2.png"

filt_toptable_dxr_144r

Past versions of "DOX24T DEGs logFC all Conditions-3.png"

#Perform RUVs Correction ##Set Up Data and Perform PCAs on Un-normalized Data

filt_gene_list <- rownames(filcpm_matrix)
#14319 genes as usual

#in order to make this match with annot later down the line, change the col names for counts_raw_matrix to match final_sample_names in annot

#i'll also want to make sure I keep the replicate for this set

colnames(counts_raw_matrix) <- Metadata$Final_sample_name

RUV_filt_counts <- counts_raw_matrix %>% 
  as.data.frame() %>% 
  dplyr::filter(., row.names(.)%in% filt_gene_list)

#add in the annotation files
ind_num <- readRDS("data/new/ind_num.RDS")
annot <- read.csv("data/new/Metadata.csv")

#  counts need to be integer values and in a numeric matrix
# note: the log transformation needs to be accounted for in the isLog argument in RUVs function.
counts <- as.matrix(RUV_filt_counts)
#saveRDS(counts, "data/new/RUV/filt_counts_matrix.RDS")

# Create a DataFrame for the phenoData
phenoData <- DataFrame(annot)

# Now create the RangedSummarizedExperiment necessary for RUVs input
# looks like it did need both the phenodata and the counts.
set <- SummarizedExperiment(assays =  counts, metadata = phenoData)

# Generate a background matrix
# The column "Cond" holds the comparisons that you actually want to make. DOX_24, DMSO_24,5FU_24, DOX_3,etc.
scIdx <-RUVSeq::makeGroups(phenoData$Condition)
scIdx

     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]    6   12   18   24   30   36   42
[2,]    4   10   16   22   28   34   40
[3,]    2    8   14   20   26   32   38
[4,]    5   11   17   23   29   35   41
[5,]    3    9   15   21   27   33   39
[6,]    1    7   13   19   25   31   37

#now I've made all of the data I need for this - they are located in each section for k values

#DO NOT USE THESE COUNTS FOR LINEAR MODELING

#colors for all of the plots
fill_col_ind <- c("#66C2A5", "#FC8D62", "#1F78B4", "#E78AC3", "#A6D854", "#FFD92A", "#8B3E9B")

fill_col_ind_dark <- c("#003F5C", "#45AE91",  "#58508D", "#BC4099", "#8B3E9B", "#FF6361", "#FF2362")

fill_col_tx <- c("#63666D", "#499FBD", "#DCACED")

fill_col_txtime <- c("#003F5C", "#45AE91",  "#58508D", "#BC4099", "#8B3E9B", "#FF6361", "#FF2362", "#A6D854", "#FC8D62")

# before ruv (counts PCA)
prcomp_res_counts <- prcomp(t(counts), scale. = FALSE, center = TRUE)
annot_prcomp_res <- prcomp_res_counts$x %>% cbind(., annot)

group_2 <- annot$Condition

#now plot my PCA for filtered counts
####PC1/PC2####
ggplot2::autoplot(prcomp_res_counts, 
                  data = annot, 
                  colour = "Condition", 
                  shape = "Time", 
                  size =4, 
                  x=1, 
                  y=2) +
  ggrepel::geom_text_repel(label=ind_num, 
                           vjust = -0.5, 
                           max.overlaps = 30) +
  scale_color_manual(values=cond_col) +
  ggtitle(expression("PCA of Filtered Counts")) +
  theme_bw()

Past versions of "RUVs Correction Data Setup - no RUV yet-1.png"

####PC2/PC3####
ggplot2::autoplot(prcomp_res_counts, 
                  data = annot, 
                  colour = "Condition", 
                  shape = "Time", 
                  size =4, 
                  x=2, 
                  y=3) +
  ggrepel::geom_text_repel(label=ind_num, 
                           vjust = -0.5, 
                           max.overlaps = 30) +
  scale_color_manual(values=cond_col) +
  ggtitle(expression("PCA of Filtered Counts")) +
  theme_bw()

Past versions of "RUVs Correction Data Setup - no RUV yet-2.png"

#go ahead and plot PCA of log2cpm to compare (somewhat) later since the norm counts output isn't possible with these data since they don't undergo correction

prcomp_res_cpm <- prcomp(t(filcpm_matrix %>% as.matrix()), center =  TRUE)

####PC1/PC2####
ggplot2::autoplot(prcomp_res_cpm, 
                  data = annot, 
                  colour = "Condition", 
                  shape = "Time", 
                  size =4, 
                  x=1, 
                  y=2) +
  ggrepel::geom_text_repel(label=ind_num, 
                           vjust = -0.5, 
                           max.overlaps = 30) +
  scale_color_manual(values=cond_col) +
  ggtitle(expression("PCA of log2cpm no RUVs")) +
  theme_bw()

Past versions of "RUVs Correction Data Setup - no RUV yet-3.png"

####PC2/PC3####
ggplot2::autoplot(prcomp_res_cpm, 
                  data = annot, 
                  colour = "Condition", 
                  shape = "Time", 
                  size =4, 
                  x=2, 
                  y=3) +
  ggrepel::geom_text_repel(label=ind_num, 
                           vjust = -0.5, 
                           max.overlaps = 30) +
  scale_color_manual(values=cond_col) +
  ggtitle(expression("PCA of log2cpm no RUVs")) +
  theme_bw()

Past versions of "RUVs Correction Data Setup - no RUV yet-4.png"

####new PCA plots no correction####
#PCA plots for each value of k attached in each section

#####Now start performing RUV 1-3###

##RUVs with k = 1

#Apply RUVs function from RUVSeq
#"k" will be iteratively adjusted over time depending on your PCA.
set1 <- RUVSeq::RUVs(x = counts, k =1, scIdx = scIdx, isLog = FALSE)

#get the ruv weights to put into the linear model. n weights = k.
#k=1
RUV_df1 <- set1$W %>% as.data.frame()
RUV_df1$Names <- rownames(RUV_df1)

#Check that the names match
#k=1
RUV_df_rm1 <- RUV_df1[RUV_df1$Names %in% annot$Final_sample_name, ] 
RUV_1 <-  RUV_df_rm1$W_1

#saveRDS(RUV_df_rm1, "data/new/RUV_df_rm1.RDS")
#saveRDS(RUV_1, "data/new/RUV_1.RDS")

#PCA checks
#k=1
prcomp_res_1 <- prcomp(t(set1$normalizedCounts), scale. = FALSE, center = TRUE)
annot_prcomp_res_1 <- prcomp_res_1$x %>% cbind(., annot)

ggplot2::autoplot(prcomp_res_1, data = annot, colour = "Condition", shape = "Time", size =4, x = 1, y = 2)+
  theme_bw()+
  scale_color_manual(values=cond_col)+
  ggrepel::geom_text_repel(label= ind_num, vjust = -0.5, max.overlaps = 30)+
  ggtitle("RUVs Correction k=1 NormCounts")

Past versions of "RUVs k1-1.png"

ggplot2::autoplot(prcomp_res_1, data = annot, colour = "Condition", shape = "Time", size =4, x=2, y=3)+
  theme_bw()+
  scale_color_manual(values=cond_col)+
  ggrepel::geom_text_repel(label= ind_num, vjust = -0.5, max.overlaps = 30)+
  ggtitle("RUVs Correction k=1 NormCounts")

Past versions of "RUVs k1-2.png"

#also try by converting these values to log2cpm

RUV_df_rm1_cpm <- cpm(set1$normalizedCounts, log = TRUE)

prcomp_res_1_cpm <- prcomp(t(RUV_df_rm1_cpm), scale. = FALSE, center = TRUE)
annot_prcomp_res_1_cpm <- prcomp_res_1_cpm$x %>% cbind(., annot)

##PC1/2
ggplot2::autoplot(prcomp_res_1_cpm, 
                  data = annot, 
                  colour = "Condition", 
                  shape = "Time", 
                  size = 4,
                  x=1,
                  y=2)+
  theme_bw()+
  scale_color_manual(values=cond_col)+
  ggrepel::geom_text_repel(label= ind_num, 
                           vjust = -0.5, 
                           max.overlaps = 30)+
  ggtitle("RUVs Correction k=1 log2cpm")

Past versions of "RUVs k1-3.png"

###PC2/3
ggplot2::autoplot(prcomp_res_1_cpm, 
                  data = annot, 
                  colour = "Condition", 
                  shape = "Time", 
                  size = 4, 
                  x=2, 
                  y=3)+
  theme_bw()+
  scale_color_manual(values=cond_col)+
  ggrepel::geom_text_repel(label= ind_num, 
                           vjust = -0.5, 
                           max.overlaps = 30)+
  ggtitle("RUVs Correction k=1 log2cpm")

Past versions of "RUVs k1-4.png"

##RUVs with k = 2

#Apply RUVs function from RUVSeq
#"k" will be iteratively adjusted over time depending on your PCA.
set2 <- RUVSeq::RUVs(x = counts, k =2, scIdx = scIdx, isLog = FALSE)

#get the ruv weights to put into the linear model. n weights = k.
#k=2
RUV_df2 <- set2$W %>% as.data.frame()
RUV_df2$Names <- rownames(RUV_df2)

#Check that the names match
#k=2
RUV_df_rm2 <- RUV_df2[RUV_df2$Names %in% annot$Final_sample_name, ] 
RUV_2 <-  RUV_df_rm2$W_2

#PCA checks
#k=2
prcomp_res_2 <- prcomp(t(set2$normalizedCounts), scale. = FALSE, center = TRUE)
annot_prcomp_res_2 <- prcomp_res_2$x %>% cbind(., annot)

##PC1/2
ggplot2::autoplot(prcomp_res_2, 
                  data = annot, 
                  colour = "Condition", 
                  shape = "Time", 
                  size =4, 
                  x=1,
                  y=2)+
  theme_bw()+
  scale_color_manual(values=cond_col)+
  ggrepel::geom_text_repel(label= ind_num, 
                           vjust = -0.5, 
                           max.overlaps = 30)+
  ggtitle("RUVs Correction k=2 NormCounts")

Warning: ggrepel: 1 unlabeled data points (too many overlaps). Consider
increasing max.overlaps

Past versions of "RUVs k2-1.png"

###PC2/3
ggplot2::autoplot(prcomp_res_2, 
                  data = annot, 
                  colour = "Condition", 
                  shape = "Time", 
                  size =4, 
                  x=2, 
                  y=3)+
  theme_bw()+
  scale_color_manual(values=cond_col)+
  ggrepel::geom_text_repel(label= ind_num, 
                           vjust = -0.5, 
                           max.overlaps = 30)+
  ggtitle("RUVs Correction k=2 NormCounts")

Warning: ggrepel: 13 unlabeled data points (too many overlaps). Consider
increasing max.overlaps

Past versions of "RUVs k2-2.png"

#now convert this to log2cpm using the normalized counts from set2

RUV_df_rm2_cpm <- cpm(set2$normalizedCounts, log = TRUE)

prcomp_res_2_cpm <- prcomp(t(RUV_df_rm2_cpm), scale. = FALSE, center = TRUE)
annot_prcomp_res_2_cpm <- prcomp_res_2_cpm$x %>% cbind(., annot)

ggplot2::autoplot(prcomp_res_2_cpm, 
                  data = annot, 
                  colour = "Condition", 
                  shape = "Time", 
                  size =4, 
                  x=1, 
                  y=2)+
  theme_bw()+
  scale_color_manual(values=cond_col)+
  ggrepel::geom_text_repel(label= ind_num, 
                           vjust = -0.5, 
                           max.overlaps = 30)+
  ggtitle("RUVs Correction k=2 log2cpm")

Past versions of "RUVs k2-3.png"

###PC2/3
ggplot2::autoplot(prcomp_res_2_cpm, 
                  data = annot, 
                  colour = "Condition", 
                  shape = "Time", 
                  size =4, 
                  x=2, 
                  y=3)+
  theme_bw()+
  scale_color_manual(values=cond_col)+
  ggrepel::geom_text_repel(label= ind_num, 
                           vjust = -0.5, 
                           max.overlaps = 30)+
  ggtitle("RUVs Correction k=2 log2cpm")

Past versions of "RUVs k2-4.png"

#RUVs with k = 3

#Apply RUVs function from RUVSeq
#"k" will be iteratively adjusted over time depending on your PCA.
set3 <- RUVSeq::RUVs(x = counts, k =3, scIdx = scIdx, isLog = FALSE)

#get the ruv weights to put into the linear model. n weights = k.
#k=3
RUV_df3 <- set3$W %>% as.data.frame()
RUV_df3$Names <- rownames(RUV_df3)

#Check that the names match
#k=3
RUV_df_rm3 <- RUV_df3[RUV_df3$Names %in% annot$Final_sample_name, ] 
RUV_3 <-  RUV_df_rm3$W_3

#PCA checks
#k=3
prcomp_res_3 <- prcomp(t(set3$normalizedCounts), scale. = FALSE, center = TRUE)
annot_prcomp_res_3 <- prcomp_res_3$x %>% cbind(., annot)

###PC1/2
ggplot2::autoplot(prcomp_res_3, 
                  data = annot, 
                  colour = "Condition", 
                  shape = "Time", 
                  size =4, 
                  x=1, 
                  y=2)+
  theme_bw()+
  scale_color_manual(values=cond_col)+
  ggrepel::geom_text_repel(label= ind_num, 
                           vjust = -0.5, 
                           max.overlaps = 30)+
  ggtitle("RUVs Correction k=3 NormCounts")

Past versions of "RUVs k3-1.png"

###PC2/3
ggplot2::autoplot(prcomp_res_3, 
                  data = annot, 
                  colour = "Condition", 
                  shape = "Time", 
                  size =4, 
                  x=2, 
                  y=3)+
  theme_bw()+
  scale_color_manual(values=cond_col)+
  ggrepel::geom_text_repel(label= ind_num, 
                           vjust = -0.5, 
                           max.overlaps = 30)+
  ggtitle("RUVs Correction k=3 NormCounts")

Past versions of "RUVs k3-2.png"

#now convert this to log2cpm using the normalized counts from set3

RUV_df_rm3_cpm <- cpm(set3$normalizedCounts, log = TRUE)

prcomp_res_3_cpm <- prcomp(t(RUV_df_rm3_cpm), scale. = FALSE, center = TRUE)
annot_prcomp_res_3_cpm <- prcomp_res_3_cpm$x %>% cbind(., annot)

##PC1/2
ggplot2::autoplot(prcomp_res_3_cpm, 
                  data = annot, 
                  colour = "Condition", 
                  shape = "Time", 
                  size =4, 
                  x=1, 
                  y=2)+
  theme_bw()+
  scale_color_manual(values=cond_col)+
  ggrepel::geom_text_repel(label= ind_num, 
                           vjust = -0.5, 
                           max.overlaps = 30)+
  ggtitle("RUVs Correction k=3 log2cpm")

Past versions of "RUVs k3-3.png"

###PC2/3
ggplot2::autoplot(prcomp_res_3_cpm, 
                  data = annot, 
                  colour = "Condition", 
                  shape = "Time", 
                  size =4, 
                  x=2, 
                  y=3)+
  theme_bw()+
  scale_color_manual(values=cond_col)+
  ggrepel::geom_text_repel(label= ind_num, 
                           vjust = -0.5, 
                           max.overlaps = 30)+
  ggtitle("RUVs Correction k=3 log2cpm")

Past versions of "RUVs k3-4.png"

#Plot Spearman Heatmaps for RUVs Normalized Counts Data

#Now that I've put together the PCA plots for both normalized counts and log2cpm
#I want to make these into heatmaps

####RUVs k=1-3
#check to make sure that the column names are correct
dim(RUV_filt_counts)

[1] 14319    42

dim(set1$normalizedCounts)

[1] 14319    42

dim(set2$normalizedCounts)

[1] 14319    42

dim(set3$normalizedCounts)

[1] 14319    42

#take the normalized counts from k=1 and put together a dataframe with the correct columns
normcounts_k0 <- RUV_filt_counts
normcounts_k1 <- set1$normalizedCounts %>% as.data.frame()
#do the same with k=2 and k=3
normcounts_k2 <- set2$normalizedCounts %>% as.data.frame()
normcounts_k3 <- set3$normalizedCounts %>% as.data.frame()

#do the same with the log2cpm conversion
cpm_k0 <- cpm(normcounts_k0, log = TRUE) %>% as.data.frame()
cpm_k1 <- cpm(set1$normalizedCounts, log = TRUE) %>% as.data.frame()
cpm_k2 <- cpm(set2$normalizedCounts, log = TRUE) %>% as.data.frame()
cpm_k3 <- cpm(set3$normalizedCounts, log = TRUE) %>% as.data.frame()

#compute the correlation matrices for RUVs 1-3 with normalized counts
#k=0
cor_matrix_spmn_k0 <- cor(normcounts_k0,
                          y = NULL,
                          use = "everything", 
                          method = "spearman")
#k=1
cor_matrix_spmn_k1 <- cor(normcounts_k1,
                           y = NULL,
                           use = "everything",
                           method = "spearman")

#k=2
cor_matrix_spmn_k2 <- cor(normcounts_k2,
                          y = NULL,
                          use = "everything",
                          method = "spearman")

#k=3
cor_matrix_spmn_k3 <- cor(normcounts_k3,
                          y = NULL,
                          use = "everything",
                          method = "spearman")

#Do the same with the log2cpm converted versions
#k=0
cor_matrix_spmn_k0_cpm <- cor(cpm_k0,
                           y = NULL,
                           use = "everything",
                           method = "spearman")

#k=1
cor_matrix_spmn_k1_cpm <- cor(cpm_k1,
                           y = NULL,
                           use = "everything",
                           method = "spearman")

#k=2
cor_matrix_spmn_k2_cpm <- cor(cpm_k2,
                          y = NULL,
                          use = "everything",
                          method = "spearman")

#k=3
cor_matrix_spmn_k3_cpm <- cor(cpm_k3,
                          y = NULL,
                          use = "everything",
                          method = "spearman")

#extract metadata columns
Individual <- as.character(Metadata$Ind)
Time <- as.character(Metadata$Time)
Treatment <- as.character(Metadata$Drug)

# Define color palettes for annotations
annot_col_cor = list(Treatment = c("DOX" = "#499FBD", 
                               "DMSO" = "#BBBBBC"),
                     Individuals = c("1" = "#003F5C", 
                                     "2" = "#45AE91", 
                                     "3" = "#58209D", 
                                     "4" = "#8B3E9B", 
                                     "5" = "#FF6361", 
                                     "6" = "#BC4169", 
                                     "6R" = "#FF2362"),
                     Timepoints = c("24T" = "#238B45",
                                    "24R" = "#74C476", 
                                    "144R" = "#C7E9C0"))

tx_colors <- c("DOX" = "#499FBD", 
                 "DMSO" = "#BBBBBC")
ind_colors <- c("1" = "red",
                "2" = "orange",
                "3" = "yellow",
                "4" = "green",
                "5" = "blue",
                "6" = "violet",
                "6R" = "purple")
time_colors <- c("24T" = "#238B45",
                 "24R" = "#74C476",
                 "144R" = "#C7E9C0")

# Create annotations
top_annotation <- HeatmapAnnotation(
  Individual = Individual, 
  Time = Time,
  Treatment = Treatment,
  col = list(
    Individual = ind_colors, 
    Time = time_colors,
    Treatment = tx_colors
  )
)

####ANNOTATED HEATMAPS####
###Spearman Heatmap k=0 ####
heatmap_spmn_k0 <- Heatmap(cor_matrix_spmn_k0,
                           name = "Spearman",
                           top_annotation = top_annotation,
                           show_row_names = TRUE,
                           show_column_names = TRUE,
                           cluster_rows = TRUE,
                           cluster_columns = TRUE,
                           border = TRUE, 
                           column_title = "Filtered Counts no RUVs")

# Draw the heatmap k=0
draw(heatmap_spmn_k0)

Past versions of "RUVs HM Spearman-1.png"

####Spearman Heatmap k=1 ####
heatmap_spmn_k1 <- Heatmap(cor_matrix_spmn_k1,
                           name = "Spearman",
                           top_annotation = top_annotation,
                           show_row_names = TRUE,
                           show_column_names = TRUE,
                           cluster_rows = TRUE,
                           cluster_columns = TRUE,
                           border = TRUE, 
                           column_title = "Normalized Counts k=1")

# Draw the heatmap k=1
draw(heatmap_spmn_k1)

Past versions of "RUVs HM Spearman-2.png"

####Spearman Heatmap k=2####
heatmap_spmn_k2 <- Heatmap(cor_matrix_spmn_k2,
                           name = "Spearman",
                           top_annotation = top_annotation,
                           show_row_names = TRUE,
                           show_column_names = TRUE,
                           cluster_rows = TRUE,
                           cluster_columns = TRUE,
                           border = TRUE, 
                           column_title = "Normalized Counts k=2")

# Draw the heatmap k=2
draw(heatmap_spmn_k2)

Past versions of "RUVs HM Spearman-3.png"

####Spearman Heatmap k=3####
heatmap_spmn_k3 <- Heatmap(cor_matrix_spmn_k3,
                           name = "Spearman",
                           top_annotation = top_annotation,
                           show_row_names = TRUE,
                           show_column_names = TRUE,
                           cluster_rows = TRUE,
                           cluster_columns = TRUE,
                           border = TRUE, 
                           column_title = "Normalized Counts k=3")

# Draw the heatmap k=3
draw(heatmap_spmn_k3)

Past versions of "RUVs HM Spearman-4.png"

####Spearman Heatmap k=0 log2cpm####
heatmap_spmn_k0_cpm <- Heatmap(cor_matrix_spmn_k0_cpm,
                           name = "Spearman",
                           top_annotation = top_annotation,
                           show_row_names = TRUE,
                           show_column_names = TRUE,
                           cluster_rows = TRUE,
                           cluster_columns = TRUE,
                           border = TRUE, 
                           column_title = "log2cpm of Filtered Counts no RUVs")

# Draw the heatmap k=0 log2cpm
draw(heatmap_spmn_k0_cpm)

Past versions of "RUVs HM Spearman-5.png"

####Spearman Heatmap k=1 log2cpm####
heatmap_spmn_k1_cpm <- Heatmap(cor_matrix_spmn_k1_cpm,
                           name = "Spearman",
                           top_annotation = top_annotation,
                           show_row_names = TRUE,
                           show_column_names = TRUE,
                           cluster_rows = TRUE,
                           cluster_columns = TRUE,
                           border = TRUE, 
                           column_title = "log2cpm of Normalized Counts k=1")

# Draw the heatmap k=1 log2cpm
draw(heatmap_spmn_k1_cpm)

Past versions of "RUVs HM Spearman-6.png"

####Spearman Heatmap k=2 log2cpm####
heatmap_spmn_k2_cpm <- Heatmap(cor_matrix_spmn_k2_cpm,
                           name = "Spearman",
                           top_annotation = top_annotation,
                           show_row_names = TRUE,
                           show_column_names = TRUE,
                           cluster_rows = TRUE,
                           cluster_columns = TRUE,
                           border = TRUE, 
                           column_title = "log2cpm of Normalized Counts k=2")

# Draw the heatmap k=2 log2cpm
draw(heatmap_spmn_k2_cpm)

Past versions of "RUVs HM Spearman-7.png"

####Spearman Heatmap k=3 log2cpm####
heatmap_spmn_k3_cpm <- Heatmap(cor_matrix_spmn_k3_cpm,
                           name = "Spearman",
                           top_annotation = top_annotation,
                           show_row_names = TRUE,
                           show_column_names = TRUE,
                           cluster_rows = TRUE,
                           cluster_columns = TRUE,
                           border = TRUE, 
                           column_title = "log2cpm of Normalized Counts k=3")

# Draw the heatmap k=3 log2cpm
draw(heatmap_spmn_k3_cpm)

Past versions of "RUVs HM Spearman-8.png"

#Perform Differential Expression Analysis with k = 1 RUVs Correction

#same DGEList object as before
dge <- readRDS("data/new/dge_matrix.RDS")

#check normalization factors from TMM normalization of LIBRARIES
dge$samples

                     group lib.size norm.factors
84-1_DOX_24        DOX_24T 23393931    0.9745263
84-1_DMSO_24      DMSO_24T 22853195    0.9565797
84-1_DOX_24+24     DOX_24R 23846995    1.1659432
84-1_DMSO_24+24   DMSO_24R 21299355    0.9649641
84-1_DOX_24+144   DOX_144R 18222568    0.9913625
84-1_DMSO_24+144 DMSO_144R 28115884    0.9653464
87-1_DOX_24        DOX_24T 19935097    1.0526605
87-1_DMSO_24      DMSO_24T 21302879    0.9773889
87-1_DOX_24+24     DOX_24R 25636959    1.0751043
87-1_DMSO_24+24   DMSO_24R 26319662    0.9940323
87-1_DOX_24+144   DOX_144R 23463426    0.9003102
87-1_DMSO_24+144 DMSO_144R 25840938    0.9888449
78-1_DOX_24        DOX_24T 23085807    0.7676077
78-1_DMSO_24      DMSO_24T 25610495    1.0077383
78-1_DOX_24+24     DOX_24R 18083930    1.1682704
78-1_DMSO_24+24   DMSO_24R 24331177    0.9906872
78-1_DOX_24+144   DOX_144R 19754391    0.9941834
78-1_DMSO_24+144 DMSO_144R 22641509    1.0010734
75-1_DOX_24        DOX_24T 20583626    1.0676786
75-1_DMSO_24      DMSO_24T 28166198    1.0031906
75-1_DOX_24+24     DOX_24R 25831427    1.1530208
75-1_DMSO_24+24   DMSO_24R 26081158    1.0058953
75-1_DOX_24+144   DOX_144R 24659898    0.9261599
75-1_DMSO_24+144 DMSO_144R 25412931    0.9703454
17-3_DOX_24        DOX_24T 22518848    0.9766893
17-3_DMSO_24      DMSO_24T 24589534    0.9612345
17-3_DOX_24+24     DOX_24R 24797547    1.1703079
17-3_DMSO_24+24   DMSO_24R 25977536    0.9509690
17-3_DOX_24+144   DOX_144R 27447106    0.9422729
17-3_DMSO_24+144 DMSO_144R 24893583    0.9356377
90-1_DOX_24        DOX_24T 25187428    1.0311957
90-1_DMSO_24      DMSO_24T 25630519    1.0283437
90-1_DOX_24+24     DOX_24R 26138399    1.1183471
90-1_DMSO_24+24   DMSO_24R 24430396    0.9988688
90-1_DOX_24+144   DOX_144R 23323463    0.9496884
90-1_DMSO_24+144 DMSO_144R 25424152    0.9872926

#read in my covariate information RUV_1 dataframe + remove 6R
RUV_df_rm1 <- readRDS("data/new/RUV_df_rm1.RDS")

#filter out 6R for DE analysis
RUV_1_df_filt <- RUV_df_rm1[!grepl("6R$", rownames(RUV_df_rm1)), , drop = FALSE]

Metadata_2$W_1 <- RUV_1_df_filt$W_1

#now ensure that RUV_1 has the right number of cols after removing rep
length(RUV_1_df_filt$W_1)

[1] 36

#36 

#now make this into a list
RUV_1_DE <-  RUV_1_df_filt$W_1

#saveRDS(RUV_1_DE, "data/new/RUV_1_DE.RDS")

#create my design matrix for DE + RUVs covariate k=1
design1 <- model.matrix(~0 + RUV_1_DE + Metadata_2$Condition)
colnames(design1) <- gsub("Metadata_2\\$Condition", "", colnames(design1))

#take care that the matrix automatically sorts cols alphabetically
##currently DMSO144R, DMSO24R, DMSO24T, DOX144R, DOX24R, DOX24T

#run duplicate correlation for individual effect
corfit1 <- duplicateCorrelation(object = dge$counts, design = design1, block = Metadata_2$Ind)

#voom transformation and plot
v1 <- voom(dge, design1, block = Metadata_2$Ind, correlation = corfit1$consensus.correlation, plot = TRUE)

Past versions of "DE RUVs k1-1.png"

#fit my linear model
fit1 <- lmFit(v1, design1, block = Metadata_2$Ind, correlation = corfit1$consensus.correlation)

#make my contrast matrix to compare across tx and veh
contrast_matrix_RUV <- makeContrasts(
  V.D24T = DOX_24T - DMSO_24T,
  V.D24R = DOX_24R - DMSO_24R, 
  V.D144R = DOX_144R - DMSO_144R,
  RUV_1_24T = RUV_1_DE - DMSO_24T,
  RUV_1_24R = RUV_1_DE - DMSO_24R,
  RUV_1_144R = RUV_1_DE - DMSO_144R,  
  levels = design1
)

#apply these contrasts to compare DOX to DMSO VEH
fit_RUV <- contrasts.fit(fit1, contrast_matrix_RUV)
fit_RUV <- eBayes(fit_RUV)

#plot the mean-variance trend
plotSA(fit_RUV, main = "Mean-Variance Trend, RUVs k=1")

Past versions of "DE RUVs k1-2.png"

#look at the summary of your results
##this tells you the number of DEGs in each condition
results_summary1 <- decideTests(fit_RUV, adjust.method = "BH", p.value = 0.05)

summary(results_summary1)

       V.D24T V.D24R V.D144R RUV_1_24T RUV_1_24R RUV_1_144R
Down     4866   3727     466     13212     13224      13168
NotSig   4899   7107   13659      1014      1002       1058
Up       4554   3485     194        93        93         93

#        V.D24T V.D24R V.D144R RUV_1_24T RUV_1_24R RUV_1_144R
# Down     4866   3727     466     13212     13224      13168
# NotSig   4899   7107   13659      1014      1002       1058
# Up       4554   3485     194        93        93         93

#compare this to my previous DEGs I found
summary(results_summary)

       V.D24T V.D24R V.D144R
Down     4723   3593     359
NotSig   5076   7151   13810
Up       4520   3575     150

#        V.D24T V.D24R V.D144R
# Down     4723   3593     359
# NotSig   5076   7151   13810
# Up       4520   3575     150


#overall, there are more DEGs found after RUVs k=1 correction 
#this was an expected result as it increases tx effect w/ correction

##Generate Toptables RUVs

# Generate Top Table for Specific Comparisons

Toptable_V.D24T_k1 <- topTable(fit = fit_RUV, coef = "V.D24T", number = nrow(x), adjust.method = "BH", p.value = 1, sort.by = "none")
#write.csv(Toptable_V.D24T_k1, "data/new/DEGs/Toptable_V.D24T_k1.csv")

Toptable_V.D24R_k1 <- topTable(fit = fit_RUV, coef = "V.D24R", number = nrow(x), adjust.method = "BH", p.value = 1, sort.by = "none")
#write.csv(Toptable_V.D24R_k1, "data/new/DEGs/Toptable_V.D24R_k1.csv")

Toptable_V.D144R_k1 <- topTable(fit = fit_RUV, coef = "V.D144R", number = nrow(x), adjust.method = "BH", p.value = 1, sort.by = "none")
#write.csv(Toptable_V.D144R_k1, "data/new/DEGs/Toptable_V.D144R_k1.csv")

# save all of these toptables as R objects
# saveRDS(list(
#   V.D24T_k1 = Toptable_V.D24T_k1,
#   V.D24R_k1 = Toptable_V.D24R_k1,
#   V.D144R_k1 = Toptable_V.D144R_k1
# ), file = "data/new/Toptable_list_RUVk1.RDS")

Toptable_list_RUVk1 <- readRDS("data/new/Toptable_list_RUVk1.RDS")

#################################################################
#this section is commented out since it only needs to run once
#kept the code for posterity

# #I want to add the hgnc symbols to my toptables as well
# ####D24T####
# #load in data
# sample_toptab_24T <- read_csv("data/new/DEGs/Toptable_V.D24T_k1.csv", show_col_types = FALSE)
# # ----------------- Ensure Entrez_ID is Present and in Character Format -----------------
# # Check column names
# print(colnames(sample_toptab_24T))
# # Rename if needed (adjust if the column name is not exactly 'Entrez_ID')
# colnames(sample_toptab_24T)[1] <- "Entrez_ID"
# # Convert Entrez_ID to character
# sample_toptab_24T <- sample_toptab_24T %>%
#   mutate(Entrez_ID = as.character(Entrez_ID))
# # ----------------- Map Entrez_ID to Gene Symbol -----------------
# gene_symbols1 <- AnnotationDbi::select(
#   org.Hs.eg.db,
#   keys = sample_toptab_24T$Entrez_ID,
#   columns = c("SYMBOL"),
#   keytype = "ENTREZID"
# )
# # ----------------- Join Back to Main Data -----------------
# Toptable_RUV_24T <- left_join(sample_toptab_24T, gene_symbols1, by = c("Entrez_ID" = "ENTREZID"))
# Toptable_RUV_24T %>% column_to_rownames(., var = "Entrez_ID")
# # ----------------- Save Annotated Output -----------------
# #write_csv(Toptable_RUV_24T, "data/new/DEGs/Toptable_RUV_24T.csv")
# 
# #now do this again for my other two toptables
# 
# ####24R####
# #load in data
# sample_toptab_24R <- read_csv("data/new/DEGs/Toptable_V.D24R_k1.csv", show_col_types = FALSE)
# # ----------------- Ensure Entrez_ID is Present and in Character Format -----------------
# # Check column names
# print(colnames(sample_toptab_24R))
# # Rename if needed (adjust if the column name is not exactly 'Entrez_ID')
# colnames(sample_toptab_24R)[1] <- "Entrez_ID"
# # Convert Entrez_ID to character
# sample_toptab_24R <- sample_toptab_24R %>%
#   mutate(Entrez_ID = as.character(Entrez_ID))
# # ----------------- Map Entrez_ID to Gene Symbol -----------------
# gene_symbols2 <- AnnotationDbi::select(
#   org.Hs.eg.db,
#   keys = sample_toptab_24R$Entrez_ID,
#   columns = c("SYMBOL"),
#   keytype = "ENTREZID"
# )
# # ----------------- Join Back to Main Data -----------------
# Toptable_RUV_24R <- left_join(sample_toptab_24R, gene_symbols2, by = c("Entrez_ID" = "ENTREZID"))
# Toptable_RUV_24R %>% column_to_rownames(., var = "Entrez_ID")
# # ----------------- Save Annotated Output -----------------
# #I'll make these symbols into rownames later for my volcano plots
# #write_csv(Toptable_RUV_24R, "data/new/DEGs/Toptable_RUV_24R.csv")
# 
# 
# ####D144R####
# #load in data
# sample_toptab_144R <- read_csv("data/new/DEGs/Toptable_V.D144R_k1.csv", show_col_types = FALSE)
# # ----------------- Ensure Entrez_ID is Present and in Character Format -----------------
# # Check column names
# print(colnames(sample_toptab_144R))
# # Rename if needed (adjust if the column name is not exactly 'Entrez_ID')
# colnames(sample_toptab_144R)[1] <- "Entrez_ID"
# # Convert Entrez_ID to character
# sample_toptab_144R <- sample_toptab_144R %>%
#   mutate(Entrez_ID = as.character(Entrez_ID))
# # ----------------- Map Entrez_ID to Gene Symbol -----------------
# gene_symbols3 <- AnnotationDbi::select(
#   org.Hs.eg.db,
#   keys = sample_toptab_144R$Entrez_ID,
#   columns = c("SYMBOL"),
#   keytype = "ENTREZID"
# )
# # ----------------- Join Back to Main Data -----------------
# Toptable_RUV_144R <- left_join(sample_toptab_144R, gene_symbols3, by = c("Entrez_ID" = "ENTREZID"))
# Toptable_RUV_144R %>% column_to_rownames(., var = "Entrez_ID")
# # ----------------- Save Annotated Output -----------------
# #I'll make these symbols into rownames later for my volcano plots
# #write_csv(Toptable_RUV_144R, "data/new/DEGs/Toptable_RUV_144R.csv")

# write.csv(Toptable_RUV_24T, "data/new/DEGs/Toptable_RUV_24T_final.csv")
# write.csv(Toptable_RUV_24R, "data/new/DEGs/Toptable_RUV_24R_final.csv")
# write.csv(Toptable_RUV_144R, "data/new/DEGs/Toptable_RUV_144R_final.csv")

Toptable_RUV_24T <- read.csv("data/new/DEGs/Toptable_RUV_24T_final.csv")
Toptable_RUV_24R <- read.csv("data/new/DEGs/Toptable_RUV_24R_final.csv")
Toptable_RUV_144R <- read.csv("data/new/DEGs/Toptable_RUV_144R_final.csv")



# save all of these toptables as R objects
# saveRDS(list(
#   V.D24T_RUV = Toptable_RUV_24T,
#   V.D24R_RUV = Toptable_RUV_24R,
#   V.D144R_RUV = Toptable_RUV_144R
# ), file = "data/new/Toptable_list_RUVk1_Symbols.RDS")

Toptable_list_RUVk1_symbols <- readRDS("data/new/Toptable_list_RUVk1_Symbols.RDS")

##Volcano Plots of Original DEGs

#make a function to generate volcano plots + add gene numbers
generate_volcano_plot <- function(toptable, title) {
  
  #make significance labels
  toptable$Significance <- "Not Significant"
  toptable$Significance[toptable$logFC > 0 & toptable$adj.P.Val < 0.05] <- "Upregulated"
  toptable$Significance[toptable$logFC < 0 & toptable$adj.P.Val < 0.05] <- "Downregulated"
  
  #add number of genes for each significance label
  upgenes <- toptable %>% filter(Significance == "Upregulated") %>% nrow()
  nsgenes <- toptable %>% filter(Significance == "Not Significant") %>% nrow()
  downgenes <- toptable %>% filter(Significance == "Downregulated") %>% nrow()

  #make legend labels for no of genes
  legend_lab <- c(
    str_c("Upregulated: ", upgenes),
    str_c("Not Significant: ", nsgenes),
    str_c("Downregulated: ", downgenes)
  )
  
  #specify the colors for the legend
  legend_col <- c(
    str_c("Upregulated: " = "blue"),
    str_c("Not Significant: " = "gray"),
    str_c("Downregulated: " = "red")
  )

  #generate volcano plot w/ legend
  ggplot(toptable, aes(x = logFC, 
                       y = -log10(P.Value), 
                       color = Significance)) +
    geom_point(alpha = 0.4, size = 2) + 
    scale_color_manual(values = c("Upregulated" = "blue",
                                  "Not Significant" = "gray",
                                  "Downregulated" = "red"), 
                       labels = legend_lab) +
    xlim(-10, 10) +
    labs(title = title, 
         x = expression(x = "log"[2]*"FC"), 
         y = expression(y = "-log"[10]*"P-value")) +
    theme_bw()+
    guides(color = guide_legend(override.aes = list(color = legend_col)))+
    theme(legend.position = "right", 
          plot.title = element_text(size = rel(1.5), hjust = 0.5),
          axis.title = element_text(size = rel(1.25)))
}

#generate volcano plots across each comparison
volcano_plots <- list(
  "V.D24T" = generate_volcano_plot(Toptable_V.D24T, "Volcano Plot DOX 24hr (adj P-val<0.05)"),
  "V.D24R" = generate_volcano_plot(Toptable_V.D24R, "Volcano Plot DOX 24hr Recovery (adj P-val<0.05)"),
  "V.D144R" = generate_volcano_plot(Toptable_V.D144R, "Volcano Plot DOX 144hr Recovery (adj P-val<0.05)")
)

# Display each volcano plot
for (plot_name in names(volcano_plots)) {
  print(volcano_plots[[plot_name]])
}

Past versions of "Volcano Plots of DEGs Original-1.png"

Past versions of "Volcano Plots of DEGs Original-2.png"

Past versions of "Volcano Plots of DEGs Original-3.png"

##Top 5 DEGs after RUVs k=1 Correction

#use your three toptables so I can pull out top 5 genes from each based on adj. p val

top5_D24T_1 <- Toptable_RUV_24T[order(Toptable_RUV_24T$adj.P.Val), ][1:5,] 
top5_D24R_1 <- Toptable_RUV_24R[order(Toptable_RUV_24R$adj.P.Val), ][1:5,] 
top5_D144R_1 <- Toptable_RUV_144R[order(Toptable_RUV_144R$adj.P.Val), ][1:5,]

#now that I've pulled the top 5 DEGs from each, make a list to pull them from my log2cpm data

boxplot1 <- read.csv("data/new/filcpm_final_matrix.csv") %>% 
  as.data.frame()

#Define gene list
#these are the top 5 genes pulled from my toptables
top5_D24T_geneslist_1 <- c(top5_D24T_1$Entrez_ID)
top5_D24R_geneslist_1 <- c(top5_D24R_1$Entrez_ID)
top5_D144R_geneslist_1 <- c(top5_D144R_1$Entrez_ID)
#Add more gene symbols as needed or add more categories

#now pull these from my log2cpm matrix
top5_D24T_genes_1 <- boxplot1[boxplot1$Entrez_ID %in% top5_D24T_geneslist_1,]
dim(top5_D24T_genes_1)

[1]  5 44

#5 genes in 44 cols
top5_D24R_genes_1 <- boxplot1[boxplot1$Entrez_ID %in% top5_D24R_geneslist_1,]
dim(top5_D24R_genes_1)

[1]  5 44

#5 genes in 44 cols
top5_D144R_genes_1 <- boxplot1[boxplot1$Entrez_ID %in% top5_D144R_geneslist_1,]
dim(top5_D144R_genes_1)

[1]  5 44

#5 genes in 44 cols

#Now put in the function I want to use to generate boxplots of genes
#####D24T#####
process_top5_D24T_1 <- function(gene) {
  gene_data <- top5_D24T_genes_1 %>% filter(Entrez_ID == gene)
  long_data <- gene_data %>%
    pivot_longer(cols = -c(Entrez_ID, SYMBOL), names_to = "Sample", values_to = "log2cpm") %>%
    mutate(
      Drug = case_when(
        grepl("DOX", Sample) ~ "DOX",
        grepl("DMSO", Sample) ~ "DMSO",
        TRUE ~ NA_character_
      ),
      Timepoint = case_when(
        grepl("_24T_", Sample) ~ "24T",
        grepl("_24R_", Sample) ~ "24R",
        grepl("_144R_", Sample) ~ "144R",
        TRUE ~ NA_character_
      ),
      Indv = case_when(
        grepl("Ind1$", Sample) ~ "1",
        grepl("Ind2$", Sample) ~ "2",
        grepl("Ind3$", Sample) ~ "3",
        grepl("Ind4$", Sample) ~ "4",
        grepl("Ind5$", Sample) ~ "5",
        grepl("Ind6$", Sample) ~ "6",
        grepl("Ind6R$", Sample) ~ "6R",
        TRUE ~ NA_character_
      ),
      Condition = paste(Drug, Timepoint, sep = "_")
    )
  long_data$Condition <- factor(
    long_data$Condition,
    levels = c(
      "DOX_24T", 
      "DMSO_24T", 
      "DOX_24R", 
      "DMSO_24R", 
      "DOX_144R", 
      "DMSO_144R"
    )
  )
  return(long_data)
}


#Generate Boxplots from the above function using our gene list above
for (gene in top5_D24T_geneslist_1) {
  gene_data <- process_top5_D24T_1(gene)
  p <- ggplot(gene_data, aes(x = Condition, y = log2cpm, fill = Drug)) +
    geom_boxplot(outlier.shape = NA) +
    geom_point(aes(color = Indv), size = 2, alpha = 0.5, position = position_jitter(width = -1, height = 0)) +
    scale_fill_manual(values = c("DOX" = "#499FBD", "DMSO" = "#BBBBBC")) +
    ggtitle(paste("Log2cpm", gene, "top5DEGs D24T RUVs")) +
    labs(x = "Treatment", y = "log2cpm") +
    theme_bw() +
    theme(
      plot.title = element_text(size = rel(2), hjust = 0.5),
      axis.title = element_text(size = 15, color = "black"),
      axis.text.x = element_text(size = 10, color = "black", angle = 90, hjust = 1)
    )
  print(p)
}

Past versions of "RUV Top5 DEGs-1.png"

Past versions of "RUV Top5 DEGs-2.png"

Past versions of "RUV Top5 DEGs-3.png"

Past versions of "RUV Top5 DEGs-4.png"

Past versions of "RUV Top5 DEGs-5.png"

#####D24R#####
process_top5_D24R_1 <- function(gene) {
  gene_data <- top5_D24R_genes_1 %>% filter(Entrez_ID == gene)
  long_data <- gene_data %>%
    pivot_longer(cols = -c(Entrez_ID, SYMBOL), names_to = "Sample", values_to = "log2cpm") %>%
    mutate(
      Drug = case_when(
        grepl("DOX", Sample) ~ "DOX",
        grepl("DMSO", Sample) ~ "DMSO",
        TRUE ~ NA_character_
      ),
      Timepoint = case_when(
        grepl("_24T_", Sample) ~ "24T",
        grepl("_24R_", Sample) ~ "24R",
        grepl("_144R_", Sample) ~ "144R",
        TRUE ~ NA_character_
      ),
      Indv = case_when(
        grepl("Ind1$", Sample) ~ "1",
        grepl("Ind2$", Sample) ~ "2",
        grepl("Ind3$", Sample) ~ "3",
        grepl("Ind4$", Sample) ~ "4",
        grepl("Ind5$", Sample) ~ "5",
        grepl("Ind6$", Sample) ~ "6",
        grepl("Ind6R$", Sample) ~ "6R",
        TRUE ~ NA_character_
      ),
      Condition = paste(Drug, Timepoint, sep = "_")
    )
  long_data$Condition <- factor(
    long_data$Condition,
    levels = c(
      "DOX_24T", 
      "DMSO_24T", 
      "DOX_24R", 
      "DMSO_24R", 
      "DOX_144R", 
      "DMSO_144R"
    )
  )
  return(long_data)
}


#Generate Boxplots from the above function using our gene list above
for (gene in top5_D24R_geneslist_1) {
  gene_data <- process_top5_D24R_1(gene)
  p <- ggplot(gene_data, aes(x = Condition, y = log2cpm, fill = Drug)) +
    geom_boxplot(outlier.shape = NA) +
    geom_point(aes(color = Indv), size = 2, alpha = 0.5, position = position_jitter(width = -1, height = 0)) +
    scale_fill_manual(values = c("DOX" = "#499FBD", "DMSO" = "#BBBBBC")) +
    ggtitle(paste("Log2cpm", gene, "top5DEGs D24R RUVs")) +
    labs(x = "Treatment", y = "log2cpm") +
    theme_bw() +
    theme(
      plot.title = element_text(size = rel(2), hjust = 0.5),
      axis.title = element_text(size = 15, color = "black"),
      axis.text.x = element_text(size = 10, color = "black", angle = 90, hjust = 1)
    )
  print(p)
}

Past versions of "RUV Top5 DEGs-6.png"

Past versions of "RUV Top5 DEGs-7.png"

Past versions of "RUV Top5 DEGs-8.png"

Past versions of "RUV Top5 DEGs-9.png"

Past versions of "RUV Top5 DEGs-10.png"

#####D144R#####
process_top5_D144R_1 <- function(gene) {
  gene_data <- top5_D144R_genes_1 %>% filter(Entrez_ID == gene)
  long_data <- gene_data %>%
    pivot_longer(cols = -c(Entrez_ID, SYMBOL), names_to = "Sample", values_to = "log2cpm") %>%
    mutate(
      Drug = case_when(
        grepl("DOX", Sample) ~ "DOX",
        grepl("DMSO", Sample) ~ "DMSO",
        TRUE ~ NA_character_
      ),
      Timepoint = case_when(
        grepl("_24T_", Sample) ~ "24T",
        grepl("_24R_", Sample) ~ "24R",
        grepl("_144R_", Sample) ~ "144R",
        TRUE ~ NA_character_
      ),
      Indv = case_when(
        grepl("Ind1$", Sample) ~ "1",
        grepl("Ind2$", Sample) ~ "2",
        grepl("Ind3$", Sample) ~ "3",
        grepl("Ind4$", Sample) ~ "4",
        grepl("Ind5$", Sample) ~ "5",
        grepl("Ind6$", Sample) ~ "6",
        grepl("Ind6R$", Sample) ~ "6R",
        TRUE ~ NA_character_
      ),
      Condition = paste(Drug, Timepoint, sep = "_")
    )
  long_data$Condition <- factor(
    long_data$Condition,
    levels = c(
      "DOX_24T", 
      "DMSO_24T", 
      "DOX_24R", 
      "DMSO_24R", 
      "DOX_144R", 
      "DMSO_144R"
    )
  )
  return(long_data)
}


#Generate Boxplots from the above function using our gene list above
for (gene in top5_D144R_geneslist_1) {
  gene_data <- process_top5_D144R_1(gene)
  p <- ggplot(gene_data, aes(x = Condition, y = log2cpm, fill = Drug)) +
    geom_boxplot(outlier.shape = NA) +
    geom_point(aes(color = Indv), size = 2, alpha = 0.5, position = position_jitter(width = -1, height = 0)) +
    scale_fill_manual(values = c("DOX" = "#499FBD", "DMSO" = "#BBBBBC")) +
    ggtitle(paste("Log2cpm", gene, "top5DEGs D144R RUVs")) +
    labs(x = "Treatment", y = "log2cpm") +
    theme_bw() +
    theme(
      plot.title = element_text(size = rel(2), hjust = 0.5),
      axis.title = element_text(size = 15, color = "black"),
      axis.text.x = element_text(size = 10, color = "black", angle = 90, hjust = 1)
    )
  print(p)
}

Past versions of "RUV Top5 DEGs-11.png"

Past versions of "RUV Top5 DEGs-12.png"

Past versions of "RUV Top5 DEGs-13.png"

Past versions of "RUV Top5 DEGs-14.png"

Past versions of "RUV Top5 DEGs-15.png"

##LogFC DEGs RUVs k=1 All Conditions

# Toptable_RUV_24T
# Toptable_RUV_24R
# Toptable_RUV_144R

# Load DEGs Data
DOX_24T_1 <- read.csv("data/new/DEGs/Toptable_RUV_24T_final.csv")
DOX_24R_1 <- read.csv("data/new/DEGs/Toptable_RUV_24R_final.csv")
DOX_144R_1 <- read.csv("data/new/DEGs/Toptable_RUV_144R_final.csv")

#make a list of all of the genes in this set so I can plot the logFC in other sets
D24T_DEGs_1 <- DOX_24T_1$Entrez_ID[DOX_24T_1$adj.P.Val < 0.05]
length(D24T_DEGs_1)

[1] 9420

#9243 genes

D24R_DEGs_1 <- DOX_24R_1$Entrez_ID[DOX_24R_1$adj.P.Val < 0.05]
length(D24R_DEGs_1)

[1] 7212

#7168 genes

D144R_DEGs_1 <- DOX_144R_1$Entrez_ID[DOX_144R_1$adj.P.Val < 0.05]
length(D144R_DEGs_1)

[1] 660

#509 genes

#plot those full gene sets in logFC
D24T_DEGs_1 <- DOX_24T$Entrez_ID[DOX_24T$adj.P.Val < 0.05]
length(D24T_DEGs)

[1] 9243

D24R_DEGs_1 <- DOX_24R$Entrez_ID[DOX_24R$adj.P.Val < 0.05]
length(D24R_DEGs)

[1] 7168

D144R_DEGs_1 <- DOX_144R$Entrez_ID[DOX_144R$adj.P.Val < 0.05]
length(D144R_DEGs)

[1] 509

#now that I have the full list of genes, I want to plot the logFC across conditions


#Combine the toptables I have from pairwise analysis into a single dataframe
d24_toptable_dxr_1 <- Toptable_RUV_24T %>% 
  mutate(Time = "24")

d24r_toptable_dxr_1 <- Toptable_RUV_24R %>% 
  mutate(Time = "24R")

d144r_toptable_dxr_1 <- Toptable_RUV_144R %>% 
  mutate(Time = "144R")

combined_toptables_dxr_RUV <- bind_rows(
  d24_toptable_dxr_1,
  d24r_toptable_dxr_1,
  d144r_toptable_dxr_1)

#Filter the data based on each motif
filt_toptable_dxr_RUV <- combined_toptables_dxr_RUV %>% 
  dplyr::filter(Entrez_ID  %in% D24T_DEGs_1) %>% 
  mutate(absFC = abs(logFC)) %>% 
  mutate(time = factor(Time, levels = c("24", "24R", "144R"), labels = c("24hr", "24hr Recovery", "144hr Recovery"))) %>% 
  ggplot(., aes(x = time, y = logFC))+
  geom_boxplot(aes(fill = time))+
  scale_fill_manual(values = time_col)+
  guides(fill = guide_legend(title = "Time"))+
  theme_bw()+
  xlab(" ")+
  ylab("logFC")+
  theme_bw()+
  ggtitle("LogFC D24T DEGs Across Conditions RUVs")+
  theme(plot.title = element_text(size = rel(1.5), hjust = 0.5),
        axis.title = element_text(size = 15, colour = "black"),
        strip.background = element_rect(fill = "#CAD899"),
        axis.text.x = element_text(size = 8, colour = "white", angle = 15),
        strip.text.x = element_text(size = 12, colour = "black", face = "bold"))

#now do the same with the other conditions

filt_toptable_dxr_24r_RUV <- combined_toptables_dxr_RUV %>% 
  dplyr::filter(Entrez_ID  %in% D24R_DEGs_1) %>% 
  mutate(absFC = abs(logFC)) %>% 
  mutate(time = factor(Time, levels = c("24", "24R", "144R"), labels = c("24hr", "24hr Recovery", "144hr Recovery"))) %>% 
  ggplot(., aes(x = time, y = logFC))+
  geom_boxplot(aes(fill = time))+
  scale_fill_manual(values = time_col)+
  guides(fill = guide_legend(title = "Time"))+
  theme_bw()+
  xlab(" ")+
  ylab("logFC")+
  theme_bw()+
  ggtitle("LogFC D24R DEGs Across Conditions RUVs")+
  theme(plot.title = element_text(size = rel(1.5), hjust = 0.5),
        axis.title = element_text(size = 15, colour = "black"),
        strip.background = element_rect(fill = "#CAD899"),
        axis.text.x = element_text(size = 8, colour = "white", angle = 15),
        strip.text.x = element_text(size = 12, colour = "black", face = "bold"))

#D144R

filt_toptable_dxr_144r_RUV <- combined_toptables_dxr_RUV %>% 
  dplyr::filter(Entrez_ID  %in% D144R_DEGs_1) %>% 
  mutate(absFC = abs(logFC)) %>% 
  mutate(time = factor(Time, levels = c("24", "24R", "144R"), labels = c("24hr", "24hr Recovery", "144hr Recovery"))) %>% 
  ggplot(., aes(x = time, y = logFC))+
  geom_boxplot(aes(fill = time))+
  scale_fill_manual(values = time_col)+
  guides(fill = guide_legend(title = "Time"))+
  theme_bw()+
  xlab(" ")+
  ylab("logFC")+
  theme_bw()+
  ggtitle("LogFC D144R DEGs Across Conditions RUVs")+
  theme(plot.title = element_text(size = rel(1.5), hjust = 0.5),
        axis.title = element_text(size = 15, colour = "black"),
        strip.background = element_rect(fill = "#CAD899"),
        axis.text.x = element_text(size = 8, colour = "white", angle = 15),
        strip.text.x = element_text(size = 12, colour = "black", face = "bold"))

#now put the names of these graphs to print them
filt_toptable_dxr_RUV

Past versions of "RUV logFC DEGs all Conditions-1.png"

filt_toptable_dxr_24r_RUV

Past versions of "RUV logFC DEGs all Conditions-2.png"

filt_toptable_dxr_144r_RUV

Past versions of "RUV logFC DEGs all Conditions-3.png"

##Volcano Plots of RUVs Corrected Data k=1

#make a function to generate volcano plots + add gene numbers
#ensure the gene symbols are in a col so I can plot names of top 15
generate_volcano_plot_RUV <- function(toptable, title) {
  
  #make significance labels
  toptable$Significance <- "Not Significant"
  toptable$Significance[toptable$logFC > 0 & toptable$adj.P.Val < 0.05] <- "Upregulated"
  toptable$Significance[toptable$logFC < 0 & toptable$adj.P.Val < 0.05] <- "Downregulated"
  
  #add number of genes for each significance label
  upgenes <- toptable %>% filter(Significance == "Upregulated") %>% nrow()
  nsgenes <- toptable %>% filter(Significance == "Not Significant") %>% nrow()
  downgenes <- toptable %>% filter(Significance == "Downregulated") %>% nrow()

  #make legend labels for no of genes
  legend_lab <- c(
    str_c("Upregulated: ", upgenes),
    str_c("Not Significant: ", nsgenes),
    str_c("Downregulated: ", downgenes)
  )
  
  #specify the colors for the legend
  legend_col <- c(
    str_c("Upregulated: " = "blue"),
    str_c("Not Significant: " = "gray"),
    str_c("Downregulated: " = "red")
  )
  
  #add the top 15 genes by adj. p value
  top_genes <- toptable %>%
    filter(!is.na(SYMBOL)) %>%
    arrange(adj.P.Val) %>%
    slice_head(n = 15)
  
  #generate volcano plot w/ legend
  ggplot(toptable, aes(x = logFC, 
                       y = -log10(P.Value), 
                       color = Significance)) +
    geom_point(alpha = 0.4, size = 2) + 
    scale_color_manual(values = c("Upregulated" = "blue",
                                  "Not Significant" = "gray",
                                  "Downregulated" = "red"), 
                       labels = legend_lab) +
    xlim(-10, 10) +
    labs(title = title, 
         x = expression(x = "log"[2]*"FC"), 
         y = expression(y = "-log"[10]*"P-value")) +
    theme_bw()+
    guides(color = guide_legend(override.aes = list(color = legend_col)))+
    theme(legend.position = "right", 
          plot.title = element_text(size = rel(1.5), hjust = 0.5),
          axis.title = element_text(size = rel(1.25)))
}

#generate volcano plots across each comparison
volcano_plots_RUV <- list(
  "V.D24T_RUV" = generate_volcano_plot(Toptable_RUV_24T, "DOX 24T RUVs k=1 (adj P-val<0.05)"),
  "V.D24R_RUV" = generate_volcano_plot(Toptable_RUV_24R, "DOX 24R RUVs k=1 (adj P-val<0.05)"),
  "V.D144R_RUV" = generate_volcano_plot(Toptable_RUV_144R, "DOX 144R RUVs k=1 (adj P-val<0.05)")
)

# Display each volcano plot
for (plot_name in names(volcano_plots_RUV)) {
  print(volcano_plots_RUV[[plot_name]])
}

Past versions of "Volcano Plots RUVs Corrected k1-1.png"

Past versions of "Volcano Plots RUVs Corrected k1-2.png"

Past versions of "Volcano Plots RUVs Corrected k1-3.png"

##Plot DNA Damage Response Genes (DDR) with Original Data

#DDR Gene Expression Heatmap — DOX Over Recovery Time (68 genes, with categories)

# Load libraries
# library(circlize)
# library(grid)
# library(reshape2)

# Load DEG files
load_deg <- function(path) read.csv(path)

DOX_24T <- load_deg("data/new/DEGs/Toptable_V.D24T.csv")
DOX_24R <- load_deg("data/new/DEGs/Toptable_V.D24R.csv")
DOX_144R <- load_deg("data/new/DEGs/Toptable_V.D144R.csv")

# Final Entrez IDs and categories (68 genes)
entrez_category_DDR <- tribble(
  ~ENTREZID, ~Category,
  317, "Apoptosis", 355, "Apoptosis", 581, "Apoptosis", 637, "Apoptosis",
  836, "Apoptosis", 841, "Apoptosis", 842, "Apoptosis", 27113, "Apoptosis",
  5366, "Apoptosis", 54205, "Apoptosis", 55367, "Apoptosis", 8795, "Apoptosis",
  1026, "Cell Cycle / Checkpoint", 1027, "Cell Cycle / Checkpoint", 595, "Cell Cycle / Checkpoint",
  894, "Cell Cycle / Checkpoint", 896, "Cell Cycle / Checkpoint", 898, "Cell Cycle / Checkpoint",
  9133, "Cell Cycle / Checkpoint", 9134, "Cell Cycle / Checkpoint", 891, "Cell Cycle / Checkpoint",
  983, "Cell Cycle / Checkpoint", 1017, "Cell Cycle / Checkpoint", 1019, "Cell Cycle / Checkpoint",
  1020, "Cell Cycle / Checkpoint", 1021, "Cell Cycle / Checkpoint", 993, "Cell Cycle / Checkpoint",
  995, "Cell Cycle / Checkpoint", 1869, "Cell Cycle / Checkpoint", 4609, "Cell Cycle / Checkpoint",
  5925, "Cell Cycle / Checkpoint", 9874, "Cell Cycle / Checkpoint", 11011, "Cell Cycle / Checkpoint",
  1385, "Cell Cycle / Checkpoint",
  472, "Damage Sensors / Signal Transducers", 545, "Damage Sensors / Signal Transducers",
  5591, "Damage Sensors / Signal Transducers", 5810, "Damage Sensors / Signal Transducers",
  5883, "Damage Sensors / Signal Transducers", 5884, "Damage Sensors / Signal Transducers",
  6118, "Damage Sensors / Signal Transducers", 4361, "Damage Sensors / Signal Transducers",
  10111, "Damage Sensors / Signal Transducers", 4683, "Damage Sensors / Signal Transducers",
  84126, "Damage Sensors / Signal Transducers", 3014, "Damage Sensors / Signal Transducers",
  672, "DNA Repair", 2177, "DNA Repair", 5888, "DNA Repair", 5893, "DNA Repair",
  1647, "DNA Repair", 4616, "DNA Repair", 10912, "DNA Repair", 1111, "DNA Repair",
  11200, "DNA Repair", 1643, "DNA Repair", 8243, "DNA Repair", 5981, "DNA Repair",
  7157, "p53 Regulators / Targets", 4193, "p53 Regulators / Targets", 5371, "p53 Regulators / Targets",
  27244, "p53 Regulators / Targets", 50484, "p53 Regulators / Targets",
  5916, "DOX Cardiotoxicity", 7799, "DOX Cardiotoxicity", 4292, "DOX Cardiotoxicity",
  207, "Miscellaneous / Broad", 25, "Miscellaneous / Broad"
)

entrez_ids_DDR <- entrez_category_DDR$ENTREZID

# Extract relevant DEG values
extract_data_DDR <- function(df, name) {
  df %>%
    filter(Entrez_ID %in% entrez_ids_DDR) %>%
    mutate(
      Gene = mapIds(org.Hs.eg.db, as.character(Entrez_ID),
                    column = "SYMBOL", keytype = "ENTREZID", multiVals = "first"),
      Condition = name,
      Signif = ifelse(adj.P.Val < 0.05, "*", "")
    )
}

# DEG list
deg_list <- list("DOX_24T" = DOX_24T, 
                 "DOX_24R" = DOX_24R, 
                 "DOX_144R" = DOX_144R
)

# Combine all DEGs and annotate
all_data_DDR <- bind_rows(mapply(extract_data_DDR, deg_list, names(deg_list), SIMPLIFY = FALSE)) %>%
  left_join(entrez_category_DDR, by = c("Entrez_ID" = "ENTREZID"))

'select()' returned 1:1 mapping between keys and columns
'select()' returned 1:1 mapping between keys and columns
'select()' returned 1:1 mapping between keys and columns

# Create matrices
logFC_matddr <- acast(all_data_DDR, Gene ~ Condition, value.var = "logFC")
signif_matddr <- acast(all_data_DDR, Gene ~ Condition, value.var = "Signif")

# Set desired order
desired_order <- c("DOX_24T",
                   "DOX_24R",
                   "DOX_144R")

logFC_mat_DDR <- logFC_matddr[, desired_order, drop = FALSE]
signif_mat_DDR <- signif_matddr[, desired_order, drop = FALSE]

# Column annotation
meta_DDR <- str_split_fixed(colnames(logFC_mat_DDR), "_", 2)
col_annot <- HeatmapAnnotation(
  Drug = meta_DDR[, 1],
  Time = meta_DDR[, 2],
  col = list(
    Drug = c("DOX" = "#499FBD"),
    Time = c("24T" = "#238B45", 
             "24R" = "#74C476", 
             "144R" = "#C7E9C0")
  ),
  annotation_height = unit(c(1, 1, 1), "cm")
)

# Row annotation
gene_order_df_DDR <- all_data_DDR %>%
  distinct(Gene, Category) %>%
  arrange(factor(Category, levels = sort(unique(entrez_category_DDR$Category))), Gene)

ordered_genes_DDR <- gene_order_df_DDR$Gene
logFC_mat_DDR <- logFC_mat_DDR[ordered_genes_DDR, ]
signif_mat_DDR <- signif_mat_DDR[ordered_genes_DDR, ]

category_colors_DDR <- structure(
  c("darkorange", "steelblue", "darkgreen", "firebrick", "gold", "mediumpurple", "gray60"),
  names = sort(unique(entrez_category_DDR$Category))
)

ha_left_DDR <- rowAnnotation(
  Category = gene_order_df_DDR$Category,
  col = list(Category = category_colors_DDR),
  annotation_name_side = "top"
)

# Final Heatmap
Heatmap(logFC_mat_DDR,
        name = "logFC",
        top_annotation = col_annot,
        left_annotation = ha_left_DDR,
        cluster_columns = FALSE,
        cluster_rows = FALSE,
        show_row_names = TRUE,
        show_column_names = FALSE,
        row_names_gp = gpar(fontsize = 10),
        column_title = "DDR Gene Expression Response (n = 68)\n DOX Recovery",
        column_title_gp = gpar(fontsize = 14, fontface = "bold"),
        cell_fun = function(j, i, x, y, width, height, fill) {
          grid.text(signif_mat_DDR[i, j], x, y, gp = gpar(fontsize = 9))
        }
)

Past versions of "DDR Genes logFC-1.png"

##Plot DOX Cardiotox Genes logFC HM Original Data

#plot a list of 29 functionally validated DIC genes
# Load DEG files
load_deg <- function(path) read.csv(path)

DOX_24T <- load_deg("data/new/DEGs/Toptable_V.D24T.csv")
DOX_24R <- load_deg("data/new/DEGs/Toptable_V.D24R.csv")
DOX_144R <- load_deg("data/new/DEGs/Toptable_V.D144R.csv")

#the data I have is in hgnc_symbols, I want to convert this to entrez_id like my df
DIC_genes <- tribble(
  ~SYMBOL, ~Category,
  "CAT", "ROS Generation / Handling",
  "CBR1", "ROS Generation / Handling",
  "CBR3", "ROS Generation / Handling",
  "ERBB2", "ROS Generation / Handling",
  "GPX3", "ROS Generation / Handling",
  "GSTM1", "ROS Generation / Handling",
  "GSTP", "ROS Generation / Handling",
  "HAS3", "ROS Generation / Handling",
  "NOS3", "ROS Generation / Handling",
  "PLCE1", "ROS Generation / Handling",
  "RAC2", "ROS Generation / Handling",
  "SPG7", "ROS Generation / Handling",
  "PRDM2", "DNA Damage", 
  "MLH1", "DNA Damage",
  "RARG", "DNA Damage",
  "HFE", "Iron Uptake & Homeostasis",
  "SLC22A17", "DOX Uptake",
  "SLC28A1", "DOX Uptake",
  "SLC28A3", "DOX Uptake",
  "ABCB4", "DOX Efflux",
  "ABCC2", "DOX Efflux",
  "ABCC5", "DOX Efflux", 
  "ABCC9", "DOX Efflux",
  "ABCC10", "DOX Efflux",
  "CELF4", "Calcium Handling",
  "MYH7", "Calcium Handling",
  "CYP2J2", "Cardiac Electrical Activity",
  "RIN3", "Cardiac Electrical Activity",
  "ZFN521", "Cardiac Electrical Activity")

gene_df_DIC <- tibble(HGNC = DIC_genes)

gene_df_DIC <- gene_df_DIC %>% 
  mutate(Entrez_ID = mapIds(org.Hs.eg.db,
                            keys = DIC_genes$SYMBOL,
                            column = "ENTREZID", 
                            keytype = "SYMBOL", 
                            multiVals = "first")) %>% 
  unnest_wider(HGNC) %>% 
  mutate(Entrez_ID = as.character(Entrez_ID))

'select()' returned 1:1 mapping between keys and columns

#now I've put together a dataframe with the HGNC, Category, and Entrez_ID
#plus I've ensured that Entrez_ID is a character for later joining

entrez_ids_DIC <- gene_df_DIC$Entrez_ID

# Extract relevant DEG values
extract_data_DIC <- function(df, name) {
  df %>%
    filter(Entrez_ID %in% entrez_ids_DIC) %>%
    mutate(
      Gene = mapIds(org.Hs.eg.db, as.character(Entrez_ID),
                    column = "SYMBOL", keytype = "ENTREZID", multiVals = "first"),
      Condition = name,
      Signif = ifelse(adj.P.Val < 0.05, "*", "")
    )
}

# DEG list
deg_list <- list("DOX_24T" = DOX_24T, 
                 "DOX_24R" = DOX_24R, 
                 "DOX_144R" = DOX_144R
)

# Combine all DEGs and annotate
all_data_DIC <- bind_rows(mapply(extract_data_DIC, 
                             deg_list, 
                             names(deg_list), 
                             SIMPLIFY = FALSE)) %>%
  mutate(Entrez_ID = as.character(Entrez_ID)) %>% 
  left_join(gene_df_DIC, 
            by = "Entrez_ID")

'select()' returned 1:1 mapping between keys and columns
'select()' returned 1:1 mapping between keys and columns
'select()' returned 1:1 mapping between keys and columns

#I've also made sure here that Entrez_ID is a character and not an integer

# Create matrices
logFC_matdic <- acast(all_data_DIC, Gene ~ Condition, value.var = "logFC")
signif_matdic <- acast(all_data_DIC, Gene ~ Condition, value.var = "Signif")

# Set desired order
desired_order <- c("DOX_24T",
                   "DOX_24R",
                   "DOX_144R")

logFC_mat_DIC <- logFC_matdic[, desired_order, drop = FALSE]
signif_mat_DIC <- signif_matdic[, desired_order, drop = FALSE]

# Column annotation
meta_DIC <- str_split_fixed(colnames(logFC_mat_DIC), "_", 3)
meta_DIC <- str_split_fixed(colnames(logFC_mat_DIC), "_", 2)
col_annot <- HeatmapAnnotation(
  Drug = meta_DIC[, 1],
  Time = meta_DIC[, 2],
  col = list(
    Drug = c("DOX" = "#499FBD"),
    Time = c("24T" = "#238B45", 
             "24R" = "#74C476", 
             "144R" = "#C7E9C0")
  ),
  annotation_height = unit(c(2, 2, 2), "cm")
)
#for this I can leave off the DMSO as it's already from a pairwise comparison

#make an extra object with my category order as well
category_order_DIC <- c(
  "ROS Generation / Handling",
  "DNA Damage",
  "Calcium Handling",
  "DOX Uptake",
  "Iron Uptake & Homeostasis",
  "DOX Efflux",
  "Cardiac Electrical Activity"
)

# Row annotation
gene_order_df_DIC <- all_data_DIC %>%
  distinct(Gene, Category) %>%
  mutate(Category = factor(Category, levels = category_order_DIC)) %>% 
  arrange(Category, Gene)

ordered_genes_DIC <- gene_order_df_DIC$Gene
logFC_mat_DIC <- logFC_mat_DIC[ordered_genes_DIC, ]
signif_mat_DIC <- signif_mat_DIC[ordered_genes_DIC, ]

#add in your colors for each category
 category_colors_DIC <- structure(
  c("darkorange", "steelblue", "darkgreen", "firebrick", "gold", "mediumpurple", "gray60"),
  names = category_order_DIC
)

ha_left_DIC <- rowAnnotation(
  Category = gene_order_df_DIC$Category,
  col = list(Category = category_colors_DIC),
  annotation_name_side = "top"
)

# Draw heatmap
Heatmap(logFC_mat_DIC,
        name = "logFC",
        top_annotation = col_annot,
        left_annotation = ha_left_DIC,
        cluster_columns = FALSE,
        cluster_rows = FALSE,
        show_row_names = TRUE,
        show_column_names = FALSE,
        cell_fun = function(j, i, x, y, width, height, fill) {
          grid.text(signif_mat_DIC[i, j], x, y, gp = gpar(fontsize = 9))
        },
        column_title = "DIC Genes Expression (n=29)\nDOX Recovery",
        column_title_gp = gpar(fontsize = 14, fontface = "bold")
)

Past versions of "DIC Genes logFC HM noRUVs-1.png"

##Plot p53 Target Genes logFC HM Original Data

# Load libraries
library(tidyverse)
library(ComplexHeatmap)
library(circlize)
library(grid)
library(org.Hs.eg.db)
library(reshape2)

# Load DEG files
load_deg <- function(path) read.csv(path)

DOX_24T <- load_deg("data/new/DEGs/Toptable_V.D24T.csv")
DOX_24R <- load_deg("data/new/DEGs/Toptable_V.D24R.csv")
DOX_144R <- load_deg("data/new/DEGs/Toptable_V.D144R.csv")

#P53 target Entrez_ID
entrez_ids_p53 <- c(1026,50484,4193,9766,9518,7832,1643,1647,1263,57103,51065,8795,51499,64393,581,
                5228,5429,8493,55959,7508,64782,282991,355,53836,4814,10769,9050,27244,9540,94241,
                26154,57763,900,26999,55332,26263,23479,23612,29950,9618,10346,8824,134147,55294,
                22824,4254,6560,467,27113,60492,8444,60401,1969,220965,2232,3976,55191,84284,93129,
                5564,7803,83667,7779,132671,7039,51768,137695,93134,7633,10973,340485,307,27350,
                23245,3732,29965,1363,1435,196513,8507,8061,2517,51278,53354,54858,23228,5366,5912,
                6236,51222,26152,59,1907,50650,91012,780,9249,11072,144455,64787,116151,27165,2876,
                57822,55733,57722,121457,375449,85377,4851,5875,127544,29901,84958,8797,8793,441631,
                220001,54541,5889,5054,25816,25987,5111,98,317,598,604,10904,1294,80315,53944,
                1606,2770,3628,3675,3985,4035,4163,84552,29085,55367,5371,5791,54884,5980,8794,
                1462,50808,220,583,694,1056,9076,10978,54677,1612,55040,114907,2274,127707,4000,
                8079,4646,4747,27445,5143,80055,79156,5360,5364,23654,5565,5613,5625,10076,56963,
                6004,390,255488,6326,6330,23513,7869,283130,204962,83959,6548,6774,9263,10228,
                22954,10475,85363,494514,10142,79714,1006,8446,9648,79828,5507,55240,63874,25841,
                9289,84883,154810,51321,421,8553,655,119032,84280,10950,824,839,57828,857,8812,
                8837,94027,113189,22837,132864,10898,3300,81704,1847,1849,1947,9538,24139,5168,
                147965,115548,9873,23768,2632,2817,3280,3265,23308,3490,51477,182,3856,8844,144811,
                9404,4043,9848,2872,23041,740,343263,4638,26509,4792,22861,57523,55214,80025,164091,
                57060,64065,51090,5453,8496,333926,55671,5900,55544,23179,8601,389,6223,55800,6385,
                4088,6643,122809,257397,285343,7011,54790,374618,55362,51754,7157,9537,22906,7205,
                80705,219699,55245,83719,7748,25946,118738)

# Function to extract relevant data
extract_data_p53 <- function(df, name) {
  df %>%
    filter(Entrez_ID %in% entrez_ids_p53) %>%
    mutate(Gene = mapIds(org.Hs.eg.db, as.character(Entrez_ID),
                         column = "SYMBOL", keytype = "ENTREZID", multiVals = "first"),
           Condition = name,
           Signif = ifelse(adj.P.Val < 0.05, "*", ""))
}

# Collect all data
# DEG list
deg_list <- list("DOX_24T" = DOX_24T, 
                 "DOX_24R" = DOX_24R, 
                 "DOX_144R" = DOX_144R
)

# Combine all DEGs and annotate
all_data_p53 <- bind_rows(mapply(extract_data_p53, deg_list, names(deg_list), SIMPLIFY = FALSE)) 

'select()' returned 1:1 mapping between keys and columns
'select()' returned 1:1 mapping between keys and columns
'select()' returned 1:1 mapping between keys and columns

#unnecessary to join again as there are not categories on this list of entrezids
# %>%
#   left_join(entrez_ids, by = c("Entrez_ID" = "ENTREZID"))

# Create matrices
logFC_mat53 <- acast(all_data_p53, Gene ~ Condition, value.var = "logFC")
signif_mat53 <- acast(all_data_p53, Gene ~ Condition, value.var = "Signif")

# Desired column order
desired_order <- c("DOX_24T",
                   "DOX_24R",
                   "DOX_144R")

logFC_mat_p53 <- logFC_mat53[, desired_order]
signif_mat_p53 <- signif_mat53[, desired_order]

# Column annotation
meta_p53 <- str_split_fixed(colnames(logFC_mat_p53), "_", 3)
meta_p53 <- str_split_fixed(colnames(logFC_mat_p53), "_", 2)
col_annot <- HeatmapAnnotation(
  Drug = meta_p53[, 1],
  Time = meta_p53[, 2],
  col = list(
    Drug = c("DOX" = "#499FBD"),
    Time = c("24T" = "#238B45", 
             "24R" = "#74C476", 
             "144R" = "#C7E9C0")
  ),
  annotation_height = unit(c(1, 1, 1), "cm")
)

# Draw heatmap
Heatmap(logFC_mat_p53,
        name = "logFC",
        top_annotation = col_annot,
        cluster_columns = FALSE,
        cluster_rows = TRUE,
        show_row_names = TRUE,
        show_column_names = FALSE,
        layer_fun = function(j, i, x, y, width, height, fill) {
          grid.text(signif_mat_p53[cbind(i, j)], x, y, gp = gpar(fontsize = 9))
        },
        column_title = "P53 Target Genes Expression\nDOX Recovery",
        column_title_gp = gpar(fontsize = 14, fontface = "bold")
)

Past versions of "p53 Target Genes logFC HM-1.png"

##Create a Venn Diagram of DEGs from Original DEA

#plot a venn diagram with all of your conditions from your toptables

# Load DEGs Data
DOX_24T <- read.csv("data/new/DEGs/Toptable_V.D24T.csv")
DOX_24R <- read.csv("data/new/DEGs/Toptable_V.D24R.csv")
DOX_144R <- read.csv("data/new/DEGs/Toptable_V.D144R.csv")

# Extract Significant DEGs
DEG1 <- DOX_24T$Entrez_ID[DOX_24T$adj.P.Val < 0.05]
DEG2 <- DOX_24R$Entrez_ID[DOX_24R$adj.P.Val < 0.05]
DEG3 <- DOX_144R$Entrez_ID[DOX_144R$adj.P.Val < 0.05]


venntest <- list(DEG1, DEG2, DEG3)
ggVennDiagram(
  venntest,
  category.names = c("DOX_24T", "DOX_24R", "DOX_144R")
) + ggtitle("DXR Specific and Shared DEGs")+
  theme(
    plot.title = element_text(size = 16, face = "bold"),  # Increase title size
    text = element_text(size = 16)  # Increase text size globally
  )

Past versions of "Overlap of DEGs Original-1.png"

#Now that I've made my venn diagram, I want to compare these DEGs
#set 1 : 4362 DOX24T specific genes
#set2 : 4362 + 4550 + 50 + 272 genes shared across DOX24T (all genes)
#how many of these are downregulated and how many are upregulated?

# Extract Significant DEGs
# Create a list of DEGs for each sample

# Example gene sets

DEG1 <- DOX_24T$Entrez_ID[DOX_24T$adj.P.Val < 0.05]
DEG2 <- DOX_24R$Entrez_ID[DOX_24R$adj.P.Val < 0.05]
DEG3 <- DOX_144R$Entrez_ID[DOX_144R$adj.P.Val < 0.05]

#try and use the VennDetail package to extract the genes from each condition
#Set 1 - DOX_24T only genes
#Set 2 - DOX_24T shared genes
plot.new()
venn_test <- venndetail(venntest)
plot(venn_test)

Past versions of "Overlap of DEGs Original-2.png"

detail(venn_test)

=== Here is the detail of Venndiagram ===
Total results:  11612 x 2 
Total sets is: 7 
  Subset Detail
1 Shared 643837
2 Shared   1870
3 Shared  84958
4 Shared   2537
5 Shared 192670
6 Shared   1298
... with 11606 more rows ...

#now that I have the genes subsetted by condition, I can pull out the genes I want to look at

venn_DOX_24T <- getSet(object = venn_test, subset = c("Group 1"))

Warning: `filter_()` was deprecated in dplyr 0.7.0.
ℹ Please use `filter()` instead.
ℹ See vignette('programming') for more help
ℹ The deprecated feature was likely used in the VennDetail package.
  Please report the issue to the authors.
This warning is displayed once every 8 hours.
Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
generated.

dim(venn_DOX_24T)

[1] 4362    2

#4362 genes in DOX_24T only

venn_DOX_shared <- getSet(object = venn_test, subset = c("Shared", "Group 1", "Group 1_Group 2", "Group 1_Group 3"))
dim(venn_DOX_shared)

[1] 9243    2

#total of 9243 genes:
  #4559 DOX_24T only 
  #272 shared all
  #50 DOX_24T + DOX_144R
  #4362 DOX_24T + DOX_24R

venn_DOX144R_shared <- getSet(object = venn_test, subset = c("Shared", "Group 3", "Group 2_Group 3", "Group 1_Group 3"))
dim(venn_DOX144R_shared)

[1] 509   2

#total of 509 genes:
 #272 shared all
 #32 DOX144R Specific
 #50 DOX24T + DOX144R
 #155 DOX24R + DOX144R

#now I can look at these sets to see which ones are up and down regulated in each
#after that, run GO analysis
venn_shared_DEGs <- venn_DOX_shared %>% 
  as.data.frame() %>% 
  column_to_rownames(., var = "Detail")

venn_sharedD144R_DEGs <- venn_DOX144R_shared %>% 
  as.data.frame() %>% 
  column_to_rownames(., var = "Detail")

venn_DOX24T_DEGs <- venn_DOX_24T %>% 
  as.data.frame() %>% 
  column_to_rownames(., var = "Detail")


#pull these genes out of my DEG matrix with logFC
DOX_24T_DEGs <- Toptable_list$V.D24T[row.names(venn_DOX24T_DEGs),]
DOX_24T_shared_DEGs <- Toptable_list$V.D24T[row.names(venn_shared_DEGs),]
DOX_144R_shared_DEGs <- Toptable_list$V.D144R[row.names(venn_sharedD144R_DEGs),]
#I want to go ahead and do this for every condition so I can see the genes in there
#I also want to filter these by their logFC being up or down for GO/KEGG

DOX24T_DEGs_GO <- DOX_24T_DEGs %>% 
  dplyr::filter(., adj.P.Val < 0.05)

DOX24T_DEGs_GO_up <- DOX24T_DEGs_GO %>% 
  dplyr::filter(., logFC > 0)  %>% 
  rownames_to_column(., var = "entrezgene_ID") %>% 
  dplyr::select("entrezgene_ID")
#has 4029 genes

DOX24T_DEGs_GO_down <- DOX24T_DEGs_GO %>% 
  dplyr::filter(., logFC < 0)  %>% 
  rownames_to_column(., var = "entrezgene_ID") %>% 
  dplyr::select("entrezgene_ID")
#has 4305 genes

#shared genes venn diagram
DOX24Tshare_DEGs_GO <- DOX_24T_shared_DEGs %>% 
  dplyr::filter(., adj.P.Val < 0.05) 

DOX24T_share_DEGs_GO_plot <- DOX24Tshare_DEGs_GO %>% 
  rownames_to_column(., var = "Entrez_ID") %>% 
  dplyr::select("Entrez_ID")

DOX24Tshare_DEGs_GO_up <- DOX24Tshare_DEGs_GO %>% 
  dplyr::filter(., logFC > 0)  %>% 
  rownames_to_column(., var = "entrezgene_ID") %>% 
  dplyr::select("entrezgene_ID")
#has 4731 genes

DOX24Tshare_DEGs_GO_down <- DOX24Tshare_DEGs_GO %>% 
  dplyr::filter(., logFC < 0)  %>% 
  rownames_to_column(., var = "entrezgene_ID") %>% 
  dplyr::select("entrezgene_ID")

DOX144Rshare_DEGs_GO <- DOX_144R_shared_DEGs %>% 
  dplyr::filter(., adj.P.Val < 0.05)

DOX144Rshare_DEGs_GO_plot <- DOX144Rshare_DEGs_GO %>% 
  rownames_to_column(., var = "Entrez_ID") %>% 
  dplyr::select("Entrez_ID")


#now go ahead and do this for each condition as well

#DOX24T
venn_DOX_24T <- getSet(object = venn_test, subset = c("Group 1"))
dim(venn_DOX_24T)

[1] 4362    2

#4362 genes in DOX_24T only

venn_DOX24T_DEGs <- venn_DOX_24T %>% 
  as.data.frame() %>% 
  column_to_rownames(., var = "Detail")

#pull these genes out of my DEG matrix with logFC
DOX_24T_DEGs <- Toptable_list$V.D24T[row.names(venn_DOX24T_DEGs),]

DOX24T_DEGs_GO <- DOX_24T_DEGs %>% 
  rownames_to_column(., var = "Entrez_ID") %>% 
  dplyr::filter(., adj.P.Val < 0.05) %>% 
  dplyr::select("Entrez_ID")

#DOX24R
venn_DOX_24R <- getSet(object = venn_test, subset = c("Group 2"))
dim(venn_DOX_24R)

[1] 2182    2

#2182 genes in DOX_24R only

venn_DOX24R_DEGs <- venn_DOX_24R %>% 
  as.data.frame() %>% 
  column_to_rownames(., var = "Detail")

#pull these genes out of my DEG matrix with logFC
DOX_24R_DEGs <- Toptable_list$V.D24R[row.names(venn_DOX24R_DEGs),]

DOX24R_DEGs_GO <- DOX_24R_DEGs %>% 
  rownames_to_column(., var = "Entrez_ID") %>% 
  dplyr::filter(., adj.P.Val < 0.05) %>% 
  dplyr::select("Entrez_ID")


#DOX144R
venn_DOX_144R <- getSet(object = venn_test, subset = c("Group 3"))
dim(venn_DOX_144R)

[1] 32  2

#32 genes in DOX_144R only

venn_DOX144R_DEGs <- venn_DOX_144R %>% 
  as.data.frame() %>% 
  column_to_rownames(., var = "Detail")

#pull these genes out of my DEG matrix with logFC
DOX_144R_DEGs <- Toptable_list$V.D144R[row.names(venn_DOX144R_DEGs),]

DOX144R_DEGs_GO <- DOX_144R_DEGs %>% 
  rownames_to_column(., var = "Entrez_ID") %>% 
  dplyr::filter(., adj.P.Val < 0.05) %>% 
  dplyr::select("Entrez_ID")

#this set doesn't yield a GO/KEGG plot as there are too few genes
#instead, I pulled all genes associated with DOX144R

##Venn Diagrams of DEG Overlap after RUVs Correction k=1

#plot a venn diagram with all of your conditions from your toptables

# Load DEGs Data
DOX_24T_1 <- read.csv("data/new/DEGs/TTBL_RUV_24T_OLD.csv")
DOX_24R_1 <- read.csv("data/new/DEGs/TTBL_RUV_24R_OLD.csv")
DOX_144R_1 <- read.csv("data/new/DEGs/TTBL_RUV_144R_OLD.csv")

# Extract Significant DEGs
DEG1_RUV <- DOX_24T_1$Entrez_ID[DOX_24T_1$adj.P.Val < 0.05]
DEG2_RUV <- DOX_24R_1$Entrez_ID[DOX_24R_1$adj.P.Val < 0.05]
DEG3_RUV <- DOX_144R_1$Entrez_ID[DOX_144R_1$adj.P.Val < 0.05]

venntest_RUV <- list(DEG1_RUV, DEG2_RUV, DEG3_RUV)
ggVennDiagram(
  venntest_RUV,
  category.names = c("DOX_24T", "DOX_24R", "DOX_144R")
) + ggtitle("DXR Specific and Shared DEGs RUVs")+
  theme(
    plot.title = element_text(size = 16, face = "bold"),  # Increase title size
    text = element_text(size = 16)  # Increase text size globally
  )

Past versions of "RUVs Overlap of DEGs-1.png"

#try and use the VennDetail package to extract the genes from each condition
#Set 1 - DOX_24T only genes
#Set 2 - DOX_24T shared genes
#Set 3 - DOX_144R shared genes
plot.new()
venn_test_RUV <- venndetail(venntest_RUV)
plot(venn_test_RUV)

Past versions of "RUVs Overlap of DEGs-2.png"

detail(venn_test_RUV)

=== Here is the detail of Venndiagram ===
Total results:  11760 x 2 
Total sets is: 7 
  Subset Detail
1 Shared 643837
2 Shared  60672
3 Shared  55450
4 Shared   1870
5 Shared  57822
6 Shared  84958
... with 11754 more rows ...

#now that I have the genes subsetted by condition, I can pull out the genes I want to look at

venn_DOX_24T_RUV <- getSet(object = venn_test_RUV, subset = c("Group 1"))
dim(venn_DOX_24T_RUV)

[1] 4411    2

#4411 genes in DOX_24T only

venn_DOX_shared_RUV <- getSet(object = venn_test_RUV, subset = c("Shared", "Group 1_Group 2", "Group 1_Group 3"))
dim(venn_DOX_shared_RUV)

[1] 5009    2

#total of 5009 genes in this set which does not include D24T specific genes

venn_DOX144R_shared_RUV <- getSet(object = venn_test_RUV, subset = c("Shared", "Group 3", "Group 2_Group 3", "Group 1_Group 3"))
dim(venn_DOX144R_shared_RUV)

[1] 660   2

#total of 660 genes: (original 509)
 #342 shared all
 #32 DOX144R Specific
 #85 DOX24T + DOX144R
 #181 DOX24R + DOX144R

#after that, run GO analysis
venn_shared_DEGs_RUV <- venn_DOX_shared_RUV %>% 
  as.data.frame() %>% 
  column_to_rownames(., var = "Detail")

venn_DOX24T_DEGs_RUV <- venn_DOX_24T_RUV %>% 
  as.data.frame() %>% 
  column_to_rownames(., var = "Detail")

venn_sharedD144R_DEGs_RUV <- venn_DOX144R_shared_RUV %>% 
  as.data.frame() %>% 
  column_to_rownames(., var = "Detail")

#pull these genes out of my DEG matrix with logFC
DOX_24T_DEGs_RUV <-
  Toptable_list_RUVk1_symbols$V.D24T_RUV[row.names(venn_DOX24T_DEGs_RUV),]
DOX_24T_shared_DEGs_RUV <- 
  Toptable_list_RUVk1_symbols$V.D24T_RUV[row.names(venn_shared_DEGs_RUV),]
DOX_144R_shared_DEGs_RUV <- 
  Toptable_list_RUVk1_symbols$V.D144R_RUV[row.names(venn_sharedD144R_DEGs_RUV),]

#next, make the vectors that I need for plotting after cutting off by adj. p value < 0.05
DOX24T_DEGs_GO_RUV <- DOX_24T_DEGs_RUV %>% 
  dplyr::filter(., adj.P.Val < 0.05)

#shared D24T genes venn diagram
DOX24Tshare_DEGs_GO_RUV <- DOX_24T_shared_DEGs_RUV %>% 
  dplyr::filter(., adj.P.Val < 0.05) 

DOX24T_share_DEGs_GO_plot_RUV <- DOX24Tshare_DEGs_GO_RUV %>% 
  dplyr::select("Entrez_ID")

#shared D144R genes venn diagram

DOX144Rshare_DEGs_GO_RUV <- DOX_144R_shared_DEGs_RUV %>% 
  dplyr::filter(., adj.P.Val < 0.05)

DOX144Rshare_DEGs_GO_plot_RUV <- DOX144Rshare_DEGs_GO_RUV %>% 
  dplyr::select("Entrez_ID")

#now go ahead and do this for each condition as well

#DOX24T
venn_DOX_24T_RUV <- getSet(object = venn_test_RUV, subset = c("Group 1"))
dim(venn_DOX_24T_RUV)

[1] 4411    2

#4411 genes in DOX_24T only

venn_DOX24T_DEGs_RUV <- venn_DOX_24T_RUV %>% 
  as.data.frame() %>% 
  column_to_rownames(., var = "Detail")

#pull these genes out of my DEG matrix with logFC
DOX_24T_DEGs_RUV <- Toptable_list_RUVk1_symbols$V.D24T_RUV[row.names(venn_DOX24T_DEGs_RUV),]

DOX24T_DEGs_GO_RUV <- DOX_24T_DEGs_RUV %>% 
  dplyr::filter(., adj.P.Val < 0.05) %>% 
  dplyr::select("Entrez_ID")

#DOX24R
venn_DOX_24R_RUV <- getSet(object = venn_test_RUV, subset = c("Group 2"))
dim(venn_DOX_24R_RUV)

[1] 2107    2

#2107 genes in DOX_24R only

venn_DOX24R_DEGs_RUV <- venn_DOX_24R_RUV %>% 
  as.data.frame() %>% 
  column_to_rownames(., var = "Detail")

#pull these genes out of my DEG matrix with logFC
DOX_24R_DEGs_RUV <- Toptable_list_RUVk1_symbols$V.D24R_RUV[row.names(venn_DOX24R_DEGs_RUV),]

DOX24R_DEGs_GO_RUV <- DOX_24R_DEGs_RUV %>% 
  dplyr::filter(., adj.P.Val < 0.05) %>% 
  dplyr::select("Entrez_ID")


#DOX144R
venn_DOX_144R_RUV <- getSet(object = venn_test_RUV, subset = c("Group 3"))
dim(venn_DOX_144R_RUV)

[1] 52  2

#52 genes in DOX_144R only

venn_DOX144R_DEGs_RUV <- venn_DOX_144R_RUV %>% 
  as.data.frame() %>% 
  column_to_rownames(., var = "Detail")

#pull these genes out of my DEG matrix with logFC
DOX_144R_DEGs_RUV <- 
  Toptable_list_RUVk1_symbols$V.D144R_RUV[row.names(venn_DOX144R_DEGs_RUV),]

DOX144R_DEGs_GO_RUV <- DOX_144R_DEGs_RUV %>% 
  dplyr::filter(., adj.P.Val < 0.05) %>% 
  dplyr::select("Entrez_ID")

###GO KEGG Overlap of DEGs From Original Data

#####DOX24 Genes#####
library(gprofiler2)

D24_DEGs_mat <- as.matrix(DOX24T_DEGs_GO)

DOX_24_dxr_gene <- gost(query = D24_DEGs_mat,
                      organism = "hsapiens",
                      ordered_query = FALSE,
                      measure_underrepresentation = FALSE,
                      evcodes = FALSE,
                      user_threshold = 0.05,
                      correction_method = c("fdr"),
                      sources = c("GO:BP", "KEGG"))

DOX_24_gost_genes <- gostplot(DOX_24_dxr_gene, capped = FALSE, interactive = TRUE)
DOX_24_gost_genes

table_DOX24_genes <- DOX_24_dxr_gene$result %>% 
  dplyr::select(c(source, term_id, term_name, intersection_size, term_size, p_value))

table_DOX24_genes %>% 
  mutate_at(.vars = 6, .funs = scales::label_scientific(digits=4)) %>% 
  kableExtra::kable(.,) %>% 
  kableExtra::kable_paper("striped", full_width = FALSE) %>% 
  kableExtra::kable_styling(full_width = FALSE, position = "left", bootstrap_options = c("striped", "hover")) %>% 
  kableExtra::scroll_box(width = "100%", height = "400px")