src/common_functions.R src/06_compute_tdd_index.R: add fille containing command function

src/06_compute_tdd_index.R

@@ -9,6 +9,9 @@ library(gridExtra)
 library(DHARMa)
 library(emmeans)
 
+
+source("src/common_functions.R")
+
 #' compute the TDD of every genes for the experiment treatment
 #'
 #' @param treatment The treatment for which we want to compute the TDD
@@ -99,10 +102,12 @@ create_density_rep <- function(rep, tdd_tibble, treatment = "BRAF") {
 #' @param tdd_tibble A tablle containing the TDD of each replicate
 #' @param treatment The treatment to which the tdd_tibblerefers
 #' @param output_folder The folder where the figure will be created
+#' @param ercc a bolean indicating if the count data was normalized using ercc
 #' @import tidyverse
 #' @import gridExtra
 create_density_figs <- function(tdd_tibble, treatment = "BRAF",
-                                output_folder = "./results/TDD_analysis") {
+                                output_folder = "./results/TDD_analysis",
+                                ercc = F) {
     dir.create(output_folder, showWarnings = F)
     rep <- tdd_tibble %>%
         select(-gene) %>%
@@ -117,7 +122,11 @@ create_density_figs <- function(tdd_tibble, treatment = "BRAF",
     my_plots[[length(my_plots) + 1]] <- create_density_rep(
         "mean", tdd_tibble, treatment
     )
-    pdf(paste0(output_folder, "/", treatment, "TDD_density.pdf"),
+    ercc_name <- ""
+    if (ercc == T) {
+        ercc_name = "_ercc"
+    }
+    pdf(paste0(output_folder, "/", treatment, ercc_name, "TDD_density.pdf"),
         height = 8, width = 12
     )
     print(do.call("grid.arrange", c(my_plots, ncol = 2)))
@@ -143,7 +152,7 @@ get_tdd_table <- function(treatment,
                           output_folder = "./results/TDD_analysis",
                           ercc = F) {
     tdd_table <- compute_tdd(treatment, mean_gene_count, ercc)
-    create_density_figs(tdd_table, treatment, output_folder)
+    create_density_figs(tdd_table, treatment, output_folder, ercc)
     write_tsv(
         tdd_table,
         paste0(output_folder, "/", treatment, "_TDD_table.txt")
@@ -356,7 +365,7 @@ create_figures <- function(mean_gene_count = 10, ercc = F,
 t_test_func <- function(row) {
     x1 <- as.numeric(row[1:3])
     x2 <- as.numeric(row[4:6])
-    test <- t.test(x1, x2)
+    test <- t.test(x1, x2, paired = T)
     return(test$p.value) # nolint
 }
 
@@ -389,21 +398,6 @@ get_cor_value <- function(tdd_full, my_group = "BRAF_DOWN") {
 }
 
 
-#' Add a column group to a table containing a gene column
-#'
-#' @param tdd_full a table of tdd with a gene columne
-#' @return The table with a column indicating if the genes are up, down or not
-#' regulated by BRAF
-get_group_columns <- function(tdd_full) {
-    down_braf <- read.table(paste0(dir, "BRAF_DOWN_1.5_name.txt"))$V1
-    up_braf <- read.table(paste0(dir, "BRAF_UP_1.5_name.txt"))$V1
-    tdd_full$group <- rep("CTRL", times = nrow(tdd_full))
-    tdd_full[tdd_full$gene %in% down_braf, "group"] <- "BRAF_DOWN"
-    tdd_full[tdd_full$gene %in% up_braf, "group"] <- "BRAF_UP"
-    return(tdd_full)
-}
-
-
 #' Prepare a dataframe used to compute the correlation between BRAF and DMSO TDD
 #'
 #' @param mean_gene_count The minimum mean gene counts to compute their TDD
@@ -439,6 +433,14 @@ prepare_df_for_cor <- function(mean_gene_count = 10, ercc = F,
         filter(group != "CTRL") %>%
         arrange(kind)
     write_tsv(tdd_full, paste0(output_folder, "/TDD_correlation_table.txt")) # nolint
+    tdd_full %>% 
+        filter(kind == "BRAF_DOWN_sig", mean_BRAF_TDD > mean_DMSO_TDD) %>%
+        select(gene) %>% 
+        write_tsv(paste0(output_folder, "/TDD_BRAF_BRAF_DOWN.txt"))
+    tdd_full %>% 
+        filter(kind == "BRAF_UP_sig", mean_DMSO_TDD > mean_BRAF_TDD ) %>%
+        select(gene) %>% 
+        write_tsv(paste0(output_folder, "/TDD_DMSO_BRAF_UP.txt"))
     return(tdd_full)
 }
 
@@ -479,11 +481,14 @@ produce_tdd_cor_figure <- function(mean_gene_count = 10, ercc = F,
         ))
     p_down <- ggplot(tdd_full %>%
         filter(group %in% c("BRAF_DOWN", "BRAF_DOWN_sig")),
-    mapping = aes(x = mean_DMSO_TDD, y = mean_BRAF_TDD, color = kind)
+    mapping = aes(x = mean_DMSO_TDD, y = mean_BRAF_TDD, color = kind, label = gene)
     ) +
         geom_point(aes(fill = kind), colour = "white", pch = 21, size = 3) +
         scale_fill_manual(values = c("black", "red")) +
         coord_cartesian(xlim = c(0, 1), ylim = c(0, 1)) +
+        geom_label_repel(data = tdd_full %>%
+        filter(kind == "BRAF_DOWN_sig"), size = 1,
+            mapping = aes(mean_DMSO_TDD, mean_BRAF_TDD, label = gene)) +
         geom_abline(
             intercept = downv$intercept, slope = downv$slope,
             color = "#4169e1"
@@ -502,6 +507,9 @@ produce_tdd_cor_figure <- function(mean_gene_count = 10, ercc = F,
         geom_point(aes(fill = kind), colour = "white", pch = 21, size = 3) +
         scale_fill_manual(values = c("black", "red")) +
         coord_cartesian(xlim = c(0, 1), ylim = c(0, 1)) +
+        geom_label_repel(data = tdd_full %>%
+        filter(kind == "BRAF_UP_sig"), size = 1,
+            mapping = aes(mean_DMSO_TDD, mean_BRAF_TDD, label = gene)) +
         geom_abline(
             intercept = downv$intercept,
             slope = downv$slope, color = "#CD5C5C"
@@ -637,7 +645,8 @@ create_pic_figs <- function(tdd_full, gene_type = "BRAF_DOWN",
 #' @param mean_gene_count The minimum gene count of a gene to be analysed
 #' @param ercc Indicate if the normalisation with ercc should be used to compute
 #' the TDD index
-get_tdd_pic_table <- function(peak_file, mean_gene_count = 10, ercc = F) {
+get_tdd_pic_table <- function(peak_file, mean_gene_count = 10, ercc = F,
+                              output_folder) {
     peak_list <- read.table(peak_file)$V1
     peak_table <- tibble(
         gene = names(table(peak_list)),
@@ -646,14 +655,14 @@ get_tdd_pic_table <- function(peak_file, mean_gene_count = 10, ercc = F) {
     tdd_full <- get_final_tdd_table(mean_gene_count, ercc)
     tdd_full <- tdd_full[, 
         c("gene", grep("X...", colnames(tdd_full), value = T))]
+    tdd_braf_braf_down <- read_tsv(paste0(output_folder, 
+        "/TDD_BRAF_BRAF_DOWN.txt"))
+    tdd_dmso_braf_up <- read_tsv(paste0(output_folder, 
+        "/TDD_DMSO_BRAF_UP.txt"))
     tdd_full <- tdd_full %>%
         mutate(
-            TDD_BRAF = X276_BRAF_TDD - X276_DMSO_TDD > 0.2 &
-                X277_BRAF_TDD - X277_DMSO_TDD > 0.2 &
-                X278_BRAF_TDD - X278_DMSO_TDD > 0.2,
-            TDD_DMSO = X276_DMSO_TDD - X276_BRAF_TDD > 0.2 &
-                X277_DMSO_TDD - X277_BRAF_TDD > 0.2 &
-                X278_DMSO_TDD - X278_BRAF_TDD > 0.2
+            TDD_BRAF = gene %in% tdd_braf_braf_down$gene,
+            TDD_DMSO = gene %in% tdd_dmso_braf_up$gene
         ) %>%
         select(gene, TDD_BRAF, TDD_DMSO) #nolint
     tdd_full <- tdd_full %>% left_join(peak_table) #nolint
@@ -675,7 +684,7 @@ build_peaks_fig <- function(peak_file, mean_gene_count = 10, ercc = F,
                             gene_type = "BRAF_DOWN",
                             output_folder = "./results/TDD_analysis",
                             peak_type = "RiboRNApeak") {
-    tdd_table <- get_tdd_pic_table(peak_file, mean_gene_count, ercc)
+    tdd_table <- get_tdd_pic_table(peak_file, mean_gene_count, ercc, output_folder)
     tdd_table <- get_group_columns(tdd_table)
     create_pic_figs(tdd_table, gene_type, output_folder, peak_type = peak_type)
 }
@@ -685,8 +694,8 @@ create_figures(mean_gene_count = 10, ercc = F, output_folder = "./results/TDD_an
 create_figures(mean_gene_count = 20, ercc = F, output_folder = "./results/TDD_analysis")
 create_figures(mean_gene_count = 50, ercc = F, output_folder = "./results/TDD_analysis")
 create_figures(mean_gene_count = 0, ercc = T, output_folder = "./results/TDD_analysis")
+create_figures(mean_gene_count = 10, ercc = T, output_folder = "./results/TDD_analysis")
 produce_tdd_cor_figure()
-build_peaks_fig()
 dir <- "./data/gene_lists/"
 files <- c("BRAF_vs_DMSO_RepMean.merged.bed", "DMSO_vs_BRAF_RepMean.merged.bed",
            "BRAF_sorted_26-40.merged.bed", "DMSO_sorted_26-40.merged.bed")

src/common_functions.R

+#!/bin/Rscript
+
+
+royalblue <- "#4169e1"
+indianred <- "#CD5C5C"
+
+
+#' Add a column group to a table containing a gene column
+#'
+#' @param tdd_full a table of tdd with a gene columne
+#' @return The table with a column indicating if the genes are up, down or not
+#' regulated by BRAF
+get_group_columns <- function(tdd_full, dir = "./data/gene_lists/") {
+    down_braf <- read.table(paste0(dir, "BRAF_DOWN_1.5_name.txt"))$V1
+    up_braf <- read.table(paste0(dir, "BRAF_UP_1.5_name.txt"))$V1
+    tdd_full$group <- rep("CTRL", times = nrow(tdd_full))
+    tdd_full[tdd_full$gene %in% down_braf, "group"] <- "BRAF_DOWN"
+    tdd_full[tdd_full$gene %in% up_braf, "group"] <- "BRAF_UP"
+    return(tdd_full)
+}