calculate_interaction_zscores.R 53 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273747576777879808182838485868788899091929394959697989910010110210310410510610710810911011111211311411511611711811912012112212312412512612712812913013113213313413513613713813914014114214314414514614714814915015115215315415515615715815916016116216316416516616716816917017117217317417517617717817918018118218318418518618718818919019119219319419519619719819920020120220320420520620720820921021121221321421521621721821922022122222322422522622722822923023123223323423523623723823924024124224324424524624724824925025125225325425525625725825926026126226326426526626726826927027127227327427527627727827928028128228328428528628728828929029129229329429529629729829930030130230330430530630730830931031131231331431531631731831932032132232332432532632732832933033133233333433533633733833934034134234334434534634734834935035135235335435535635735835936036136236336436536636736836937037137237337437537637737837938038138238338438538638738838939039139239339439539639739839940040140240340440540640740840941041141241341441541641741841942042142242342442542642742842943043143243343443543643743843944044144244344444544644744844945045145245345445545645745845946046146246346446546646746846947047147247347447547647747847948048148248348448548648748848949049149249349449549649749849950050150250350450550650750850951051151251351451551651751851952052152252352452552652752852953053153253353453553653753853954054154254354454554654754854955055155255355455555655755855956056156256356456556656756856957057157257357457557657757857958058158258358458558658758858959059159259359459559659759859960060160260360460560660760860961061161261361461561661761861962062162262362462562662762862963063163263363463563663763863964064164264364464564664764864965065165265365465565665765865966066166266366466566666766866967067167267367467567667767867968068168268368468568668768868969069169269369469569669769869970070170270370470570670770870971071171271371471571671771871972072172272372472572672772872973073173273373473573673773873974074174274374474574674774874975075175275375475575675775875976076176276376476576676776876977077177277377477577677777877978078178278378478578678778878979079179279379479579679779879980080180280380480580680780880981081181281381481581681781881982082182282382482582682782882983083183283383483583683783883984084184284384484584684784884985085185285385485585685785885986086186286386486586686786886987087187287387487587687787887988088188288388488588688788888989089189289389489589689789889990090190290390490590690790890991091191291391491591691791891992092192292392492592692792892993093193293393493593693793893994094194294394494594694794894995095195295395495595695795895996096196296396496596696796896997097197297397497597697797897998098198298398498598698798898999099199299399499599699799899910001001100210031004100510061007100810091010101110121013101410151016101710181019102010211022102310241025102610271028102910301031103210331034103510361037103810391040104110421043104410451046104710481049105010511052105310541055105610571058105910601061106210631064106510661067106810691070107110721073107410751076107710781079108010811082108310841085108610871088108910901091109210931094109510961097109810991100110111021103110411051106110711081109111011111112111311141115111611171118111911201121112211231124112511261127112811291130113111321133113411351136113711381139114011411142114311441145114611471148114911501151115211531154115511561157115811591160116111621163116411651166116711681169117011711172117311741175117611771178117911801181118211831184118511861187118811891190119111921193119411951196119711981199120012011202120312041205120612071208120912101211121212131214121512161217121812191220122112221223122412251226122712281229123012311232123312341235123612371238123912401241124212431244124512461247124812491250125112521253125412551256125712581259126012611262126312641265126612671268126912701271127212731274127512761277127812791280128112821283128412851286128712881289129012911292129312941295129612971298129913001301130213031304130513061307130813091310131113121313131413151316131713181319132013211322132313241325132613271328132913301331133213331334133513361337133813391340134113421343134413451346134713481349135013511352135313541355135613571358135913601361136213631364136513661367136813691370137113721373137413751376137713781379138013811382138313841385138613871388138913901391139213931394139513961397139813991400140114021403140414051406
  1. suppressMessages({
  2. library("ggplot2")
  3. library("plotly")
  4. library("htmlwidgets")
  5. library("htmltools")
  6. library("dplyr")
  7. library("rlang")
  8. library("ggthemes")
  9. library("data.table")
  10. library("gridExtra")
  11. library("future")
  12. library("furrr")
  13. library("purrr")
  14. })
  15. # These parallelization libraries are very noisy
  16. suppressPackageStartupMessages({
  17. library("future")
  18. library("furrr")
  19. library("purrr")
  20. })
  21. options(warn = 2)
  22. # Constants for configuration
  23. plot_width <- 14
  24. plot_height <- 9
  25. base_size <- 14
  26. parse_arguments <- function() {
  27. args <- if (interactive()) {
  28. c(
  29. "/home/bryan/documents/develop/hartmanlab/qhtcp-workflow/out/20240116_jhartman2_DoxoHLD/20240116_jhartman2_DoxoHLD",
  30. "/home/bryan/documents/develop/hartmanlab/qhtcp-workflow/apps/r/SGD_features.tab",
  31. "/home/bryan/documents/develop/hartmanlab/qhtcp-workflow/out/20240116_jhartman2_DoxoHLD/easy/20240116_jhartman2_DoxoHLD/results_std.txt",
  32. "/home/bryan/documents/develop/hartmanlab/qhtcp-workflow/out/20240116_jhartman2_DoxoHLD/20240822_jhartman2_DoxoHLD/exp1",
  33. "Experiment 1: Doxo versus HLD",
  34. 3,
  35. "/home/bryan/documents/develop/hartmanlab/qhtcp-workflow/out/20240116_jhartman2_DoxoHLD/20240822_jhartman2_DoxoHLD/exp2",
  36. "Experiment 2: HLD versus Doxo",
  37. 3
  38. )
  39. } else {
  40. commandArgs(trailingOnly = TRUE)
  41. }
  42. out_dir <- normalizePath(args[1], mustWork = FALSE)
  43. sgd_gene_list <- normalizePath(args[2], mustWork = FALSE)
  44. easy_results_file <- normalizePath(args[3], mustWork = FALSE)
  45. # The remaining arguments should be in groups of 3
  46. exp_args <- args[-(1:3)]
  47. if (length(exp_args) %% 3 != 0) {
  48. stop("Experiment arguments should be in groups of 3: path, name, sd.")
  49. }
  50. # Extract the experiments into a list
  51. experiments <- list()
  52. for (i in seq(1, length(exp_args), by = 3)) {
  53. exp_name <- exp_args[i + 1]
  54. experiments[[exp_name]] <- list(
  55. path = normalizePath(exp_args[i], mustWork = FALSE),
  56. sd = as.numeric(exp_args[i + 2])
  57. )
  58. }
  59. # Extract the trailing number from each path
  60. trailing_numbers <- sapply(experiments, function(x) {
  61. path <- x$path
  62. nums <- gsub("[^0-9]", "", basename(path))
  63. as.integer(nums)
  64. })
  65. # Sort the experiments based on the trailing numbers
  66. sorted_experiments <- experiments[order(trailing_numbers)]
  67. list(
  68. out_dir = out_dir,
  69. sgd_gene_list = sgd_gene_list,
  70. easy_results_file = easy_results_file,
  71. experiments = sorted_experiments
  72. )
  73. }
  74. args <- parse_arguments()
  75. # Should we keep output in exp dirs or combine in the study output dir?
  76. # dir.create(file.path(args$out_dir, "zscores"), showWarnings = FALSE)
  77. # dir.create(file.path(args$out_dir, "zscores", "qc"), showWarnings = FALSE)
  78. theme_publication <- function(base_size = 14, base_family = "sans", legend_position = NULL) {
  79. # Ensure that legend_position has a valid value or default to "none"
  80. legend_position <- if (is.null(legend_position) || length(legend_position) == 0) "none" else legend_position
  81. theme_foundation <- ggthemes::theme_foundation(base_size = base_size, base_family = base_family)
  82. theme_foundation %+replace%
  83. theme(
  84. plot.title = element_text(face = "bold", size = rel(1.6), hjust = 0.5),
  85. text = element_text(),
  86. panel.background = element_blank(),
  87. plot.background = element_blank(),
  88. panel.border = element_blank(),
  89. axis.title = element_text(face = "bold", size = rel(1.4)),
  90. axis.title.y = element_text(angle = 90, vjust = 2),
  91. axis.text = element_text(size = rel(1.2)),
  92. axis.line = element_line(colour = "black"),
  93. panel.grid.major = element_line(colour = "#f0f0f0"),
  94. panel.grid.minor = element_blank(),
  95. legend.key = element_rect(colour = NA),
  96. legend.position = legend_position,
  97. legend.direction =
  98. if (legend_position == "right") {
  99. "vertical"
  100. } else if (legend_position == "bottom") {
  101. "horizontal"
  102. } else {
  103. NULL # No legend direction if position is "none" or other values
  104. },
  105. legend.spacing = unit(0, "cm"),
  106. legend.title = element_text(face = "italic", size = rel(1.3)),
  107. legend.text = element_text(size = rel(1.2)),
  108. plot.margin = unit(c(10, 5, 5, 5), "mm")
  109. )
  110. }
  111. scale_fill_publication <- function(...) {
  112. discrete_scale("fill", "Publication", manual_pal(values = c(
  113. "#386cb0", "#fdb462", "#7fc97f", "#ef3b2c", "#662506",
  114. "#a6cee3", "#fb9a99", "#984ea3", "#ffff33"
  115. )), ...)
  116. }
  117. scale_colour_publication <- function(...) {
  118. discrete_scale("colour", "Publication", manual_pal(values = c(
  119. "#386cb0", "#fdb462", "#7fc97f", "#ef3b2c", "#662506",
  120. "#a6cee3", "#fb9a99", "#984ea3", "#ffff33"
  121. )), ...)
  122. }
  123. # Load the initial dataframe from the easy_results_file
  124. load_and_filter_data <- function(easy_results_file, sd = 3) {
  125. df <- read.delim(easy_results_file, skip = 2, as.is = TRUE, row.names = 1, strip.white = TRUE)
  126. df <- df %>%
  127. filter(!(.[[1]] %in% c("", "Scan"))) %>%
  128. filter(!is.na(ORF) & ORF != "" & !Gene %in% c("BLANK", "Blank", "blank") & Drug != "BMH21") %>%
  129. # Rename columns
  130. rename(L = l, num = Num., AUC = AUC96, scan = Scan, last_bg = LstBackgrd, first_bg = X1stBackgrd) %>%
  131. mutate(
  132. across(c(Col, Row, num, L, K, r, scan, AUC, last_bg, first_bg), as.numeric),
  133. delta_bg = last_bg - first_bg,
  134. delta_bg_tolerance = mean(delta_bg, na.rm = TRUE) + (sd * sd(delta_bg, na.rm = TRUE)),
  135. NG = if_else(L == 0 & !is.na(L), 1, 0),
  136. DB = if_else(delta_bg >= delta_bg_tolerance, 1, 0),
  137. SM = 0,
  138. OrfRep = if_else(ORF == "YDL227C", "YDL227C", OrfRep), # should these be hardcoded?
  139. conc_num = as.numeric(gsub("[^0-9\\.]", "", Conc)),
  140. conc_num_factor = as.numeric(as.factor(conc_num)) - 1, # for legacy purposes
  141. conc_num_factor_factor = as.factor(conc_num)
  142. )
  143. return(df)
  144. }
  145. # Update Gene names using the SGD gene list
  146. update_gene_names <- function(df, sgd_gene_list) {
  147. # Load SGD gene list
  148. genes <- read.delim(file = sgd_gene_list,
  149. quote = "", header = FALSE,
  150. colClasses = c(rep("NULL", 3), rep("character", 2), rep("NULL", 11)))
  151. # Create a named vector for mapping ORF to GeneName
  152. gene_map <- setNames(genes$V5, genes$V4)
  153. # Vectorized match to find the GeneName from gene_map
  154. mapped_genes <- gene_map[df$ORF]
  155. # Replace NAs in mapped_genes with original Gene names (preserves existing Gene names if ORF is not found)
  156. updated_genes <- ifelse(is.na(mapped_genes) | df$OrfRep == "YDL227C", df$Gene, mapped_genes)
  157. # Ensure Gene is not left blank or incorrectly updated to "OCT1"
  158. df <- df %>%
  159. mutate(Gene = ifelse(updated_genes == "" | updated_genes == "OCT1", OrfRep, updated_genes))
  160. return(df)
  161. }
  162. calculate_summary_stats <- function(df, variables, group_vars) {
  163. summary_stats <- df %>%
  164. group_by(across(all_of(group_vars))) %>%
  165. summarise(
  166. N = n(),
  167. across(all_of(variables),
  168. list(
  169. mean = ~ mean(.x, na.rm = TRUE),
  170. median = ~ median(.x, na.rm = TRUE),
  171. max = ~ ifelse(all(is.na(.x)), NA, max(.x, na.rm = TRUE)),
  172. min = ~ ifelse(all(is.na(.x)), NA, min(.x, na.rm = TRUE)),
  173. sd = ~ sd(.x, na.rm = TRUE),
  174. se = ~ sd(.x, na.rm = TRUE) / sqrt(n() - 1)
  175. ),
  176. .names = "{.fn}_{.col}"
  177. ),
  178. .groups = "drop"
  179. )
  180. # Create a cleaned version of df that doesn't overlap with summary_stats
  181. cleaned_df <- df %>%
  182. select(-any_of(setdiff(intersect(names(df), names(summary_stats)), group_vars)))
  183. df_joined <- left_join(cleaned_df, summary_stats, by = group_vars)
  184. return(list(summary_stats = summary_stats, df_with_stats = df_joined))
  185. }
  186. calculate_interaction_scores <- function(df, max_conc, bg_stats,
  187. group_vars = c("OrfRep", "Gene", "num")) {
  188. # Calculate total concentration variables
  189. total_conc_num <- length(unique(df$conc_num))
  190. calculations <- df %>%
  191. group_by(across(all_of(group_vars))) %>%
  192. mutate(
  193. NG = sum(NG, na.rm = TRUE),
  194. DB = sum(DB, na.rm = TRUE),
  195. SM = sum(SM, na.rm = TRUE),
  196. num_non_removed_concs = total_conc_num - sum(DB, na.rm = TRUE) - 1,
  197. # Calculate raw data
  198. Raw_Shift_L = first(mean_L) - bg_stats$mean_L,
  199. Raw_Shift_K = first(mean_K) - bg_stats$mean_K,
  200. Raw_Shift_r = first(mean_r) - bg_stats$mean_r,
  201. Raw_Shift_AUC = first(mean_AUC) - bg_stats$mean_AUC,
  202. Z_Shift_L = first(Raw_Shift_L) / bg_stats$sd_L,
  203. Z_Shift_K = first(Raw_Shift_K) / bg_stats$sd_K,
  204. Z_Shift_r = first(Raw_Shift_r) / bg_stats$sd_r,
  205. Z_Shift_AUC = first(Raw_Shift_AUC) / bg_stats$sd_AUC,
  206. Exp_L = WT_L + Raw_Shift_L,
  207. Exp_K = WT_K + Raw_Shift_K,
  208. Exp_r = WT_r + Raw_Shift_r,
  209. Exp_AUC = WT_AUC + Raw_Shift_AUC,
  210. Delta_L = mean_L - Exp_L,
  211. Delta_K = mean_K - Exp_K,
  212. Delta_r = mean_r - Exp_r,
  213. Delta_AUC = mean_AUC - Exp_AUC,
  214. Delta_L = if_else(NG == 1, mean_L - WT_L, Delta_L),
  215. Delta_K = if_else(NG == 1, mean_K - WT_K, Delta_K),
  216. Delta_r = if_else(NG == 1, mean_r - WT_r, Delta_r),
  217. Delta_AUC = if_else(NG == 1, mean_AUC - WT_AUC, Delta_AUC),
  218. Delta_L = if_else(SM == 1, mean_L - WT_L, Delta_L),
  219. # Calculate Z-scores
  220. Zscore_L = Delta_L / WT_sd_L,
  221. Zscore_K = Delta_K / WT_sd_K,
  222. Zscore_r = Delta_r / WT_sd_r,
  223. Zscore_AUC = Delta_AUC / WT_sd_AUC,
  224. # Fit linear models and store in list-columns
  225. gene_lm_L = list(lm(Delta_L ~ conc_num_factor, data = pick(everything()))),
  226. gene_lm_K = list(lm(Delta_K ~ conc_num_factor, data = pick(everything()))),
  227. gene_lm_r = list(lm(Delta_r ~ conc_num_factor, data = pick(everything()))),
  228. gene_lm_AUC = list(lm(Delta_AUC ~ conc_num_factor, data = pick(everything()))),
  229. # Extract coefficients using purrr::map_dbl
  230. lm_intercept_L = map_dbl(gene_lm_L, ~ coef(.x)[1]),
  231. lm_slope_L = map_dbl(gene_lm_L, ~ coef(.x)[2]),
  232. lm_intercept_K = map_dbl(gene_lm_K, ~ coef(.x)[1]),
  233. lm_slope_K = map_dbl(gene_lm_K, ~ coef(.x)[2]),
  234. lm_intercept_r = map_dbl(gene_lm_r, ~ coef(.x)[1]),
  235. lm_slope_r = map_dbl(gene_lm_r, ~ coef(.x)[2]),
  236. lm_intercept_AUC = map_dbl(gene_lm_AUC, ~ coef(.x)[1]),
  237. lm_slope_AUC = map_dbl(gene_lm_AUC, ~ coef(.x)[2]),
  238. # Calculate lm_Score_* based on coefficients
  239. lm_Score_L = max_conc * lm_slope_L + lm_intercept_L,
  240. lm_Score_K = max_conc * lm_slope_K + lm_intercept_K,
  241. lm_Score_r = max_conc * lm_slope_r + lm_intercept_r,
  242. lm_Score_AUC = max_conc * lm_slope_AUC + lm_intercept_AUC,
  243. # Calculate R-squared values
  244. R_Squared_L = map_dbl(gene_lm_L, ~ summary(.x)$r.squared),
  245. R_Squared_K = map_dbl(gene_lm_K, ~ summary(.x)$r.squared),
  246. R_Squared_r = map_dbl(gene_lm_r, ~ summary(.x)$r.squared),
  247. R_Squared_AUC = map_dbl(gene_lm_AUC, ~ summary(.x)$r.squared)
  248. ) %>%
  249. ungroup()
  250. # Calculate overall mean and SD for lm_Score_* variables
  251. lm_means_sds <- calculations %>%
  252. summarise(
  253. lm_mean_L = mean(lm_Score_L, na.rm = TRUE),
  254. lm_sd_L = sd(lm_Score_L, na.rm = TRUE),
  255. lm_mean_K = mean(lm_Score_K, na.rm = TRUE),
  256. lm_sd_K = sd(lm_Score_K, na.rm = TRUE),
  257. lm_mean_r = mean(lm_Score_r, na.rm = TRUE),
  258. lm_sd_r = sd(lm_Score_r, na.rm = TRUE),
  259. lm_mean_AUC = mean(lm_Score_AUC, na.rm = TRUE),
  260. lm_sd_AUC = sd(lm_Score_AUC, na.rm = TRUE)
  261. )
  262. calculations <- calculations %>%
  263. mutate(
  264. Z_lm_L = (lm_Score_L - lm_means_sds$lm_mean_L) / lm_means_sds$lm_sd_L,
  265. Z_lm_K = (lm_Score_K - lm_means_sds$lm_mean_K) / lm_means_sds$lm_sd_K,
  266. Z_lm_r = (lm_Score_r - lm_means_sds$lm_mean_r) / lm_means_sds$lm_sd_r,
  267. Z_lm_AUC = (lm_Score_AUC - lm_means_sds$lm_mean_AUC) / lm_means_sds$lm_sd_AUC
  268. )
  269. # Summarize some of the stats
  270. interactions <- calculations %>%
  271. group_by(across(all_of(group_vars))) %>%
  272. mutate(
  273. # Calculate raw shifts
  274. Raw_Shift_L = first(Raw_Shift_L),
  275. Raw_Shift_K = first(Raw_Shift_K),
  276. Raw_Shift_r = first(Raw_Shift_r),
  277. Raw_Shift_AUC = first(Raw_Shift_AUC),
  278. # Calculate Z-shifts
  279. Z_Shift_L = first(Z_Shift_L),
  280. Z_Shift_K = first(Z_Shift_K),
  281. Z_Shift_r = first(Z_Shift_r),
  282. Z_Shift_AUC = first(Z_Shift_AUC),
  283. # Sum Z-scores
  284. Sum_Z_Score_L = sum(Zscore_L),
  285. Sum_Z_Score_K = sum(Zscore_K),
  286. Sum_Z_Score_r = sum(Zscore_r),
  287. Sum_Z_Score_AUC = sum(Zscore_AUC),
  288. # Calculate Average Z-scores
  289. Avg_Zscore_L = Sum_Z_Score_L / num_non_removed_concs,
  290. Avg_Zscore_K = Sum_Z_Score_K / num_non_removed_concs,
  291. Avg_Zscore_r = Sum_Z_Score_r / (total_conc_num - 1),
  292. Avg_Zscore_AUC = Sum_Z_Score_AUC / (total_conc_num - 1)
  293. ) %>%
  294. arrange(desc(Z_lm_L), desc(NG)) %>%
  295. ungroup()
  296. # Declare column order for output
  297. calculations <- calculations %>%
  298. select(
  299. "OrfRep", "Gene", "num", "N",
  300. "conc_num", "conc_num_factor", "conc_num_factor_factor",
  301. "mean_L", "mean_K", "mean_r", "mean_AUC",
  302. "median_L", "median_K", "median_r", "median_AUC",
  303. "sd_L", "sd_K", "sd_r", "sd_AUC",
  304. "se_L", "se_K", "se_r", "se_AUC",
  305. "Raw_Shift_L", "Raw_Shift_K", "Raw_Shift_r", "Raw_Shift_AUC",
  306. "Z_Shift_L", "Z_Shift_K", "Z_Shift_r", "Z_Shift_AUC",
  307. "WT_L", "WT_K", "WT_r", "WT_AUC",
  308. "WT_sd_L", "WT_sd_K", "WT_sd_r", "WT_sd_AUC",
  309. "Exp_L", "Exp_K", "Exp_r", "Exp_AUC",
  310. "Delta_L", "Delta_K", "Delta_r", "Delta_AUC",
  311. "Zscore_L", "Zscore_K", "Zscore_r", "Zscore_AUC",
  312. "NG", "SM", "DB"
  313. )
  314. interactions <- interactions %>%
  315. select(
  316. "OrfRep", "Gene", "num", "NG", "DB", "SM",
  317. "conc_num", "conc_num_factor", "conc_num_factor_factor",
  318. "Raw_Shift_L", "Raw_Shift_K", "Raw_Shift_r", "Raw_Shift_AUC",
  319. "Z_Shift_L", "Z_Shift_K", "Z_Shift_r", "Z_Shift_AUC",
  320. "Sum_Z_Score_L", "Sum_Z_Score_K", "Sum_Z_Score_r", "Sum_Z_Score_AUC",
  321. "Avg_Zscore_L", "Avg_Zscore_K", "Avg_Zscore_r", "Avg_Zscore_AUC",
  322. "lm_Score_L", "lm_Score_K", "lm_Score_r", "lm_Score_AUC",
  323. "R_Squared_L", "R_Squared_K", "R_Squared_r", "R_Squared_AUC",
  324. "Z_lm_L", "Z_lm_K", "Z_lm_r", "Z_lm_AUC"
  325. )
  326. # Clean the original dataframe by removing overlapping columns
  327. cleaned_df <- df %>%
  328. select(-any_of(
  329. setdiff(intersect(names(df), names(calculations)),
  330. c("OrfRep", "Gene", "num", "conc_num", "conc_num_factor", "conc_num_factor_factor"))))
  331. # Join the original dataframe with calculations
  332. df_with_calculations <- left_join(cleaned_df, calculations,
  333. by = c("OrfRep", "Gene", "num", "conc_num", "conc_num_factor", "conc_num_factor_factor"))
  334. # Remove overlapping columns between df_with_calculations and interactions
  335. df_with_calculations_clean <- df_with_calculations %>%
  336. select(-any_of(
  337. setdiff(intersect(names(df_with_calculations), names(interactions)),
  338. c("OrfRep", "Gene", "num", "conc_num", "conc_num_factor", "conc_num_factor_factor"))))
  339. # Join with interactions to create the full dataset
  340. full_data <- left_join(df_with_calculations_clean, interactions,
  341. by = c("OrfRep", "Gene", "num", "conc_num", "conc_num_factor", "conc_num_factor_factor"))
  342. return(list(
  343. calculations = calculations,
  344. interactions = interactions,
  345. full_data = full_data
  346. ))
  347. }
  348. generate_and_save_plots <- function(out_dir, filename, plot_configs) {
  349. message("Generating ", filename, ".pdf and ", filename, ".html")
  350. static_plots <- list()
  351. plotly_plots <- list()
  352. for (i in seq_along(plot_configs$plots)) {
  353. config <- plot_configs$plots[[i]]
  354. df <- config$df
  355. # Define aes_mapping, ensuring y_var is only used when it's not NULL
  356. aes_mapping <- switch(config$plot_type,
  357. "bar" = if (!is.null(config$color_var)) {
  358. aes(x = .data[[config$x_var]], fill = .data[[config$color_var]], color = .data[[config$color_var]])
  359. } else {
  360. aes(x = .data[[config$x_var]])
  361. },
  362. "density" = if (!is.null(config$color_var)) {
  363. aes(x = .data[[config$x_var]], color = .data[[config$color_var]])
  364. } else {
  365. aes(x = .data[[config$x_var]])
  366. },
  367. # For other plot types, only include y_var if it's not NULL
  368. if (!is.null(config$y_var) && !is.null(config$color_var)) {
  369. aes(x = .data[[config$x_var]], y = .data[[config$y_var]], color = .data[[config$color_var]])
  370. } else if (!is.null(config$y_var)) {
  371. aes(x = .data[[config$x_var]], y = .data[[config$y_var]])
  372. } else {
  373. aes(x = .data[[config$x_var]]) # no y_var needed for density and bar plots
  374. }
  375. )
  376. plot <- ggplot(df, aes_mapping) + theme_publication(legend_position = config$legend_position)
  377. # Apply appropriate plot function
  378. plot <- switch(config$plot_type,
  379. "scatter" = generate_scatter_plot(plot, config),
  380. "box" = generate_box_plot(plot, config),
  381. "density" = plot + geom_density(),
  382. "bar" = plot + geom_bar(),
  383. plot # default (unused)
  384. )
  385. # Add titles and labels
  386. if (!is.null(config$title)) plot <- plot + ggtitle(config$title)
  387. if (!is.null(config$x_label)) plot <- plot + xlab(config$x_label)
  388. if (!is.null(config$y_label)) plot <- plot + ylab(config$y_label)
  389. if (!is.null(config$coord_cartesian)) plot <- plot + coord_cartesian(ylim = config$coord_cartesian)
  390. # Convert ggplot to plotly, skipping subplot
  391. if (!is.null(config$tooltip_vars)) {
  392. plotly_plot <- suppressWarnings(plotly::ggplotly(plot, tooltip = config$tooltip_vars))
  393. } else {
  394. plotly_plot <- suppressWarnings(plotly::ggplotly(plot))
  395. }
  396. if (!is.null(plotly_plot[["frames"]])) {
  397. plotly_plot[["frames"]] <- NULL
  398. }
  399. # Adjust legend position in plotly
  400. if (!is.null(config$legend_position) && config$legend_position == "bottom") {
  401. plotly_plot <- plotly_plot %>% layout(legend = list(orientation = "h"))
  402. }
  403. static_plots[[i]] <- plot
  404. plotly_plots[[i]] <- plotly_plot
  405. }
  406. # Save static PDF plots
  407. pdf(file.path(out_dir, paste0(filename, ".pdf")), width = 16, height = 9)
  408. if (is.null(plot_configs$grid_layout)) {
  409. # Print each plot on a new page if grid_layout is not set
  410. for (plot in static_plots) {
  411. print(plot)
  412. }
  413. } else {
  414. # Use grid.arrange if grid_layout is set
  415. grid_nrow <- ifelse(is.null(plot_configs$grid_layout$nrow), length(plot_configs$plots), plot_configs$grid_layout$nrow)
  416. grid_ncol <- ifelse(is.null(plot_configs$grid_layout$ncol), 1, plot_configs$grid_layout$ncol)
  417. grid.arrange(grobs = static_plots, ncol = grid_ncol, nrow = grid_nrow)
  418. }
  419. dev.off()
  420. # Save combined HTML plot
  421. out_html_file <- file.path(out_dir, paste0(filename, ".html"))
  422. message("Saving combined HTML file: ", out_html_file)
  423. htmltools::save_html(
  424. htmltools::tagList(plotly_plots),
  425. file = out_html_file
  426. )
  427. }
  428. generate_scatter_plot <- function(plot, config) {
  429. # Define the points
  430. shape <- if (!is.null(config$shape)) config$shape else 3
  431. size <- if (!is.null(config$size)) config$size else 1.5
  432. position <-
  433. if (!is.null(config$position) && config$position == "jitter") {
  434. position_jitter(width = 0.1, height = 0)
  435. } else {
  436. "identity"
  437. }
  438. plot <- plot + geom_point(
  439. shape = shape,
  440. size = size,
  441. position = position
  442. )
  443. if (!is.null(config$cyan_points) && config$cyan_points) {
  444. plot <- plot + geom_point(
  445. aes(x = .data[[config$x_var]], y = .data[[config$y_var]]),
  446. color = "cyan",
  447. shape = 3,
  448. size = 0.5
  449. )
  450. }
  451. # Add Smooth Line if specified
  452. if (!is.null(config$smooth) && config$smooth) {
  453. smooth_color <- if (!is.null(config$smooth_color)) config$smooth_color else "blue"
  454. if (!is.null(config$lm_line)) {
  455. plot <- plot +
  456. geom_abline(
  457. intercept = config$lm_line$intercept,
  458. slope = config$lm_line$slope,
  459. color = smooth_color
  460. )
  461. } else {
  462. plot <- plot +
  463. geom_smooth(
  464. method = "lm",
  465. se = FALSE,
  466. color = smooth_color
  467. )
  468. }
  469. }
  470. # Add SD Bands if specified
  471. if (!is.null(config$sd_band)) {
  472. plot <- plot +
  473. annotate(
  474. "rect",
  475. xmin = -Inf, xmax = Inf,
  476. ymin = config$sd_band, ymax = Inf,
  477. fill = ifelse(!is.null(config$fill_positive), config$fill_positive, "#542788"),
  478. alpha = ifelse(!is.null(config$alpha_positive), config$alpha_positive, 0.3)
  479. ) +
  480. annotate(
  481. "rect",
  482. xmin = -Inf, xmax = Inf,
  483. ymin = -config$sd_band, ymax = -Inf,
  484. fill = ifelse(!is.null(config$fill_negative), config$fill_negative, "orange"),
  485. alpha = ifelse(!is.null(config$alpha_negative), config$alpha_negative, 0.3)
  486. ) +
  487. geom_hline(
  488. yintercept = c(-config$sd_band, config$sd_band),
  489. color = ifelse(!is.null(config$hl_color), config$hl_color, "gray")
  490. )
  491. }
  492. # Add Rectangles if specified
  493. if (!is.null(config$rectangles)) {
  494. for (rect in config$rectangles) {
  495. plot <- plot + annotate(
  496. "rect",
  497. xmin = rect$xmin,
  498. xmax = rect$xmax,
  499. ymin = rect$ymin,
  500. ymax = rect$ymax,
  501. fill = ifelse(is.null(rect$fill), NA, rect$fill),
  502. color = ifelse(is.null(rect$color), "black", rect$color),
  503. alpha = ifelse(is.null(rect$alpha), 0.1, rect$alpha)
  504. )
  505. }
  506. }
  507. # Add error bars if specified
  508. if (!is.null(config$error_bar) && config$error_bar && !is.null(config$y_var)) {
  509. if (!is.null(config$error_bar_params)) {
  510. plot <- plot +
  511. geom_errorbar(
  512. aes(
  513. ymin = config$error_bar_params$ymin,
  514. ymax = config$error_bar_params$ymax
  515. ),
  516. alpha = 0.3,
  517. linewidth = 0.5
  518. )
  519. } else {
  520. y_mean_col <- paste0("mean_", config$y_var)
  521. y_sd_col <- paste0("sd_", config$y_var)
  522. plot <- plot +
  523. geom_errorbar(
  524. aes(
  525. ymin = !!sym(y_mean_col) - !!sym(y_sd_col),
  526. ymax = !!sym(y_mean_col) + !!sym(y_sd_col)
  527. ),
  528. alpha = 0.3,
  529. linewidth = 0.5
  530. )
  531. }
  532. }
  533. # Customize X-axis if specified
  534. if (!is.null(config$x_breaks) && !is.null(config$x_labels) && !is.null(config$x_label)) {
  535. plot <- plot +
  536. scale_x_discrete(
  537. name = config$x_label,
  538. breaks = config$x_breaks,
  539. labels = config$x_labels
  540. )
  541. }
  542. # Set Y-axis limits if specified
  543. if (!is.null(config$ylim_vals)) {
  544. plot <- plot + scale_y_continuous(limits = config$ylim_vals)
  545. }
  546. # Add annotations if specified
  547. if (!is.null(config$annotations)) {
  548. for (annotation in config$annotations) {
  549. plot <- plot +
  550. annotate(
  551. "text",
  552. x = annotation$x,
  553. y = annotation$y,
  554. label = annotation$label,
  555. hjust = ifelse(is.null(annotation$hjust), 0.5, annotation$hjust),
  556. vjust = ifelse(is.null(annotation$vjust), 0.5, annotation$vjust),
  557. size = ifelse(is.null(annotation$size), 6, annotation$size),
  558. color = ifelse(is.null(annotation$color), "black", annotation$color)
  559. )
  560. }
  561. }
  562. return(plot)
  563. }
  564. generate_box_plot <- function(plot, config) {
  565. # Convert x_var to a factor within aes mapping
  566. plot <- plot + geom_boxplot(aes(x = factor(.data[[config$x_var]])))
  567. # Apply scale_x_discrete for breaks, labels, and axis label if provided
  568. if (!is.null(config$x_breaks) && !is.null(config$x_labels) && !is.null(config$x_label)) {
  569. plot <- plot + scale_x_discrete(
  570. name = config$x_label,
  571. breaks = config$x_breaks,
  572. labels = config$x_labels
  573. )
  574. }
  575. return(plot)
  576. }
  577. generate_plate_analysis_plot_configs <- function(variables, df_before = NULL, df_after = NULL,
  578. plot_type = "scatter", stages = c("before", "after")) {
  579. plots <- list()
  580. for (var in variables) {
  581. for (stage in stages) {
  582. df_plot <- if (stage == "before") df_before else df_after
  583. # Check for non-finite values in the y-variable
  584. df_plot_filtered <- df_plot %>% filter(is.finite(!!sym(var)))
  585. # Adjust settings based on plot_type
  586. config <- list(
  587. df = df_plot_filtered,
  588. x_var = "scan",
  589. y_var = var,
  590. plot_type = plot_type,
  591. title = paste("Plate analysis by Drug Conc for", var, stage, "quality control"),
  592. color_var = "conc_num_factor_factor",
  593. position = if (plot_type == "scatter") "jitter" else NULL,
  594. size = 0.2,
  595. error_bar = (plot_type == "scatter")
  596. )
  597. # Add config to plots list
  598. plots <- append(plots, list(config))
  599. }
  600. }
  601. return(list(grid_layout = list(ncol = 1, nrow = length(plots)), plots = plots))
  602. }
  603. generate_interaction_plot_configs <- function(df, limits_map = NULL, plot_type = "reference") {
  604. if (is.null(limits_map)) {
  605. limits_map <- list(
  606. L = c(0, 130),
  607. K = c(-20, 160),
  608. r = c(0, 1),
  609. AUC = c(0, 12500),
  610. Delta_L = c(-60, 60),
  611. Delta_K = c(-60, 60),
  612. Delta_r = c(-0.6, 0.6),
  613. Delta_AUC = c(-6000, 6000)
  614. )
  615. }
  616. group_vars <- if (plot_type == "reference") c("OrfRep", "Gene", "num") else c("OrfRep", "Gene")
  617. df_filtered <- df %>%
  618. mutate(OrfRepCombined = if (plot_type == "reference") paste(OrfRep, Gene, num, sep = "_") else paste(OrfRep, Gene, sep = "_"))
  619. # Separate the plots into two groups: overall variables and delta comparisons
  620. overall_plots <- list()
  621. delta_plots <- list()
  622. for (var in c("L", "K", "r", "AUC")) {
  623. y_limits <- limits_map[[var]]
  624. config <- list(
  625. df = df_filtered,
  626. plot_type = "scatter",
  627. x_var = "conc_num_factor_factor",
  628. y_var = var,
  629. x_label = unique(df_filtered$Drug)[1],
  630. title = sprintf("Scatter RF for %s with SD", var),
  631. coord_cartesian = y_limits,
  632. error_bar = TRUE,
  633. x_breaks = unique(df_filtered$conc_num_factor_factor),
  634. x_labels = as.character(unique(df_filtered$conc_num)),
  635. position = "jitter",
  636. smooth = TRUE
  637. )
  638. overall_plots <- append(overall_plots, list(config))
  639. }
  640. unique_groups <- df_filtered %>% select(all_of(group_vars)) %>% distinct()
  641. for (i in seq_len(nrow(unique_groups))) {
  642. group <- unique_groups[i, ]
  643. group_data <- df_filtered %>% semi_join(group, by = group_vars)
  644. OrfRep <- as.character(group$OrfRep)
  645. Gene <- if ("Gene" %in% names(group)) as.character(group$Gene) else ""
  646. num <- if ("num" %in% names(group)) as.character(group$num) else ""
  647. for (var in c("Delta_L", "Delta_K", "Delta_r", "Delta_AUC")) {
  648. y_limits <- limits_map[[var]]
  649. y_span <- y_limits[2] - y_limits[1]
  650. # Error bars
  651. WT_sd_var <- paste0("WT_sd_", sub("Delta_", "", var))
  652. WT_sd_value <- group_data[[WT_sd_var]][1]
  653. error_bar_ymin <- 0 - (2 * WT_sd_value)
  654. error_bar_ymax <- 0 + (2 * WT_sd_value)
  655. # Annotations
  656. Z_Shift_value <- round(group_data[[paste0("Z_Shift_", sub("Delta_", "", var))]][1], 2)
  657. Z_lm_value <- round(group_data[[paste0("Z_lm_", sub("Delta_", "", var))]][1], 2)
  658. NG_value <- group_data$NG[1]
  659. DB_value <- group_data$DB[1]
  660. SM_value <- group_data$SM[1]
  661. annotations <- list(
  662. list(x = 1, y = y_limits[2] - 0.2 * y_span, label = paste("ZShift =", Z_Shift_value)),
  663. list(x = 1, y = y_limits[2] - 0.3 * y_span, label = paste("lm ZScore =", Z_lm_value)),
  664. list(x = 1, y = y_limits[1] + 0.2 * y_span, label = paste("NG =", NG_value)),
  665. list(x = 1, y = y_limits[1] + 0.1 * y_span, label = paste("DB =", DB_value)),
  666. list(x = 1, y = y_limits[1], label = paste("SM =", SM_value))
  667. )
  668. config <- list(
  669. df = group_data,
  670. plot_type = "scatter",
  671. x_var = "conc_num_factor_factor",
  672. y_var = var,
  673. x_label = unique(group_data$Drug)[1],
  674. title = paste(OrfRep, Gene, sep = " "),
  675. coord_cartesian = y_limits,
  676. annotations = annotations,
  677. error_bar = TRUE,
  678. error_bar_params = list(
  679. ymin = error_bar_ymin,
  680. ymax = error_bar_ymax
  681. ),
  682. smooth = TRUE,
  683. x_breaks = unique(group_data$conc_num_factor_factor),
  684. x_labels = as.character(unique(group_data$conc_num)),
  685. ylim_vals = y_limits
  686. )
  687. delta_plots <- append(delta_plots, list(config))
  688. }
  689. }
  690. return(list(
  691. overall_plots = list(grid_layout = list(ncol = 2, nrow = 2), plots = overall_plots),
  692. delta_plots = list(grid_layout = list(ncol = 4, nrow = 3), plots = delta_plots)
  693. ))
  694. }
  695. generate_rank_plot_configs <- function(df, variables, is_lm = FALSE, adjust = FALSE, overlap_color = FALSE) {
  696. sd_bands <- c(1, 2, 3)
  697. configs <- list()
  698. # Adjust (if necessary) and rank columns
  699. for (variable in variables) {
  700. if (adjust) {
  701. df[[paste0("Avg_Zscore_", variable)]] <- ifelse(is.na(df[[paste0("Avg_Zscore_", variable)]]), 0.001, df[[paste0("Avg_Zscore_", variable)]])
  702. df[[paste0("Z_lm_", variable)]] <- ifelse(is.na(df[[paste0("Z_lm_", variable)]]), 0.001, df[[paste0("Z_lm_", variable)]])
  703. }
  704. df[[paste0("Rank_", variable)]] <- rank(df[[paste0("Avg_Zscore_", variable)]], na.last = "keep")
  705. df[[paste0("Rank_lm_", variable)]] <- rank(df[[paste0("Z_lm_", variable)]], na.last = "keep")
  706. }
  707. # Helper function to create a plot configuration
  708. create_plot_config <- function(variable, rank_var, zscore_var, y_label, sd_band, with_annotations = TRUE) {
  709. num_enhancers <- sum(df[[zscore_var]] >= sd_band, na.rm = TRUE)
  710. num_suppressors <- sum(df[[zscore_var]] <= -sd_band, na.rm = TRUE)
  711. # Default plot config
  712. plot_config <- list(
  713. df = df,
  714. x_var = rank_var,
  715. y_var = zscore_var,
  716. plot_type = "scatter",
  717. title = paste(y_label, "vs. Rank for", variable, "above", sd_band),
  718. sd_band = sd_band,
  719. fill_positive = "#542788",
  720. fill_negative = "orange",
  721. alpha_positive = 0.3,
  722. alpha_negative = 0.3,
  723. annotations = NULL,
  724. shape = 3,
  725. size = 0.1,
  726. y_label = y_label,
  727. x_label = "Rank",
  728. legend_position = "none"
  729. )
  730. if (with_annotations) {
  731. # Add specific annotations for plots with annotations
  732. plot_config$annotations <- list(
  733. list(
  734. x = median(df[[rank_var]], na.rm = TRUE),
  735. y = max(df[[zscore_var]], na.rm = TRUE) * 0.9,
  736. label = paste("Deletion Enhancers =", num_enhancers)
  737. ),
  738. list(
  739. x = median(df[[rank_var]], na.rm = TRUE),
  740. y = min(df[[zscore_var]], na.rm = TRUE) * 0.9,
  741. label = paste("Deletion Suppressors =", num_suppressors)
  742. )
  743. )
  744. }
  745. return(plot_config)
  746. }
  747. # Generate plots for each variable
  748. for (variable in variables) {
  749. rank_var <- if (is_lm) paste0("Rank_lm_", variable) else paste0("Rank_", variable)
  750. zscore_var <- if (is_lm) paste0("Z_lm_", variable) else paste0("Avg_Zscore_", variable)
  751. y_label <- if (is_lm) paste("Int Z score", variable) else paste("Avg Z score", variable)
  752. # Loop through SD bands
  753. for (sd_band in sd_bands) {
  754. # Create plot with annotations
  755. configs[[length(configs) + 1]] <- create_plot_config(variable, rank_var, zscore_var, y_label, sd_band, with_annotations = TRUE)
  756. # Create plot without annotations
  757. configs[[length(configs) + 1]] <- create_plot_config(variable, rank_var, zscore_var, y_label, sd_band, with_annotations = FALSE)
  758. }
  759. }
  760. return(list(grid_layout = list(ncol = 3, nrow = 2), plots = configs))
  761. }
  762. generate_correlation_plot_configs <- function(df, highlight_cyan = FALSE) {
  763. relationships <- list(
  764. list(x = "Z_lm_L", y = "Z_lm_K", label = "Interaction L vs. Interaction K"),
  765. list(x = "Z_lm_L", y = "Z_lm_r", label = "Interaction L vs. Interaction r"),
  766. list(x = "Z_lm_L", y = "Z_lm_AUC", label = "Interaction L vs. Interaction AUC"),
  767. list(x = "Z_lm_K", y = "Z_lm_r", label = "Interaction K vs. Interaction r"),
  768. list(x = "Z_lm_K", y = "Z_lm_AUC", label = "Interaction K vs. Interaction AUC"),
  769. list(x = "Z_lm_r", y = "Z_lm_AUC", label = "Interaction r vs. Interaction AUC")
  770. )
  771. plots <- list()
  772. for (rel in relationships) {
  773. lm_model <- lm(as.formula(paste(rel$y, "~", rel$x)), data = df)
  774. r_squared <- summary(lm_model)$r.squared
  775. config <- list(
  776. df = df,
  777. x_var = rel$x,
  778. y_var = rel$y,
  779. plot_type = "scatter",
  780. title = rel$label,
  781. annotations = list(
  782. list(x = mean(df[[rel$x]], na.rm = TRUE),
  783. y = mean(df[[rel$y]], na.rm = TRUE),
  784. label = paste("R-squared =", round(r_squared, 3)))
  785. ),
  786. smooth = TRUE,
  787. smooth_color = "tomato3",
  788. lm_line = list(intercept = coef(lm_model)[1], slope = coef(lm_model)[2]),
  789. shape = 3,
  790. size = 0.5,
  791. color_var = "Overlap",
  792. cyan_points = highlight_cyan
  793. )
  794. plots <- append(plots, list(config))
  795. }
  796. return(list(grid_layout = list(ncol = 3, nrow = 2), plots = plots))
  797. }
  798. main <- function() {
  799. lapply(names(args$experiments), function(exp_name) {
  800. exp <- args$experiments[[exp_name]]
  801. exp_path <- exp$path
  802. exp_sd <- exp$sd
  803. out_dir <- file.path(exp_path, "zscores")
  804. out_dir_qc <- file.path(exp_path, "zscores", "qc")
  805. dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
  806. dir.create(out_dir_qc, recursive = TRUE, showWarnings = FALSE)
  807. summary_vars <- c("L", "K", "r", "AUC", "delta_bg") # fields to filter and calculate summary stats across
  808. interaction_vars <- c("L", "K", "r", "AUC") # fields to calculate interaction z-scores
  809. message("Loading and filtering data for experiment: ", exp_name)
  810. df <- load_and_filter_data(args$easy_results_file, sd = exp_sd) %>%
  811. update_gene_names(args$sgd_gene_list) %>%
  812. as_tibble()
  813. # Filter rows above delta background tolerance
  814. df_above_tolerance <- df %>% filter(DB == 1)
  815. df_na <- df %>% mutate(across(all_of(summary_vars), ~ ifelse(DB == 1, NA, .))) # formerly X
  816. df_no_zeros <- df_na %>% filter(L > 0) # formerly X_noZero
  817. # Save some constants
  818. max_conc <- max(df$conc_num_factor)
  819. message("Calculating summary statistics before quality control")
  820. df_stats <- calculate_summary_stats(
  821. df = df,
  822. variables = summary_vars,
  823. group_vars = c("conc_num", "conc_num_factor", "conc_num_factor_factor"))$df_with_stats
  824. message("Calculating summary statistics after quality control")
  825. ss <- calculate_summary_stats(
  826. df = df_na,
  827. variables = summary_vars,
  828. group_vars = c("conc_num", "conc_num_factor", "conc_num_factor_factor"))
  829. df_na_ss <- ss$summary_stats
  830. df_na_stats <- ss$df_with_stats
  831. write.csv(df_na_ss, file = file.path(out_dir, "summary_stats_all_strains.csv"), row.names = FALSE)
  832. # For plotting (ggplot warns on NAs)
  833. df_na_stats_filtered <- df_na_stats %>% filter(if_all(all_of(summary_vars), is.finite))
  834. df_na_stats <- df_na_stats %>%
  835. mutate(
  836. WT_L = mean_L,
  837. WT_K = mean_K,
  838. WT_r = mean_r,
  839. WT_AUC = mean_AUC,
  840. WT_sd_L = sd_L,
  841. WT_sd_K = sd_K,
  842. WT_sd_r = sd_r,
  843. WT_sd_AUC = sd_AUC
  844. )
  845. # Pull the background means and standard deviations from zero concentration for interactions
  846. bg_stats <- df_na_stats %>%
  847. filter(conc_num == 0) %>%
  848. summarise(
  849. mean_L = first(mean_L),
  850. mean_K = first(mean_K),
  851. mean_r = first(mean_r),
  852. mean_AUC = first(mean_AUC),
  853. sd_L = first(sd_L),
  854. sd_K = first(sd_K),
  855. sd_r = first(sd_r),
  856. sd_AUC = first(sd_AUC)
  857. )
  858. message("Calculating summary statistics after quality control excluding zero values")
  859. df_no_zeros_stats <- calculate_summary_stats(
  860. df = df_no_zeros,
  861. variables = summary_vars,
  862. group_vars = c("conc_num", "conc_num_factor", "conc_num_factor_factor")
  863. )$df_with_stats
  864. message("Filtering by 2SD of K")
  865. df_na_within_2sd_k <- df_na_stats %>%
  866. filter(K >= (mean_K - 2 * sd_K) & K <= (mean_K + 2 * sd_K))
  867. df_na_outside_2sd_k <- df_na_stats %>%
  868. filter(K < (mean_K - 2 * sd_K) | K > (mean_K + 2 * sd_K))
  869. message("Calculating summary statistics for L within 2SD of K")
  870. # TODO We're omitting the original z_max calculation, not sure if needed?
  871. ss <- calculate_summary_stats(df_na_within_2sd_k, "L", group_vars = c("conc_num", "conc_num_factor", "conc_num_factor_factor"))$summary_stats
  872. write.csv(ss,
  873. file = file.path(out_dir_qc, "max_observed_L_vals_for_spots_within_2sd_K.csv"),
  874. row.names = FALSE)
  875. message("Calculating summary statistics for L outside 2SD of K")
  876. ss <- calculate_summary_stats(df_na_outside_2sd_k, "L", group_vars = c("conc_num", "conc_num_factor", "conc_num_factor_factor"))
  877. df_na_l_outside_2sd_k_stats <- ss$df_with_stats
  878. write.csv(ss$summary_stats,
  879. file = file.path(out_dir, "max_observed_L_vals_for_spots_outside_2sd_K.csv"),
  880. row.names = FALSE)
  881. # Each plots list corresponds to a file
  882. l_vs_k_plot_configs <- list(
  883. plots = list(
  884. list(
  885. df = df,
  886. x_var = "L",
  887. y_var = "K",
  888. plot_type = "scatter",
  889. tooltip_vars = c("OrfRep", "Gene", "delta_bg"),
  890. title = "Raw L vs K before quality control",
  891. color_var = "conc_num_factor_factor",
  892. error_bar = FALSE,
  893. legend_position = "right"
  894. )
  895. )
  896. )
  897. frequency_delta_bg_plot_configs <- list(
  898. plots = list(
  899. list(
  900. df = df_stats,
  901. x_var = "delta_bg",
  902. y_var = NULL,
  903. plot_type = "density",
  904. title = "Density plot for Delta Background by [Drug] (All Data)",
  905. color_var = "conc_num_factor_factor",
  906. x_label = "Delta Background",
  907. y_label = "Density",
  908. error_bar = FALSE,
  909. legend_position = "right"
  910. ),
  911. list(
  912. df = df_stats,
  913. x_var = "delta_bg",
  914. y_var = NULL,
  915. plot_type = "bar",
  916. title = "Bar plot for Delta Background by [Drug] (All Data)",
  917. color_var = "conc_num_factor_factor",
  918. x_label = "Delta Background",
  919. y_label = "Count",
  920. error_bar = FALSE,
  921. legend_position = "right"
  922. )
  923. )
  924. )
  925. above_threshold_plot_configs <- list(
  926. plots = list(
  927. list(
  928. df = df_above_tolerance,
  929. x_var = "L",
  930. y_var = "K",
  931. plot_type = "scatter",
  932. tooltip_vars = c("OrfRep", "Gene", "delta_bg"),
  933. title = paste("Raw L vs K for strains above Delta Background threshold of",
  934. round(df_above_tolerance$delta_bg_tolerance[[1]], 3), "or above"),
  935. color_var = "conc_num_factor_factor",
  936. position = "jitter",
  937. annotations = list(
  938. list(
  939. x = median(df_above_tolerance$L, na.rm = TRUE) / 2,
  940. y = median(df_above_tolerance$K, na.rm = TRUE) / 2,
  941. label = paste("# strains above Delta Background tolerance =", nrow(df_above_tolerance))
  942. )
  943. ),
  944. error_bar = FALSE,
  945. legend_position = "right"
  946. )
  947. )
  948. )
  949. plate_analysis_plot_configs <- generate_plate_analysis_plot_configs(
  950. variables = summary_vars,
  951. df_before = df_stats,
  952. df_after = df_na_stats_filtered
  953. )
  954. plate_analysis_boxplot_configs <- generate_plate_analysis_plot_configs(
  955. variables = summary_vars,
  956. df_before = df_stats,
  957. df_after = df_na_stats_filtered,
  958. plot_type = "box"
  959. )
  960. plate_analysis_no_zeros_plot_configs <- generate_plate_analysis_plot_configs(
  961. variables = summary_vars,
  962. stages = c("after"), # Only after QC
  963. df_after = df_no_zeros_stats
  964. )
  965. plate_analysis_no_zeros_boxplot_configs <- generate_plate_analysis_plot_configs(
  966. variables = summary_vars,
  967. stages = c("after"), # Only after QC
  968. df_after = df_no_zeros_stats,
  969. plot_type = "box"
  970. )
  971. l_outside_2sd_k_plot_configs <- list(
  972. plots = list(
  973. list(
  974. df = df_na_l_outside_2sd_k_stats,
  975. x_var = "L",
  976. y_var = "K",
  977. plot_type = "scatter",
  978. title = "Raw L vs K for strains falling outside 2SD of the K mean at each Conc",
  979. color_var = "conc_num_factor_factor",
  980. position = "jitter", # Apply jitter for better visibility
  981. tooltip_vars = c("OrfRep", "Gene", "delta_bg"),
  982. annotations = list(
  983. list(
  984. x = median(df_na_l_outside_2sd_k_stats$L, na.rm = TRUE) / 2,
  985. y = median(df_na_l_outside_2sd_k_stats$K, na.rm = TRUE) / 2,
  986. label = paste("Total strains:", nrow(df_na_l_outside_2sd_k_stats))
  987. )
  988. ),
  989. error_bar = FALSE,
  990. legend_position = "right"
  991. )
  992. )
  993. )
  994. delta_bg_outside_2sd_k_plot_configs <- list(
  995. plots = list(
  996. list(
  997. df = df_na_l_outside_2sd_k_stats,
  998. x_var = "delta_bg",
  999. y_var = "K",
  1000. plot_type = "scatter",
  1001. title = "Delta Background vs K for strains falling outside 2SD of the K mean at each Conc",
  1002. color_var = "conc_num_factor_factor",
  1003. position = "jitter", # Apply jitter for better visibility
  1004. tooltip_vars = c("OrfRep", "Gene", "delta_bg"),
  1005. annotations = list(
  1006. list(
  1007. x = median(df_na_l_outside_2sd_k_stats$delta_bg, na.rm = TRUE) / 2,
  1008. y = median(df_na_l_outside_2sd_k_stats$K, na.rm = TRUE) / 2,
  1009. label = paste("Total strains:", nrow(df_na_l_outside_2sd_k_stats))
  1010. )
  1011. ),
  1012. error_bar = FALSE,
  1013. legend_position = "right"
  1014. )
  1015. )
  1016. )
  1017. message("Generating quality control plots in parallel")
  1018. # # future::plan(future::multicore, workers = parallel::detectCores())
  1019. future::plan(future::multisession, workers = 3) # generate 3 plots in parallel
  1020. plot_configs <- list(
  1021. list(out_dir = out_dir_qc, filename = "L_vs_K_before_quality_control",
  1022. plot_configs = l_vs_k_plot_configs),
  1023. list(out_dir = out_dir_qc, filename = "frequency_delta_background",
  1024. plot_configs = frequency_delta_bg_plot_configs),
  1025. list(out_dir = out_dir_qc, filename = "L_vs_K_above_threshold",
  1026. plot_configs = above_threshold_plot_configs),
  1027. list(out_dir = out_dir_qc, filename = "plate_analysis",
  1028. plot_configs = plate_analysis_plot_configs),
  1029. list(out_dir = out_dir_qc, filename = "plate_analysis_boxplots",
  1030. plot_configs = plate_analysis_boxplot_configs),
  1031. list(out_dir = out_dir_qc, filename = "plate_analysis_no_zeros",
  1032. plot_configs = plate_analysis_no_zeros_plot_configs),
  1033. list(out_dir = out_dir_qc, filename = "plate_analysis_no_zeros_boxplots",
  1034. plot_configs = plate_analysis_no_zeros_boxplot_configs),
  1035. list(out_dir = out_dir_qc, filename = "L_vs_K_for_strains_2SD_outside_mean_K",
  1036. plot_configs = l_outside_2sd_k_plot_configs),
  1037. list(out_dir = out_dir_qc, filename = "delta_background_vs_K_for_strains_2sd_outside_mean_K",
  1038. plot_configs = delta_bg_outside_2sd_k_plot_configs)
  1039. )
  1040. # Generating quality control plots in parallel
  1041. furrr::future_map(plot_configs, function(config) {
  1042. generate_and_save_plots(config$out_dir, config$filename, config$plot_configs)
  1043. }, .options = furrr_options(seed = TRUE))
  1044. # Process background strains
  1045. bg_strains <- c("YDL227C")
  1046. lapply(bg_strains, function(strain) {
  1047. message("Processing background strain: ", strain)
  1048. # Handle missing data by setting zero values to NA
  1049. # and then removing any rows with NA in L col
  1050. df_bg <- df_na %>%
  1051. filter(OrfRep == strain) %>%
  1052. mutate(
  1053. L = if_else(L == 0, NA, L),
  1054. K = if_else(K == 0, NA, K),
  1055. r = if_else(r == 0, NA, r),
  1056. AUC = if_else(AUC == 0, NA, AUC)
  1057. ) %>%
  1058. filter(!is.na(L))
  1059. # Recalculate summary statistics for the background strain
  1060. message("Calculating summary statistics for background strain")
  1061. ss_bg <- calculate_summary_stats(df_bg, summary_vars, group_vars = c("OrfRep", "conc_num", "conc_num_factor", "conc_num_factor_factor"))
  1062. summary_stats_bg <- ss_bg$summary_stats
  1063. write.csv(summary_stats_bg,
  1064. file = file.path(out_dir, paste0("summary_stats_background_strain_", strain, ".csv")),
  1065. row.names = FALSE)
  1066. # Set the missing values to the highest theoretical value at each drug conc for L
  1067. # Leave other values as 0 for the max/min
  1068. df_reference <- df_na_stats %>% # formerly X2_RF
  1069. filter(OrfRep == strain) %>%
  1070. filter(!is.na(L)) %>%
  1071. group_by(conc_num, conc_num_factor, conc_num_factor_factor) %>%
  1072. mutate(
  1073. max_l_theoretical = max(max_L, na.rm = TRUE),
  1074. L = ifelse(L == 0 & !is.na(L) & conc_num > 0, max_l_theoretical, L),
  1075. SM = ifelse(L >= max_l_theoretical & !is.na(L) & conc_num > 0, 1, 0),
  1076. L = ifelse(L >= max_l_theoretical & !is.na(L) & conc_num > 0, max_l_theoretical, L)) %>%
  1077. ungroup()
  1078. # Ditto for deletion strains
  1079. df_deletion <- df_na_stats %>% # formerly X2
  1080. filter(OrfRep != strain) %>%
  1081. filter(!is.na(L)) %>%
  1082. mutate(SM = 0) %>%
  1083. group_by(conc_num, conc_num_factor, conc_num_factor_factor) %>%
  1084. mutate(
  1085. max_l_theoretical = max(max_L, na.rm = TRUE),
  1086. L = ifelse(L == 0 & !is.na(L) & conc_num > 0, max_l_theoretical, L),
  1087. SM = ifelse(L >= max_l_theoretical & !is.na(L) & conc_num > 0, 1, SM),
  1088. L = ifelse(L >= max_l_theoretical & !is.na(L) & conc_num > 0, max_l_theoretical, L)) %>%
  1089. ungroup()
  1090. message("Calculating reference strain interaction scores")
  1091. df_reference_stats <- calculate_summary_stats(
  1092. df = df_reference,
  1093. variables = interaction_vars,
  1094. group_vars = c("OrfRep", "Gene", "num", "conc_num", "conc_num_factor", "conc_num_factor_factor")
  1095. )$df_with_stats
  1096. reference_results <- calculate_interaction_scores(df_reference_stats, max_conc, bg_stats, group_vars = c("OrfRep", "Gene", "num"))
  1097. zscore_calculations_reference <- reference_results$calculations
  1098. zscore_interactions_reference <- reference_results$interactions
  1099. zscore_interactions_reference_joined <- reference_results$full_data
  1100. message("Calculating deletion strain(s) interactions scores")
  1101. df_deletion_stats <- calculate_summary_stats(
  1102. df = df_deletion,
  1103. variables = interaction_vars,
  1104. group_vars = c("OrfRep", "Gene", "conc_num", "conc_num_factor", "conc_num_factor_factor")
  1105. )$df_with_stats
  1106. deletion_results <- calculate_interaction_scores(df_deletion_stats, max_conc, bg_stats, group_vars = c("OrfRep", "Gene"))
  1107. zscore_calculations <- deletion_results$calculations
  1108. zscore_interactions <- deletion_results$interactions
  1109. zscore_interactions_joined <- deletion_results$full_data
  1110. # Writing Z-Scores to file
  1111. write.csv(zscore_calculations_reference, file = file.path(out_dir, "zscore_calculations_reference.csv"), row.names = FALSE)
  1112. write.csv(zscore_calculations, file = file.path(out_dir, "zscore_calculations.csv"), row.names = FALSE)
  1113. write.csv(zscore_interactions_reference, file = file.path(out_dir, "zscore_interactions_reference.csv"), row.names = FALSE)
  1114. write.csv(zscore_interactions, file = file.path(out_dir, "zscore_interactions.csv"), row.names = FALSE)
  1115. # Create interaction plots
  1116. message("Generating reference interaction plots")
  1117. reference_plot_configs <- generate_interaction_plot_configs(zscore_interactions_reference_joined, plot_type = "reference")
  1118. generate_and_save_plots(out_dir, "interaction_plots_reference", reference_plot_configs)
  1119. message("Generating deletion interaction plots")
  1120. deletion_plot_configs <- generate_interaction_plot_configs(zscore_interactions_joined, plot_type = "deletion")
  1121. generate_and_save_plots(out_dir, "interaction_plots", deletion_plot_configs)
  1122. # Define conditions for enhancers and suppressors
  1123. # TODO Add to study config?
  1124. threshold <- 2
  1125. enhancer_condition_L <- zscore_interactions$Avg_Zscore_L >= threshold
  1126. suppressor_condition_L <- zscore_interactions$Avg_Zscore_L <= -threshold
  1127. enhancer_condition_K <- zscore_interactions$Avg_Zscore_K >= threshold
  1128. suppressor_condition_K <- zscore_interactions$Avg_Zscore_K <= -threshold
  1129. # Subset data
  1130. enhancers_L <- zscore_interactions[enhancer_condition_L, ]
  1131. suppressors_L <- zscore_interactions[suppressor_condition_L, ]
  1132. enhancers_K <- zscore_interactions[enhancer_condition_K, ]
  1133. suppressors_K <- zscore_interactions[suppressor_condition_K, ]
  1134. # Save enhancers and suppressors
  1135. message("Writing enhancer/suppressor csv files")
  1136. write.csv(enhancers_L, file = file.path(out_dir, "zscore_interactions_deletion_enhancers_L.csv"), row.names = FALSE)
  1137. write.csv(suppressors_L, file = file.path(out_dir, "zscore_interactions_deletion_suppressors_L.csv"), row.names = FALSE)
  1138. write.csv(enhancers_K, file = file.path(out_dir, "zscore_interactions_deletion_enhancers_K.csv"), row.names = FALSE)
  1139. write.csv(suppressors_K, file = file.path(out_dir, "zscore_interactions_deletion_suppressors_K.csv"), row.names = FALSE)
  1140. # Combine conditions for enhancers and suppressors
  1141. enhancers_and_suppressors_L <- zscore_interactions[enhancer_condition_L | suppressor_condition_L, ]
  1142. enhancers_and_suppressors_K <- zscore_interactions[enhancer_condition_K | suppressor_condition_K, ]
  1143. # Save combined enhancers and suppressors
  1144. write.csv(enhancers_and_suppressors_L,
  1145. file = file.path(out_dir, "zscore_interactions_deletion_enhancers_and_suppressors_L.csv"), row.names = FALSE)
  1146. write.csv(enhancers_and_suppressors_K,
  1147. file = file.path(out_dir, "zscore_interaction_deletion_enhancers_and_suppressors_K.csv"), row.names = FALSE)
  1148. # Handle linear model based enhancers and suppressors
  1149. lm_threshold <- 2 # TODO add to study config?
  1150. enhancers_lm_L <- zscore_interactions[zscore_interactions$Z_lm_L >= lm_threshold, ]
  1151. suppressors_lm_L <- zscore_interactions[zscore_interactions$Z_lm_L <= -lm_threshold, ]
  1152. enhancers_lm_K <- zscore_interactions[zscore_interactions$Z_lm_K >= lm_threshold, ]
  1153. suppressors_lm_K <- zscore_interactions[zscore_interactions$Z_lm_K <= -lm_threshold, ]
  1154. # Save linear model based enhancers and suppressors
  1155. message("Writing linear model enhancer/suppressor csv files")
  1156. write.csv(enhancers_lm_L,
  1157. file = file.path(out_dir, "zscore_interactions_deletion_enhancers_lm_L.csv"), row.names = FALSE)
  1158. write.csv(suppressors_lm_L,
  1159. file = file.path(out_dir, "zscore_interactions_deletion_suppressors_lm_L.csv"), row.names = FALSE)
  1160. write.csv(enhancers_lm_K,
  1161. file = file.path(out_dir, "zscore_interactions_deletion_enhancers_lm_K.csv"), row.names = FALSE)
  1162. write.csv(suppressors_lm_K,
  1163. file = file.path(out_dir, "zscore_interactions_deletion_suppressors_lm_K.csv"), row.names = FALSE)
  1164. message("Generating rank plots")
  1165. rank_plot_configs <- generate_rank_plot_configs(
  1166. df = zscore_interactions_joined,
  1167. variables = interaction_vars,
  1168. is_lm = FALSE,
  1169. adjust = TRUE
  1170. )
  1171. generate_and_save_plots(out_dir = out_dir, filename = "rank_plots",
  1172. plot_configs = rank_plot_configs)
  1173. message("Generating ranked linear model plots")
  1174. rank_lm_plot_configs <- generate_rank_plot_configs(
  1175. df = zscore_interactions_joined,
  1176. variables = interaction_vars,
  1177. is_lm = TRUE,
  1178. adjust = TRUE
  1179. )
  1180. generate_and_save_plots(out_dir = out_dir, filename = "rank_plots_lm",
  1181. plot_configs = rank_lm_plot_configs)
  1182. message("Filtering and reranking plots")
  1183. interaction_threshold <- 2 # TODO add to study config?
  1184. # Formerly X_NArm
  1185. zscore_interactions_filtered <- zscore_interactions_joined %>%
  1186. filter(!is.na(Z_lm_L) & !is.na(Avg_Zscore_L)) %>%
  1187. mutate(
  1188. Overlap = case_when(
  1189. Z_lm_L >= interaction_threshold & Avg_Zscore_L >= interaction_threshold ~ "Deletion Enhancer Both",
  1190. Z_lm_L <= -interaction_threshold & Avg_Zscore_L <= -interaction_threshold ~ "Deletion Suppressor Both",
  1191. Z_lm_L >= interaction_threshold & Avg_Zscore_L <= interaction_threshold ~ "Deletion Enhancer lm only",
  1192. Z_lm_L <= interaction_threshold & Avg_Zscore_L >= interaction_threshold ~ "Deletion Enhancer Avg Zscore only",
  1193. Z_lm_L <= -interaction_threshold & Avg_Zscore_L >= -interaction_threshold ~ "Deletion Suppressor lm only",
  1194. Z_lm_L >= -interaction_threshold & Avg_Zscore_L <= -interaction_threshold ~ "Deletion Suppressor Avg Zscore only",
  1195. Z_lm_L >= interaction_threshold & Avg_Zscore_L <= -interaction_threshold ~ "Deletion Enhancer lm, Deletion Suppressor Avg Z score",
  1196. Z_lm_L <= -interaction_threshold & Avg_Zscore_L >= interaction_threshold ~ "Deletion Suppressor lm, Deletion Enhancer Avg Z score",
  1197. TRUE ~ "No Effect"
  1198. ),
  1199. lm_R_squared_L = summary(lm(Z_lm_L ~ Avg_Zscore_L))$r.squared,
  1200. lm_R_squared_K = summary(lm(Z_lm_K ~ Avg_Zscore_K))$r.squared,
  1201. lm_R_squared_r = summary(lm(Z_lm_r ~ Avg_Zscore_r))$r.squared,
  1202. lm_R_squared_AUC = summary(lm(Z_lm_AUC ~ Avg_Zscore_AUC))$r.squared
  1203. )
  1204. message("Generating filtered ranked plots")
  1205. rank_plot_filtered_configs <- generate_rank_plot_configs(
  1206. df = zscore_interactions_filtered,
  1207. variables = interaction_vars,
  1208. is_lm = FALSE,
  1209. adjust = FALSE,
  1210. overlap_color = TRUE
  1211. )
  1212. generate_and_save_plots(
  1213. out_dir = out_dir,
  1214. filename = "RankPlots_na_rm",
  1215. plot_configs = rank_plot_filtered_configs)
  1216. message("Generating filtered ranked linear model plots")
  1217. rank_plot_lm_filtered_configs <- generate_rank_plot_configs(
  1218. df = zscore_interactions_filtered,
  1219. variables = interaction_vars,
  1220. is_lm = TRUE,
  1221. adjust = FALSE,
  1222. overlap_color = TRUE
  1223. )
  1224. generate_and_save_plots(
  1225. out_dir = out_dir,
  1226. filename = "rank_plots_lm_na_rm",
  1227. plot_configs = rank_plot_lm_filtered_configs)
  1228. message("Generating correlation curve parameter pair plots")
  1229. correlation_plot_configs <- generate_correlation_plot_configs(zscore_interactions_filtered)
  1230. generate_and_save_plots(
  1231. out_dir = out_dir,
  1232. filename = "correlation_cpps",
  1233. plot_configs = correlation_plot_configs,
  1234. )
  1235. })
  1236. })
  1237. }
  1238. main()