From 723a99a9016922fb59a1f0c670c5d1a27f750385 Mon Sep 17 00:00:00 2001
From: Fontrodona Nicolas <nicolas.fontrodona@ens-lyon.fr>
Date: Fri, 30 Sep 2022 11:21:36 +0200
Subject: [PATCH] src/06_compute_tdd_index.R: changes in some figure functions
---
src/06_compute_tdd_index.R | 278 +++++++++++++++----------------------
src/07_ribosome_density.py | 2 +-
2 files changed, 115 insertions(+), 165 deletions(-)
diff --git a/src/06_compute_tdd_index.R b/src/06_compute_tdd_index.R
index f86cef0..2f9a82a 100644
--- a/src/06_compute_tdd_index.R
+++ b/src/06_compute_tdd_index.R
@@ -119,7 +119,7 @@ create_density_rep <- function(rep, tdd_tibble, treatment = "BRAF",
#' Create density figure for every replicate and the mean of replicate
#'
-#' @param tdd_tibble A tablle containing the TDD of each replicate
+#' @param tdd_tibble A table 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
@@ -435,20 +435,21 @@ t_test_func <- function(row) {
#' and BRAF condition + a group column indicating if the gene is up
#' of down regulated by BRAF
#' @param my_group The group of gene of interest to select
+#' @param my_peak the kind of peaks of interest
#' @return a list with the correlation value of TDD, its intercept and
#' slope
-get_cor_value <- function(tdd_full, my_group = "BRAF_DOWN") {
- tdd <- tdd_full %>% filter(group == my_group) ## nolint
+get_cor_value <- function(tdd_full, my_group = "BRAF_DOWN", my_peak = "TC_peaks") {
+ tdd <- tdd_full %>% filter(group == my_group, peak == my_peak) ## nolint
c <- cor(
- tdd %>% filter(group == my_group) %>% # nolint
+ tdd %>% filter(group == my_group, peak == my_peak) %>% # nolint
dplyr::select(mean_DMSO_TDD) %>% unlist(), # nolint
tdd %>%
- filter(group == my_group) %>% # nolint
+ filter(group == my_group, peak == my_peak) %>% # nolint
dplyr::select(mean_BRAF_TDD) %>%
unlist() # nolint
)
res <- coef(lm(mean_BRAF_TDD ~ mean_DMSO_TDD, tdd %>%
- filter(group == my_group)))
+ filter(group == my_group, peak == my_peak)))
return(list(cor = round(c, 3), intercept = res[1], slope = res[2]))
}
@@ -502,28 +503,34 @@ prepare_df_for_cor <- function(mean_gene_count = 10, ercc = F,
dir <- "./data/gene_lists/"
down_braf <- read.table(paste0(dir, "BRAF_DOWN_1.5.txt"))$V1
up_braf <- read.table(paste0(dir, "BRAF_UP_1.5.txt"))$V1
+ tc_peaks <- read.table(paste0(dir, "mRNADOWN_riboRNApeakBRAF.txt"))$V1
+ cc_peaks <- read.table(paste0(dir, "mRNAUP_riboRNApeakDMSO.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"
+ tdd_full$peak <- rep("CTRL", times = nrow(tdd_full))
+ tdd_full[tdd_full$gene %in% tc_peaks, "peak"] <- "DOWN_TC_peaks"
+ tdd_full[tdd_full$gene %in% cc_peaks, "peak"] <- "UP_CC_peaks"
tdd_full$sig <- rep("", times = nrow(tdd_full))
tdd_full[tdd_full$p_value <= 0.05, "sig"] <- "_sig"
- tdd_full <- tdd_full %>% mutate(kind = paste0(group, sig)) # nolint
+ tdd_full <- tdd_full %>% mutate(kind = paste0(peak, sig),
+ gkind = paste0(group, sig)) # nolint
tdd_full <- tdd_full %>%
arrange(kind)
write_tsv(tdd_full, paste0(output_folder, "/TDD_correlation_table.txt")) # nolint
mean_tdd <- tdd_full %>%
select(-time_point) %>%
- group_by(gene, group, sig, kind) %>%
+ group_by(gene, group, peak, sig, kind) %>%
summarize_all(mean) %>%
ungroup()
write_tsv(mean_tdd, paste0(output_folder, "/TDD_correlation_table_avg.txt")) # nolint
mean_tdd %>%
- filter(kind == "BRAF_DOWN_sig", mean_BRAF_TDD > mean_DMSO_TDD) %>%
+ filter(gkind == "BRAF_DOWN_sig", mean_BRAF_TDD > mean_DMSO_TDD) %>%
dplyr::select(gene) %>%
distinct(gene) %>%
write_tsv(paste0(output_folder, "/TDD_BRAF_BRAF_DOWN.txt"))
mean_tdd %>%
- filter(kind == "BRAF_UP_sig", mean_DMSO_TDD > mean_BRAF_TDD) %>%
+ filter(gkind == "BRAF_UP_sig", mean_DMSO_TDD > mean_BRAF_TDD) %>%
dplyr::select(gene) %>%
distinct(gene) %>%
write_tsv(paste0(output_folder, "/TDD_DMSO_BRAF_UP.txt"))
@@ -548,12 +555,12 @@ prepare_df_for_cor <- function(mean_gene_count = 10, ercc = F,
filter_on_tp <- function(tdd_table, tp = 5) {
if (tp == "3_5") {
tmp3 <- tdd_table %>% filter(time_point == 3) %>% # nolint
- dplyr::select(gene, mean_DMSO_TDD, mean_BRAF_TDD, group, kind)
+ dplyr::select(gene, mean_DMSO_TDD, mean_BRAF_TDD, group, kind, peak)
tmp5 <- tdd_table %>% filter(time_point == 5) %>% # nolint
- dplyr::select(gene, mean_DMSO_TDD, mean_BRAF_TDD, group, kind)
+ dplyr::select(gene, mean_DMSO_TDD, mean_BRAF_TDD, group, kind, peak)
tdd_full <- left_join(
tmp3, tmp5,
- by = c("gene", "group", "kind"),
+ by = c("gene", "group", "kind", "peak"),
suffix = c("_3", "_5")
) %>%
rowwise() %>%
@@ -576,93 +583,69 @@ filter_on_tp <- function(tdd_table, tp = 5) {
produce_tdd_cor_figure <- function(mean_gene_count = 10, ercc = F,
output_folder = "./results/TDD_analysis/correlation") {
tdd_complete <- prepare_df_for_cor(mean_gene_count, ercc, output_folder)
- for (tp in c(3, 5, "3_5")) {
+ for (tp in c(3, 5, "3_5")) {
tdd_full <- filter_on_tp(tdd_complete, tp)
- downv <- get_cor_value(tdd_full, my_group = "BRAF_DOWN")
- upv <- get_cor_value(tdd_full, my_group = "BRAF_UP")
- ctrlv <- get_cor_value(tdd_full, my_group = "CTRL")
- p <- ggplot(tdd_full, mapping = aes(
- x = mean_DMSO_TDD,
- y = mean_BRAF_TDD, color = kind
- )) +
- geom_point(aes(fill = kind), colour = "white", pch = 21, size = 3) +
- scale_fill_manual(values = c(
- "#4169e1", "darkblue",
- "#CD5C5C", "darkred",
- "dimgrey", "#2b2b2b"
- )) +
- coord_cartesian(xlim = c(0, 1), ylim = c(0, 1)) +
- geom_abline(
- intercept = downv$intercept, slope = downv$slope,
- color = "#4169e1"
- ) +
- geom_abline(
- intercept = upv$intercept, slope = upv$slope,
- color = "#CD5C5C"
- ) +
- ggtitle(paste0(
- "TDD of genes with BRAF or DMSO treatment\n",
- "(R gene BRAF_down = ", downv$cor, ")(R BRAF_up = ",
- upv$cor, ")"
- ))
+ downv <- get_cor_value(tdd_full, my_group = "BRAF_DOWN", my_peak = "DOWN_TC_peaks")
+ upv <- get_cor_value(tdd_full, my_group = "BRAF_UP", my_peak = "UP_CC_peaks")
p_down <- ggplot(tdd_full %>%
- filter(group %in% c("BRAF_DOWN")),
- mapping = aes(x = mean_DMSO_TDD, y = mean_BRAF_TDD, color = kind, label = gene)
+ filter(peak == "DOWN_TC_peaks"),
+ mapping = aes(x = mean_DMSO_TDD, y = mean_BRAF_TDD, color = kind)
) +
- geom_point(aes(fill = kind), colour = "white", pch = 21, size = 3) +
- scale_fill_manual(values = c("black", "red")) +
+ geom_point(aes(fill = kind), colour = "white", pch = 21, size = 5) +
+ scale_fill_manual(values = c("grey", "black")) +
coord_cartesian(xlim = c(0, 1), ylim = c(0, 1)) +
- geom_abline(
- intercept = downv$intercept, slope = downv$slope,
- color = "#4169e1"
+ # geom_abline(
+ # intercept = downv$intercept, slope = downv$slope,
+ # color = "#4169e1"
+ # ) +
+ theme_classic() +
+ theme(
+ text = element_text(size = 26),
+ panel.grid.major.x = element_line(),
+ panel.grid.major.y = element_line(),
) +
geom_abline(slope = 1, color = "black") +
ggtitle(paste0(
- "TDD of genes down-regulated by BRAF",
+ "TDD of genes down-regulated w/ TC peaks by BRAF",
"\n(R gene BRAF_down = ", downv$cor, ")"
))
- p_up <- ggplot(tdd_full %>% filter(group %in% c("BRAF_UP")),
+ p_up <- ggplot(tdd_full %>% filter(peak == "UP_CC_peaks"),
mapping = aes(
x = mean_DMSO_TDD, y = mean_BRAF_TDD,
color = kind
)
) +
- geom_point(aes(fill = kind), colour = "white", pch = 21, size = 3) +
- scale_fill_manual(values = c("black", "red")) +
+ geom_point(aes(fill = kind), colour = "white", pch = 21, size = 5) +
+ scale_fill_manual(values = c("grey", "black")) +
coord_cartesian(xlim = c(0, 1), ylim = c(0, 1)) +
- geom_abline(
- intercept = upv$intercept,
- slope = upv$slope, color = "#CD5C5C"
+ theme_classic() +
+ theme(
+ text = element_text(size = 26),
+ panel.grid.major.x = element_line(),
+ panel.grid.major.y = element_line(),
) +
+ # geom_abline(
+ # intercept = upv$intercept,
+ # slope = upv$slope, color = "#CD5C5C"
+ # ) +
geom_abline(slope = 1, color = "black") +
ggtitle(paste0(
- "TDD of genes up-regulated by BRAF",
+ "TDD of genes up-regulated by BRAF w/ CC peaks",
"\n(R gene BRAF_up = ", upv$cor, ")"
))
-
- p_ctrl <- ggplot(tdd_full %>% filter(group %in% c("CTRL")),
- mapping = aes(
- x = mean_DMSO_TDD, y = mean_BRAF_TDD,
- color = kind
- )
- ) +
- 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_abline(
- intercept = ctrlv$intercept,
- slope = ctrlv$slope, color = "dimgrey"
- ) +
- geom_abline(slope = 1, color = "black") +
- ggtitle(paste0(
- "TDD of CTRL genes ",
- "\n(R gene CTRL = ", ctrlv$cor, ")"
- ))
- pdf(paste0(output_folder, "/TDD_correlation_figures", tp, "h.pdf"),
- width = 12,
- height = 17
+ pdf(paste0(output_folder,
+ "/TDD_correlation_figures_downtc_", tp, "h.pdf"),
+ width = 15,
+ height = 10
+ )
+ print(p_down)
+ dev.off()
+ pdf(paste0(output_folder,
+ "/TDD_correlation_figures_upcc_", tp, "h.pdf"),
+ width = 15,
+ height = 10
)
- grid.arrange(p, p_down, p_up, p_ctrl, ncol = 1) # nolint
+ print(p_up)
dev.off()
}
}
@@ -688,7 +671,7 @@ glm_nb_pvalue <- function(tab, pic_col,
family="nbinom2") # nolint
pval <- sprintf(summary(mod)$coeff$cond[2, 4], fmt = '%#.2e')
} else {
- mod <- glmmTMB(log1p(peak_surface) ~ group, data = tab, ziformula = ~1,
+ mod <- glmmTMB(log1p(peak_surface) ~ group, data = tab, ziformula = ~1,
family = "gaussian") # nolint
pval <- sprintf(summary(mod)$coeff$cond[2, 4], fmt = '%#.2e')
}
@@ -716,7 +699,7 @@ glm_log_pvalue <- function(tab, pic_col,
output_diag <- paste0(output_folder, "/diagnotics")
dir.create(output_diag, showWarnings = F)
ptype_name <- tab$ptype[1]
- output_fig <- paste0(output_diag, "/glm_nb_", pic_col, "_",
+ output_fig <- paste0(output_diag, "/glm_nb_", pic_col, "_",
ptype_name, ".pdf")
pdf(output_fig, width = 12, height = 8)
plot(simulateResiduals(mod)) # nolint
@@ -757,15 +740,15 @@ create_peaks_barplots <- function(tdd_filt, pval, output_fig,
if (data_col == "peak") {
title <- paste0(
"Average proportion of BRAF/DMSO peaks in genes ", gene_type,
- " and control genes\n",
- "(p_BRAF_peak = ", pval$peak_braf, ", p_DMSO_peak = ",
+ " and control genes\n",
+ "(p_BRAF_peak = ", pval$peak_braf, ", p_DMSO_peak = ",
pval$peak_dmso, ")"
)
} else {
title <- paste0(
"Average surface of BRAF/DMSO peaks in genes ", gene_type,
"and control genes\n",
- "(p_BRAF_peak = ", pval$peak_braf, ", p_DMSO_peak = ",
+ "(p_BRAF_peak = ", pval$peak_braf, ", p_DMSO_peak = ",
pval$peak_dmso, ")"
)
}
@@ -788,8 +771,8 @@ create_peaks_barplots <- function(tdd_filt, pval, output_fig,
tdd_filt %>% filter(ptype == "DMSO"), has_col, dirname(output_fig)))
ntitle <- paste0(
"Average proportion of genes ", gene_type,
- " and control genes habing at least one BRAF/DMSO peaks\n",
- "(p_BRAF_peak = ", pval_log$pbraf, ", p_DMSO_peak = ",
+ " and control genes habing at least one BRAF/DMSO peaks\n",
+ "(p_BRAF_peak = ", pval_log$pbraf, ", p_DMSO_peak = ",
pval_log$pdmso, ")"
)
p2 <- ggplot(tdd_fig, mapping = aes_string(
@@ -839,7 +822,7 @@ create_peaks_density <- function(tdd_filt, pval, outfile, data_col,
if (data_col == "peak") {
title <- paste0(
xlab, " in BRAF condition in genes ", gene_type,
- " and control genes\n",
+ " and control genes\n",
"(p_BRAF_peak = ", pval$peak_braf, ")"
)
p1 <- ggplot(tdd_filt %>% filter(ptype == "BRAF"), mapping = aes(x = !!as.symbol(data_col))) +
@@ -852,7 +835,7 @@ create_peaks_density <- function(tdd_filt, pval, outfile, data_col,
labs(x = "Number of ribosome peaks")
title <- paste0(
xlab, " in DMSO condition in genes ", gene_type,
- " and control genes\n",
+ " and control genes\n",
"(p_DMSO_peak = ", pval$peak_dmso, ")"
)
p2 <- ggplot(tdd_filt %>% filter(ptype == "DMSO"), mapping = aes(x = !!as.symbol(data_col))) +
@@ -866,10 +849,10 @@ create_peaks_density <- function(tdd_filt, pval, outfile, data_col,
} else {
title <- paste0(
"Distribution of ", xlab, " in BRAF condition in genes ", gene_type,
- " and control genes\n",
+ " and control genes\n",
"(p_BRAF_peak = ", pval$peak_braf, ")"
)
- p1 <- ggplot(tdd_filt %>% filter(ptype == "BRAF"),
+ p1 <- ggplot(tdd_filt %>% filter(ptype == "BRAF"),
mapping = aes_string(x = xcol, fill = "group")) +
geom_density(mapping = aes(fill = group), bw = bw, alpha = 0.5) +
scale_fill_manual(values = valcol) +
@@ -877,12 +860,12 @@ create_peaks_density <- function(tdd_filt, pval, outfile, data_col,
labs(x = xlab)
title <- paste0(
"Distribution of ", xlab, " in DMSO condition in genes ", gene_type,
- " and control genes\n",
+ " and control genes\n",
"(p_DMSO_peak = ", pval$peak_dmso, ")"
)
- p2 <- ggplot(tdd_filt %>% filter(ptype == "DMSO"),
+ p2 <- ggplot(tdd_filt %>% filter(ptype == "DMSO"),
mapping = aes_string(x = xcol, fill = "group")) +
- geom_density(mapping = aes(fill = group), bw = bw, alpha = 0.5) +
+ geom_density(mapping = aes(fill = group), bw = bw, alpha = 0.5) +
scale_fill_manual(values = valcol) +
ggtitle(title) +
labs(x = xlab)
@@ -909,9 +892,9 @@ create_pic_figs <- function(tdd_full, gene_type = "TDD_BRAF",
res_col <- paste0("mean_", data_col)
has_col <- paste0("has_", data_col)
has_res_col <- paste0("mean_", has_col)
- pval <- list(peak_braf = glm_nb_pvalue(tdd_full %>% filter(ptype == "BRAF"),
+ pval <- list(peak_braf = glm_nb_pvalue(tdd_full %>% filter(ptype == "BRAF"),
data_col, output_folder),
- peak_dmso = glm_nb_pvalue(tdd_full %>% filter(ptype == "DMSO"),
+ peak_dmso = glm_nb_pvalue(tdd_full %>% filter(ptype == "DMSO"),
data_col, output_folder))
output_fig <- paste0(
output_folder, "/barplot_peak_analysis_", gene_type, "_",
@@ -942,7 +925,7 @@ get_peak_df <- function(kind = "peak", table_gene) {
stop(paste0("kind must be either 'peak' or 'peak_surface' not ", kind))
}
pfiles <- c(
- "./data/gene_lists/BRAF_vs_DMSO.merged.bed",
+ "./data/gene_lists/BRAF_vs_DMSO.merged.bed",
"./data/gene_lists/DMSO_vs_BRAF.merged.bed")
if (kind == "peak") {
@@ -987,14 +970,14 @@ get_peak_df <- function(kind = "peak", table_gene) {
#' @param ercc Indicate if the normalisation with ercc should be used to compute
#' the TDD index
get_tdd_pic_table <- function(mean_gene_count = 10, ercc = F,
- output_folder, gene_type = "TDD_BRAF",
+ output_folder, gene_type = "TDD_BRAF",
kind = "peak") {
output_folder <- paste0(output_folder, "/correlation")
tdd_full <- get_final_tdd_table(mean_gene_count, ercc)
tdd_full <- tdd_full %>%
filter(time_point == 5) %>%
select(gene)
-
+
tdd_gn <- read_tsv(paste0(
output_folder,
"/", gene_type, ".txt"
@@ -1033,105 +1016,72 @@ build_peaks_fig <- function(mean_gene_count = 10, ercc = F,
}
#' Create TDD violin figure
-#'
-#' @param tp The time step chosen
-#' @param group "up_down" to display TDD for gene UP and DOWN regulated by BRAF
-#' "peak" to display the genes down-regulated with a BRAF peak and the gene
-#' up regulated with a DMSO peak
-create_tdd_violin <- function(tp = 5, groupn = "up_down",
+create_tdd_violin <- function(kind = "UP_CC_peaks",
tdd_file = paste0(
"./results/TDD_analysis/correlation/",
"TDD_correlation_table.txt"
),
output = "./results/TDD_analysis/tdd_violin") {
tdd_table <- read_tsv(tdd_file) # nolint
- tdd_table <- filter_on_tp(tdd_table, tp)
- title <- "Average TDD index comparison between the BRAF and DMSO conditions"
- if (groupn != "up_down") {
- tdd_table$group <- NULL
- dir <- "./data/gene_lists/"
- gene_files <- c(
- paste0(dir, "mRNADOWN_riboRNApeakBRAF.txt"),
- paste0(dir, "mRNAUP_riboRNApeakDMSO.txt")
- )
- gene_names <- c("BRAF_down_riboRNApeak_BRAF", "BRAF_up_riboRNApeak_DMSO")
- for (i in seq_along(gene_files)) {
- gene_selected <- read_tsv(gene_files[i], col_names = F)$X1
- tdd_table[tdd_table$gene %in% gene_selected, "group"] <- gene_names[i]
- }
- tdd_table <- tdd_table %>% filter(!is.na(group))
- title <- paste0(
- title, " for genes up by BRAFi with a dmso peak\n",
- "and genes down by BRAFi and a BRAF peak"
- )
- tdd_table$group <- factor(tdd_table$group,
- level = c("BRAF_up_riboRNApeak_DMSO", "BRAF_down_riboRNApeak_BRAF")
- )
- my_range <- c(-0.5, 1)
- } else {
- tdd_table <- tdd_table %>% filter(group != "CTRL")
- title <- paste0(title, " for genes up and down by BRAFi")
- tdd_table$group <- factor(tdd_table$group,
- level = c("BRAF_UP", "BRAF_DOWN")
- )
- my_range <- c(-1, 2)
- }
+ tdd_table <- filter_on_tp(tdd_table, "3_5")
+ tdd_table[tdd_table$peak != kind, "peak"] <- "CTRL"
tdd_table <- tdd_table %>% # nolint
dplyr::select(
gene, mean_DMSO_TDD, mean_BRAF_TDD,
- mean_BRAF_TDD, group
+ mean_BRAF_TDD, peak
) %>%
rename(DMSO = mean_DMSO_TDD, BRAF = mean_BRAF_TDD) %>%
- pivot_longer(c(-gene, -group),
+ pivot_longer(c(-gene, -peak),
names_to = "condition",
values_to = "TDD"
)
- tdd_table$condition <- factor(tdd_table$condition,
- level = c("DMSO", "BRAF")
- )
+ if (kind == "UP_CC_peaks") {
+ tdd_table <- tdd_table %>% filter(condition == "DMSO")
+ title <- paste0("Average TDD index comparisons (at 3 and 5h) in DMSO condition\n",
+ "between TCpeak-UP genes in DMSO conditions and CTRL genes"
+ )
+ } else {
+ tdd_table <- tdd_table %>% filter(condition == "BRAF")
+ title <- paste0("Average TDD index comparisons (at 3 and 5h) in BRAF condition\n",
+ "between TCpeak-DOWN genes in BRAF conditions and CTRL genes"
+ )
+ }
+ tdd_table$peak <- factor(tdd_table$peak, level=c(kind, "CTRL"))
p <- ggplot(
data = tdd_table,
- mapping = aes(x = group, y = TDD, fill = condition)
+ mapping = aes(x = peak, y = TDD, fill = peak)
) +
- scale_fill_manual(values = c(indianred, royalblue)) +
+ scale_fill_manual(values = c("grey", "white")) +
geom_violin(position = position_dodge(width = 0.8)) +
geom_boxplot(width = 0.05, position = position_dodge(width = 0.8)) +
- coord_cartesian(ylim = my_range) +
- ggtitle(paste(
- "Average TDD index comparison between the BRAF and DMSO",
- "condition on gene up and down by BRAFi"
- ))
+ theme_classic() +
+ theme(
+ text = element_text(size = 26),
+ panel.grid.major.x = element_line(),
+ panel.grid.major.y = element_line(),
+ ) +
+ coord_cartesian(ylim = c(-0.5, 1.5)) +
+ ggtitle(title)
dir.create(output, showWarnings = F)
- outname <- paste0(
- "Violin_BRAF_DMSO_stat_", groupn, "_",
- tp, "h"
- )
+ outname <- paste0("Peak_", kind, "_3_5h")
pdf(paste0(output, "/", outname, ".pdf"), height = 10, width = 17)
print(p)
dev.off()
- fit <- lm(TDD^0.3 ~ condition * group, data = tdd_table)
+ fit <- lm(TDD ^ 0.5 ~ peak, data = tdd_table)
simulationOutput <- simulateResiduals(fit, n = 250, plot = F) # nolint
dir.create(paste0(output, "/diag"), showWarnings = F)
pdf(paste0(output, "/diag/", outname, ".pdf"), width = 17, height = 10)
plot(simulationOutput)
dev.off()
- tuk <- emmeans(fit, ~ condition * group) # nolint
- res <- summary(
- contrast(tuk, method = "pairwise", by = "group"),
- adjust = "fdr",
- by = NULL
- )
+ res <- as.data.frame(summary(fit)$coefficients)
+ res$group <- rownames(res)
write_tsv(res, paste0(output, "/", outname, ".txt")) # nolint
}
#' Launch create_tdd_violins for all possible values
create_tdd_violins <- function() {
- create_tdd_violin(5, "up_down")
- create_tdd_violin("3_5", "up_down")
- create_tdd_violin(3, "up_down")
- create_tdd_violin(3, "peak")
- create_tdd_violin(5, "peak")
- create_tdd_violin("3_5", "peak")
+ create_tdd_violin(kind = "UP_CC_peaks")
+ create_tdd_violin(kind = "DOWN_TC_peaks")
}
diff --git a/src/07_ribosome_density.py b/src/07_ribosome_density.py
index b33ad7a..bb413ad 100644
--- a/src/07_ribosome_density.py
+++ b/src/07_ribosome_density.py
@@ -1,4 +1,4 @@
-#!/usr/bine/env
+#!/usr/bin/env
# -*- condig: UTF-8 -*-
--
GitLab