# import the mergedOnly dataset, provided by Dr Jimmy Breen on the 24/09/21
rawCount <- read_tsv(here::here("0_data/raw_data/allSamples_mergedOnly.featureCounts.txt"),
                     col_names = TRUE,
                     comment = "#") %>%
  dplyr::rename(CONT1 = "../2_Hisat2_merged/CONT1_ATGTCA_merged.sorted.nodup.bam",
                CONT2 = "../2_Hisat2_merged/CONT2_CGATGT_merged.sorted.nodup.bam",
                CONT4 = "../2_Hisat2_merged/CONT4_ACTTGA_merged.sorted.nodup.bam",
                INT1 = "../2_Hisat2_merged/INT1_GTCCGC_merged.sorted.nodup.bam",
                INT2 = "../2_Hisat2_merged/INT2_ACAGTG_merged.sorted.nodup.bam",
                INT3 = "../2_Hisat2_merged/INT3_GATCAG_merged.sorted.nodup.bam",
                INT4 = "../2_Hisat2_merged/INT4_CTTGTA_merged.sorted.nodup.bam",
                SVX1 = "../2_Hisat2_merged/SVX1_GTTTCG_merged.sorted.nodup.bam",
                SVX2 = "../2_Hisat2_merged/SVX2_TAGCTT_merged.sorted.nodup.bam",
                SVX3 = "../2_Hisat2_merged/SVX3_ATCACG_merged.sorted.nodup.bam",
                SVX4 = "../2_Hisat2_merged/SVX4_GCCAAT_merged.sorted.nodup.bam",
                SVX_VAS1 = "../2_Hisat2_merged/SVX_VAS1_AGTCAA_merged.sorted.nodup.bam",
                SVX_VAS2 = "../2_Hisat2_merged/SVX_VAS2_AGTTCC_merged.sorted.nodup.bam",
                SVX_VAS3 = "../2_Hisat2_merged/SVX_VAS3_TGACCA_merged.sorted.nodup.bam",
                SVX_VAS4 = "../2_Hisat2_merged/SVX_VAS4_GGCTAC_merged.sorted.nodup.bam",
                VAS1 = "../2_Hisat2_merged/VAS1_CAGATC_merged.sorted.nodup.bam",
                VAS2 = "../2_Hisat2_merged/VAS2_GTGAAA_merged.sorted.nodup.bam",
                VAS3 = "../2_Hisat2_merged/VAS3_GTGGCC_merged.sorted.nodup.bam",
                VAS4 = "../2_Hisat2_merged/VAS4_CCGTCC_merged.sorted.nodup.bam",) %>% 
  column_to_rownames("Geneid") %>%
  as.data.frame() 
  rownames(rawCount) <- gsub("\\..+$", "", rownames(rawCount))
  
# Removing the non-numerical metadata column. SVX_VAS1 may also be an outlier, it is number 18 (BTW)
rawCount<- rawCount[, c(6,7,9:25)]

#Create DGElist with rawCOunt and gene data. Remove all genes with 0 expression in all treatment groups
dge <- DGEList(counts = rawCount, 
               samples = samples,
               genes = geneMetadata,
               remove.zeros = TRUE
               ) 

libSize <- dge$samples %>% 
  
  #plot the sample with the lib.size in x and sample_group in y, colour fill for each sample_group
  ggplot(aes(
    x = lib.size, 
    y = rownames(dge$samples), 
    fill = dge$samples$group)
    ) +
  geom_col() +
  
  #draw a vertical line for the mean lib.size
  geom_vline(
    aes (xintercept = lib.size), 
    data = . %>% summarise_at(vars(lib.size), mean),
    linetype = 2
  ) +
  
  #labelling splot
  labs(
    title = "Sample Library Size",
    x = "Library Size",
    y = "Samples",
    fill = "Sample Groups"
  ) +
  
  #PUBLISHING
  theme(legend.position = "none")

libSize

#save the plot to .svg
ggsave(here::here("2_plots/qc/library_size.svg"),
       plot = libSize + pub,
       #PUBLISHING
       width = 250,
       height = 166,
       units = "mm")

beforeFiltering <- dge %>% 
  
  #transform the raw count to logCPM 
  edgeR::cpm(log = TRUE) %>% 

  #melting (reorganising) the transformed logCPM data with respect to the id variable (i.e., the row and column names). Very similar to pivot_long function
  melt %>%

  #retain all rows where the logCPM (value) is finite. All of them in this case are finite
  dplyr::filter(is.finite(value)) %>%

  #plot the long formate logCPM counts in a density plot with x the logCPM and colour by the sample_id
  ggplot(aes(
    x = value,
    colour = Var2
  )) +
  geom_density() +

  # remove the legend
  guides(colour = FALSE) +

  #add figure tittle and subtitle and labels
  ggtitle("Before Filtering Low-expressed Genes", subtitle = paste0(nrow(dge), " genes"))+
  labs(
    x = "logCPM",
    y = "Density",
    colour = "Sample Groups"
  )

#save plot
ggsave("counts_before_filtering.svg",
       plot = beforeFiltering + pub,
       width = 250,
       height = 166,
       units = "mm",
       path = here::here("2_plots/qc/"))

#the genes kept have >2 CPM for at least 3 samples
keptGenes <- (rowSums(cpm(dge) > 3) >= 3) 
afterFiltering <- dge %>% 
  
  #transform the raw count to logCPM 
  edgeR::cpm(log = TRUE) %>% 
  
  #for var1 (gene names) extract only the keptGenes and discard all other genes in the logCPM data
  magrittr::extract(keptGenes,) %>% 

  #melting (reorganising) the transformed logCPM data with respect to the id variable (i.e., the row and column names). Very similar to pivot_long function
  melt %>% 
  
  #retain all rows where the logCPM (value) is finite. All of them in this case are finite
  dplyr::filter(is.finite(value)) %>% 

  #ggplot
  ggplot(aes(
    x = value,
    colour = Var2
  )) +
  geom_density() +
  
  #remove the legend
  guides(colour = FALSE) +
  
  #add figure tittle and subtitle and labels. since keptGenes is a logic element, the second element represents the number of genes that were kept after the filtering
  ggtitle("After Filtering Low-expressed Genes", subtitle = paste0(table(keptGenes)[[2]], " genes")) +
  labs(
    x = "logCPM",
    y = "Density",
    colour = "Sample Groups"
  )

#save plot
ggsave("counts_after_filtering_3_3.svg",
       plot = afterFiltering + pub,
       width = 250,
       height = 166,
       units = "mm",
       path = here::here("2_plots/qc/"))

#display plot
# afterFiltering

#display plot side by side
cowplot::plot_grid(beforeFiltering + pub, afterFiltering + pub)

ggsave(filename = "counts_before_after_filtering_3_3.svg",
        path = here::here("2_plots/qc/"),

       # PUBLISHING
       width = 320,
       height = 180,
       units = "mm")

#after normalisation
dge <- edgeR::calcNormFactors(object = dge, 
                              method = "TMM")
knitr::kable(dge$samples, caption = "Normalised samples") %>% 
  kable_styling(bootstrap_options = c("striped", "hover")) %>% 
  scroll_box(height = "600px")

# Save DGElist object into the data/R directory 
saveRDS(object = dge, file = here::here("0_data/RDS_objects/dge.rds"))
saveRDS(object = pub, file = here::here("0_data/RDS_objects/pub.rds"))

# saveRDS(object = gg_publish, file = here::here("0_data/RDS_objects/gg_publish.rds"))