Long-read proteogenomic atlas of human neuronal differentiation reveals isoform diversity informing neurodevelopmental risk mechanisms.

1. ### MAKING FIGURES ###
2. library(ggplot2)
3. library(ggalluvial)
4. library(dplyr)
5. library(tidyr)
6. library(data.table)
7. library(stringr)
8. library(readxl)
9. library(ggpubr)
10. library(cowplot)
11. library(ggrepel)
12. library(arrow)
13. library(gprofiler2)
17. setwd("")
20. # FIGURE 3A ---------------------------------------------------------------
21. sample_data <- data.frame(timepoint=c("t00", "t04", "t30",
22. "t00", "t04", "t30",
23. "t00", "t04", "t30",
24. "t00", "t04", "t30",
25. "t00", "t04", "t30",
26. "t00", "t04", "t30",
27. "t00", "t04", "t30",
28. "t00", "t04", "t30",
29. "t00", "t04", "t30"),
30. expr=c(2.5, 2.5, 2.5,
31. 2.5, 2.5, 0.5,
32. 2.5, 2.5, 4.5,
33. 4.5, 2.5, 2.5,
34. 4.5, 2.5, 0.5,
35. 4, 2, 4,
36. 0.5, 2.5, 2.5,
37. 1, 3, 1,
38. 0.5, 2.5, 4.5),
39. cluster = c("--", "--", "--",
40. "-D", "-D", "-D",
41. "-U", "-U", "-U",
42. "D-", "D-", "D-",
43. "DD", "DD", "DD",
44. "DU", "DU", "DU",
45. "U-", "U-", "U-",
46. "UD", "UD", "UD",
47. "UU", "UU", "UU"))
48. sample_data$timepoint <- as.factor(sample_data$timepoint)
49. sample_data$cluster <- as.factor(sample_data$cluster)
53. safe_colorblind_palette <- c("#88CCEE", "#CC6677", "#DDCC77", "#117733", "#332288", "#AA4499",
54. "#44AA99", "#999933", "#882255", "#661100", "#6699CC", "#888888")
57. Fig_3A <- ggplot(data=sample_data, aes(x=timepoint, y=expr, group = 1)) +
58. geom_line(aes(colour = cluster), linewidth = 1.5) +
59. geom_point(aes(colour = cluster), size = 2.75) +
60. scale_y_continuous(limits = c(0, 5), breaks = seq(0, 5, by = 2.5)) +
61. theme_classic(base_size = 8) +
62. theme(axis.text.y = element_blank(),
63. axis.text.x = element_text(size = 8),
64. axis.line = element_line(linewidth = 1),
65. axis.ticks = element_line(linewidth = 1),
66. legend.position = "none",
67. strip.background = element_blank(),
68. strip.text = element_text(face="bold", size=8)) +
69. scale_colour_manual(values = safe_colorblind_palette[1:9]) +
70. facet_wrap(sample_data$cluster, 3, 3, axes = "all") +
71. labs(x = "Time", y = "RNA or Protein\nExpression")
72. Fig_3A
76. # RIVER PLOT --------------------------------------------------------------
77. # GENES VS PROTEINS
78. #Compare genes and proteins. Riverplots each cluster vs each cluster.
79. safe_colorblind_palette <- c("#88CCEE", "#CC6677", "#DDCC77", "#117733", "#332288", "#AA4499",
80. "#44AA99", "#999933", "#882255", "#661100", "#6699CC", "#888888")
82. gene_clusters <- read.csv("./code/expression/output/gene_clusters_v2.csv")
83. protein_clusters <- read.csv("./code/expression/output/protein_clusters_v3.csv")
86. genesum_prot_clusters <- gene_clusters[,c(1,11)]
87. colnames(genesum_prot_clusters) <- c("gene_name", "mRNA")
89. genesum_prot_clusters <- merge(genesum_prot_clusters, protein_clusters[, c("GeneSymbol", "cluster")],
90. by.x = "gene_name", by.y = "GeneSymbol", all.x = TRUE)
92. genesum_prot_clusters <- na.omit(genesum_prot_clusters) #8493 remaining out of 8878
93. colnames(genesum_prot_clusters)[3] <- "Protein"
95. genesum_prot_long <- genesum_prot_clusters %>% pivot_longer(cols = c(mRNA, Protein),
96. names_to = "gene_prot",
97. values_to = "cluster")
100. genesum_prot_long$cluster <- as.factor(genesum_prot_long$cluster)
101. genesum_prot_long$rep <- rep(1:nrow(genesum_prot_clusters), each=2)
103. label_fix <- data.frame(
104. gene_prot = c("Protein", "Protein"),
105. y = c(1131.5, 1012.5),
106. label = c("DD", "DU"),
107. x = c(1.95, 2.05))
109. genesum_prot_long$label <- as.character(genesum_prot_long$cluster)
110. genesum_prot_long$label[genesum_prot_long$cluster %in% c("DD", "DU") &
111. genesum_prot_long$gene_prot == "Protein"] <- ""
114. genesum_river <- ggplot(genesum_prot_long,
115. aes(x = gene_prot, stratum = cluster, alluvium = rep,
116. fill = cluster, label = label)) +
117. geom_alluvium(aes(fill = cluster),
118. curve_type = "sigmoid") +
119. geom_stratum() +
120. geom_text(stat = "stratum",
121. aes(label = label),
122. fontface = "bold",
123. size = 3.75) +
124. scale_fill_manual(values = safe_colorblind_palette[1:9]) +
125. scale_y_continuous(expand = c(0, 0), limits = c(0, 8600)) +
126. scale_x_discrete(expand = c(0, 0)) + #move y axis closer
127. theme_void(base_size = 8) +
128. theme(legend.position = "none",
129. axis.title.x = element_blank(),
130. axis.text.y = element_text(margin = margin(r = 5)),
131. axis.text = element_text(size = 8),
132. axis.ticks.y = element_line(size = 1),
133. axis.ticks.length.y = unit(3, "pt"),
134. axis.line.y = element_line(size = 0.75, colour = "black"),
135. plot.margin = margin(
136. t = 6,
137. r = 6,
138. b = 6,
139. l = 6,
140. unit = "pt" ))
142. genesum_river_updated <- genesum_river +
143. geom_text(data = label_fix,
144. aes(x = x, y = y, label = label),
145. inherit.aes = FALSE,
146. fontface = "bold", size = 3.75)
147. genesum_river_updated
150. #
151. table(genesum_prot_clusters$mRNA) #59% changing
152. table(genesum_prot_clusters$Protein) #42% changing
156. # EXAMPLE GENES -----------------------------------------------------------
157. ###Assigning concordant (0) and discordant 1 or 2
158. genesum_prot_clusters$mRNA <- as.character(genesum_prot_clusters$mRNA)
159. genesum_prot_clusters$Protein <- as.character(genesum_prot_clusters$Protein)
161. genesum_prot_clusters[,c("Genes_1", "Genes_2")] <- str_split_fixed(genesum_prot_clusters$mRNA, "", 2)
162. genesum_prot_clusters[,c("Proteins_1", "Proteins_2")] <- str_split_fixed(genesum_prot_clusters$Protein, "", 2)
164. genesum_prot_clusters$match1 <- ifelse((genesum_prot_clusters$Genes_1 == genesum_prot_clusters$Proteins_1), 0, 1)
165. genesum_prot_clusters$match2 <- ifelse((genesum_prot_clusters$Genes_2 == genesum_prot_clusters$Proteins_2), 0, 1)
166. genesum_prot_clusters$Level <- genesum_prot_clusters$match1 + genesum_prot_clusters$match2
168. sum(genesum_prot_clusters$Level == 0)/nrow(genesum_prot_clusters) #45.5
169. sum(genesum_prot_clusters$Level == 1)/nrow(genesum_prot_clusters) #42.5
170. sum(genesum_prot_clusters$Level == 2)/nrow(genesum_prot_clusters) #12.0
173. #Intersect the level2 with SFARI genes
174. SFARI <- read.csv("./data/SFARI-Gene_genes_04-03-2025release_04-15-2025export.csv")
176. level_SFARI <- genesum_prot_clusters[genesum_prot_clusters$gene_name %in% SFARI$gene.symbol, ]
177. level_SFARI <- level_SFARI %>% filter((mRNA != "--" | Protein != "--"))
180. #Make example plots for gene to protein fig.
181. #Prepare gene data.
182. gene_cpm <- read.csv("./code/expression/output/gene_CPM.csv")
184. gene_SFARI <- gene_cpm[gene_cpm$gene_name %in% level_SFARI$gene_name,]
185. gene_SFARI <- merge(gene_SFARI, gene_clusters[, c("gene_name", "cluster")],
186. by.x = "gene_name", by.y = "gene_name", all.x = TRUE)
188. #Pivot longer.
189. gene_summ <- gene_SFARI %>% pivot_longer(
190. cols = starts_with("sum_"), # Select columns that start with "sum_"
191. names_to = c("Type", "Group"), # Create two new columns: "Type" (iPSC, NPC, CN) and "Group"
192. names_sep = "_", # Split the names at "_"
193. values_to = "Value")
195. #Log2 transform the data, then calculate mean and 95% CI, and SD
196. gene_long <- gene_summ %>%
197. mutate(Value_log2 = log2(Value + 1))
199. gene_long$Group <- factor(gene_long$Group, levels = c("iPSC", "NPC", "CN"),
200. labels = c("t00", "t04", "t30"))
202. # Calculate mean and 95% CI, and SD
203. gene_summary <- gene_long %>%
204. group_by(gene_name, Group) %>%
205. summarise(
206. Mean = mean(Value_log2),
207. SE = sd(Value_log2) / sqrt(n()),
208. CI_Lower = Mean - qt(0.975, df = n() - 1) * SE,
209. CI_Upper = Mean + qt(0.975, df = n() - 1) * SE,
210. sd = sd(Value_log2, na.rm = TRUE) )
211. gene_summary <- merge(gene_summary, gene_clusters[, c("gene_name", "cluster")],
212. by.x = "gene_name", by.y = "gene_name", all.x = F)
214. ###Prepare protein data
215. prot_quant <- read.csv("./code/expression/output/MSstatsTMT_abundance.csv")
216. final_filtered_diff_prot <- read.csv("./code/expression/output/MSstatsTMT_diff_prot_expr_filtered.csv")
218. #Calculate median across fractions:
219. prot_quant_median <- prot_quant %>% group_by(Protein, Channel) %>%
220. summarise(Median_Abundance = median(Abundance)) %>%
221. left_join(prot_quant %>% dplyr::select(Channel, Condition) %>% distinct(Channel, Condition) , by = "Channel")
223. summary_prot_abund <- prot_quant_median %>%
224. group_by(Protein, Condition) %>%
225. summarise(
226. mean = mean(Median_Abundance, na.rm = TRUE),
227. sd = sd(Median_Abundance, na.rm = TRUE),
228. n = n()
229. )
231. summary_prot_abund <- summary_prot_abund %>%
232. left_join(final_filtered_diff_prot %>% dplyr::select(Protein, GeneSymbol) %>% distinct(Protein, GeneSymbol), by = "Protein")
233. summary_prot_abund <- summary_prot_abund %>% left_join(genesum_prot_clusters %>% dplyr::select(gene_name, Protein), by = c("GeneSymbol" = "gene_name"))
235. summary_prot_abund$Condition <- factor(summary_prot_abund$Condition, levels = c("t00", "t4", "t30"),
236. labels = c("t00", "t04", "t30"))
239. #Chosen examples
240. gene_ex <- gene_summary[gene_summary$gene_name %in% level_SFARI$gene_name,]
241. prot_ex <- summary_prot_abund[summary_prot_abund$GeneSymbol %in% level_SFARI$gene_name,]
243. #Make 6 indiv plots per gene b/c I need gene name on the left, and specific scales
244. THBS1_gene <- gene_ex %>% filter(gene_name == "THBS1") %>% ggplot(aes(x = Group, y = Mean, group = gene_name, color = cluster)) +
245. geom_line(linewidth = 1.5) +
246. geom_point(size = 2.75) +
247. geom_errorbar(aes(ymin = Mean - sd, ymax = Mean + sd), width = 0.2, linewidth = 1) +
248. scale_y_continuous(limits = c(0, 4.2)) +
249. theme_classic(base_size = 8) +
250. labs(title = "mRNA", subtitle = "log2(CPM + 1)", y = "THBS1") +
251. scale_color_manual(values=safe_colorblind_palette[4],
252. breaks = c("D-")) +
253. theme(axis.line = element_line(linewidth = 1),
254. axis.ticks = element_line(linewidth = 1, colour = "black"),
255. plot.title = element_text(hjust = 0.5, size = 8, face = "bold"),
256. plot.subtitle = element_text(hjust = 0.5, size = 6.4),
257. legend.position = "none",
258. axis.text = element_text(size = 8),
259. axis.title.x = element_blank(),
260. axis.title.y = element_text(face = "bold")) +
261. geom_text(aes(x=Inf,y=Inf,
262. hjust=c(1.5),
263. vjust=c(1.5),
264. label="D-",
265. fontface = "bold"),
266. size = 4)
267. THBS1_gene
269. MYO1E_gene <- gene_ex %>% filter(gene_name == "MYO1E") %>% ggplot(aes(x = Group, y = Mean, group = gene_name, color = cluster)) +
270. geom_line(linewidth = 1.5) +
271. geom_point(size = 2.75) +
272. geom_errorbar(aes(ymin = Mean - sd, ymax = Mean + sd), width = 0.2, linewidth = 1) +
273. scale_y_continuous(limits = c(1, 8)) +
274. theme_classic(base_size = 8) +
275. labs(y = "MYO1E") +
276. scale_color_manual(values=safe_colorblind_palette[4],
277. breaks = c("D-")) +
278. theme(axis.line = element_line(linewidth = 1),
279. axis.ticks = element_line(linewidth = 1, colour = "black"),
280. legend.position = "none",
281. axis.text = element_text(size = 8),
282. axis.title.x = element_blank(),
283. axis.title.y = element_text(face = "bold")) +
284. geom_text(aes(x=Inf,y=Inf,hjust=c(1.5),
285. vjust=c(1.5),label="D-", fontface = "bold"),
286. size = 4)
287. MYO1E_gene
289. DPP4_gene <- gene_ex %>% filter(gene_name == "DPP4") %>% ggplot(aes(x = Group, y = Mean, group = gene_name, color = cluster)) +
290. geom_line(linewidth = 1.5) +
291. geom_point(size = 2.75) +
292. geom_errorbar(aes(ymin = Mean - sd, ymax = Mean + sd), width = 0.2, linewidth = 1) +
293. scale_y_continuous(limits = c(-0.1, 4)) +
294. theme_classic(base_size = 8) +
295. labs(y = "DPP4") +
296. scale_color_manual(values=safe_colorblind_palette[4],
297. breaks = c("D-")) +
298. theme(axis.line = element_line(linewidth = 1),
299. axis.ticks = element_line(linewidth = 1, colour = "black"),
300. legend.position = "none",
301. axis.text = element_text(size = 8),
302. axis.title.x = element_blank(),
303. axis.title.y = element_text(face = "bold")) +
304. geom_text(aes(x=Inf,y=Inf,hjust=c(1.5),
305. vjust=c(1.5),label="D-", fontface = "bold"),
306. size = 4)
307. DPP4_gene
309. #Proteins
310. THBS1_prot <- prot_ex %>% filter(GeneSymbol == "THBS1") %>% ggplot(aes(x = Condition, y = mean, group = Protein.y, color = Protein.y)) +
311. geom_line(linewidth = 1.5) +
312. geom_point(size = 2.75) +
313. geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd), width = 0.2, linewidth = 1) +
314. scale_y_continuous(limits = c(3.8, 7.5)) +
315. theme_classic(base_size = 8) +
316. labs(title = "Protein", subtitle = "log2(Abundance)") +
317. scale_color_manual(values=safe_colorblind_palette[c(4)],
318. breaks = c("D-")) +
319. theme(axis.line = element_line(linewidth = 1),
320. axis.ticks = element_line(linewidth = 1, colour = "black"),
321. plot.title = element_text(hjust = 0.5, size = 8, face = "bold"),
322. plot.subtitle = element_text(hjust = 0.5, size = 6.4),
323. legend.position = "none",
324. axis.text = element_text(size = 8),
325. axis.title = element_blank()) +
326. geom_text(aes(x=Inf,y=Inf,hjust=c(1.5),
327. vjust=c(1.5),label="D-", fontface = "bold"),
328. size = 4)
329. THBS1_prot
331. MYO1E_prot <- prot_ex %>% filter(GeneSymbol == "MYO1E") %>% ggplot(aes(x = Condition, y = mean, group = Protein.y, color = Protein.y)) +
332. geom_line(linewidth = 1.5) +
333. geom_point(size = 2.75) +
334. geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd), width = 0.2, linewidth = 1) +
335. scale_y_continuous(limits = c(3.8, 7.5)) +
336. theme_classic(base_size = 8) +
337. scale_color_manual(values=safe_colorblind_palette[5],
338. breaks = c("DD")) +
339. theme(axis.line = element_line(linewidth = 1),
340. axis.ticks = element_line(linewidth = 1, colour = "black"),
341. legend.position = "none",
342. axis.text = element_text(size = 8),
343. axis.title = element_blank()) +
344. geom_text(aes(x=Inf,y=Inf,hjust=c(1.5),
345. vjust=c(1.5),label="DD", fontface = "bold"),
346. size = 4)
347. MYO1E_prot
349. DPP4_prot <- prot_ex %>% filter(GeneSymbol == "DPP4") %>% ggplot(aes(x = Condition, y = mean, group = Protein.y, color = Protein.y)) +
350. geom_line(linewidth = 1.5) +
351. geom_point(size = 2.75) +
352. geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd), width = 0.2, linewidth = 1) +
353. scale_y_continuous(limits = c(3.8, 7.5)) +
354. theme_classic(base_size = 8) +
355. scale_color_manual(values=safe_colorblind_palette[c(4,1,8)],
356. breaks = c("D-", "--", "UD")) +
357. theme(axis.line = element_line(linewidth = 1),
358. axis.ticks = element_line(linewidth = 1, colour = "black"),
359. legend.position = "none",
360. axis.text = element_text(size = 8),
361. axis.title = element_blank()) +
362. geom_text(aes(x=Inf,y=Inf,hjust=c(1.5),
363. vjust=c(1.5),label="UD", fontface = "bold"),
364. size = 4)
365. DPP4_prot
368. combined_plot <- ggarrange(THBS1_gene, THBS1_prot,
369. MYO1E_gene, MYO1E_prot,
370. DPP4_gene, DPP4_prot ,
371. ncol = 2,
372. nrow = 3,
373. heights = c(1.32, 1, 1))
374. combined_plot
377. # 9 CLUSTER ASD RISK GENE ENRICHMENT------------------------------------------
378. ####GO term enrichment for 9 protein clusters:
379. #bg proteins is all genes with mRNA and protein expression
380. bg.prot <- genesum_prot_clusters %>% pull(gene_name) %>% unique()
383. #use 0 to replace -
384. gene_lists <- list(genes_00 = genesum_prot_clusters %>% filter(Protein == "--") %>% pull(gene_name) %>% unique(),
385. genes_0D = genesum_prot_clusters %>% filter(Protein == "-D") %>% pull(gene_name) %>% unique(),
386. genes_U0 = genesum_prot_clusters %>% filter(Protein == "U-") %>% pull(gene_name) %>% unique(),
387. genes_0U = genesum_prot_clusters %>% filter(Protein == "-U") %>% pull(gene_name) %>% unique(),
388. genes_D0 = genesum_prot_clusters %>% filter(Protein == "D-") %>% pull(gene_name) %>% unique(),
389. genes_UU = genesum_prot_clusters %>% filter(Protein == "UU") %>% pull(gene_name) %>% unique(),
390. genes_UD = genesum_prot_clusters %>% filter(Protein == "UD") %>% pull(gene_name) %>% unique(),
391. genes_DU = genesum_prot_clusters %>% filter(Protein == "DU") %>% pull(gene_name) %>% unique(),
392. genes_DD = genesum_prot_clusters %>% filter(Protein == "DD") %>% pull(gene_name) %>% unique())
394. ###Run GO term search:
395. run_GO_enrichment <- function(gene_list, bg, searches) {
396. go_result <- gost(gene_list,
397. organism = "hsapiens",
398. significant = T,
399. exclude_iea = F,
400. user_threshold = 0.05,
401. correction_method = "fdr",
402. custom_bg = bg,
403. sources = searches)
404. return(go_result)
405. }
406. searches <- c("GO:BP", "GO:MF", "GO:CC", "MIRNA")
407. # Initialize an empty list to store results
408. enrichment_results <- list()
409. i = 1
410. # Loop through each gene list and perform GO enrichment
411. for (i in seq_along(gene_lists)) {
412. enrichment_results[[names(gene_lists)[i]]] <- run_GO_enrichment(gene_lists[[i]],
413. bg = bg.prot,
414. searches = searches)
415. }
417. enrichment_results$genes_00$result$query <- "--"
418. enrichment_results$genes_0D$result$query <- "-D"
419. enrichment_results$genes_U0$result$query <- "U-"
420. enrichment_results$genes_0U$result$query <- "-U"
421. enrichment_results$genes_D0$result$query <- "D-"
422. enrichment_results$genes_UU$result$query <- "UU"
423. enrichment_results$genes_UD$result$query <- "UD"
424. enrichment_results$genes_DU$result$query <- "DU"
425. enrichment_results$genes_DD$result$query <- "DD"
428. GO_list <- rbind(enrichment_results$genes_00$result,
429. enrichment_results$genes_0D$result,
430. enrichment_results$genes_U0$result,
431. enrichment_results$genes_0U$result,
432. enrichment_results$genes_D0$result,
433. enrichment_results$genes_UU$result,
434. enrichment_results$genes_UD$result,
435. enrichment_results$genes_DU$result,
436. enrichment_results$genes_DD$result)
438. GO_list$query <- as.character(GO_list$query)
439. write.csv(GO_list[,-14], "./tables/protein_clusters_GO_terms.csv",
440. row.names = F)
443. #########SFARI Gene Enrichment
444. #SFARI genes subset to only proteins expressed:
445. SFARI_prot <- SFARI %>% filter(gene.symbol %in% genesum_prot_clusters$gene_name) #805
447. # Function to compute Fisher's exact test for pairwise comparisons
448. compute_fisher_enrich <- function(setA, setB, bg) {
449. # Calculate overlap and unique genes
450. overlap <- length(intersect(setA, setB))
451. inAnotB <- length(setdiff(setA, setB))
452. inBnotA <- length(setdiff(setB, setA))
453. inNeither <- length(bg) - length(union(setA, setB))
455. # Create the contingency table
456. contingency_table <- matrix(c(overlap, inAnotB, inBnotA, inNeither), nrow = 2)
458. # Perform Fisher's Exact Test
459. fisher_result <- fisher.test(contingency_table, alternative = "greater")
461. return(c(overlap = overlap, inAnotB = inAnotB, inBnotA = inBnotA, inNeither = inNeither, p_value = fisher_result$p.value))
462. }
465. #Testing for enrichment in the each cluster:
466. results_list <- list()
468. i = 1
469. for (i in 1:9) {
470. setA <- SFARI_prot$gene.symbol
471. setB <- gene_lists[[i]]
472. result <- compute_fisher_enrich(setA, setB, bg.prot)
474. # Store the results in the list, and create the pair name
475. pair_name <- paste0("SFARI vs ", names(gene_lists[i]))
476. results_list[[pair_name]] <- result
477. }
479. results_df_SFARI_prot <- do.call(rbind, results_list)
480. results_df_SFARI_prot <- data.frame(pair = names(results_list), results_df_SFARI_prot)
481. results_df_SFARI_prot$p_adjusted <- p.adjust(results_df_SFARI_prot$p_value, method = "bonferroni")
482. results_df_SFARI_prot$gene <- "Protein"
484. SFARI_gene_enrich <- results_df_SFARI_prot
486. SFARI_gene_enrich$cluster <- gsub("0", "-", substr(SFARI_gene_enrich$pair, nchar(SFARI_gene_enrich$pair)-1, nchar(SFARI_gene_enrich$pair)))
487. SFARI_gene_enrich$cluser <- as.factor(SFARI_gene_enrich$cluster)
488. Fig3d_title <- "SFARI Gene Enrichment"
490. SFARI_gene_enrich$gene <- as.factor(SFARI_gene_enrich$gene)
493. Fig_3d <- ggplot(SFARI_gene_enrich, aes(x = cluster, y = -log10(p_adjusted)) ) +
494. geom_col(aes(x = cluster, y = -log10(p_adjusted), fill = cluster),
495. colour = "black",
496. position = position_dodge()) +
497. labs(y = expression(paste(-log[10], "(adj. ", italic("P"), ")")), title = Fig3d_title) +
498. coord_cartesian(ylim = c(0, 15.5)) +
499. scale_y_continuous(expand = c(0, 0)) +
500. geom_hline(yintercept=-log10(0.05), linetype='dashed', col = "black", linewidth = 0.7) +
501. theme_classic(base_size = 8) +
502. scale_fill_manual(values=safe_colorblind_palette) +
503. theme(axis.title.x = element_blank(),
504. axis.title.y = element_text(size = 8),
505. panel.grid.major.x = element_blank(),
506. plot.title = element_text(hjust = 0.5, size = 8, face = "bold"),
507. legend.position = "none" )
508. Fig_3d
511. #Export figure
512. Figure_3 <- ggdraw() +
513. draw_plot(Fig_3A, x = 0, y = 0.63, width = 0.5, height = 0.375) +
514. draw_plot(genesum_river_updated, x = 0.5, y = 0, width = 0.5, height = 1) +
515. draw_plot(combined_plot, x = 0, y = 0.19, width = 0.5, height = 0.44) +
516. draw_plot(Fig_3d, x = 0.0, y = 0, width = 0.5, height = 0.191) +
517. draw_plot_label(label = c("A", "B", "C", "D"),
518. size = 14,
519. x = c(0, 0.5, 0, 0),
520. y = c(1, 1, 0.63, 0.1911))
521. #Figure_3
523. pdf(paste0("./figures/Figure_3.pdf"), height = 6.5, width = 6.5)
524. print(Figure_3)
525. dev.off()
528. svg(paste0("./figures/Figure_3.svg"), height = 6.5, width = 6.5)
529. print(Figure_3)
530. dev.off()
533. # FIGURE 4 --------------------------------------------------------------
534. classification <- read_parquet("./data/final_classification.parquet") #182371
536. ###IsoformSwitches
537. aSwitchList_part2 <- readRDS("./code/IsoformSwitchAnalyzeR/output/isoformswitch_part2_v4.rds")
538. aSwitchList_part2[["isoformFeatures"]]$gene_name <- aSwitchList_part2[["isoformFeatures"]]$gene_id
543. # FIGURE 4A ---------------------------------------------------------------
544. ###DTU PROPORTIONS (bars)
545. DTU <- read.csv("./code/IsoformSwitchAnalyzeR/output/DTU_table.csv")
546. DTU_filter <- DTU %>% filter(DTU_qval < 0.05 &
547. (DTU_dIF < -0.05 | DTU_dIF > 0.05))
549. DTU_filter <- DTU_filter[,-c(1:3)]
550. DTU_filter$condition <- paste0(DTU_filter$condition_1, " vs ", DTU_filter$condition_2)
551. DTU_filter <- merge(DTU_filter, classification[, c("isoform", "structural_category")],
552. by.x = "isoform_id", by.y = "isoform", all.x = T)
554. rep_str = c('full-splice_match'='FSM','incomplete-splice_match' = 'ISM', 'novel_in_catalog'='NIC','novel_not_in_catalog'='NNC')
555. DTU_filter$structural_category <- str_replace_all(DTU_filter$structural_category, rep_str)
557. #All conditions
558. DTU_summ <- DTU_filter %>%
559. group_by(condition, structural_category) %>%
560. summarise(count = n()) %>%
561. group_by(condition) %>%
562. mutate(perc = (count/sum(count))*100) %>%
563. ungroup()
565. DTU_summ$percentage = paste0(round(DTU_summ$perc, 0), "%")
567. ####DTE
568. res_t30_vs_t00 <- read.csv("./code/expression/output/DESeq2_tr_t30_vs_t00.csv")
569. res_t30_vs_t00$condition <- "t00 vs t30"
571. res_t04_vs_t00 <- read.csv("./code/expression/output/DESeq2_tr_t04_vs_t00.csv")
572. res_t04_vs_t00$condition <- "t00 vs t04"
575. res_t30_vs_t04 <- read.csv("./code/expression/output/DESeq2_tr_t30_vs_t04.csv")
576. res_t30_vs_t04$condition <- "t04 vs t30"
579. DTE_prefilter <- rbind(res_t30_vs_t00, res_t04_vs_t00)
580. DTE_prefilter <- rbind(DTE_prefilter, res_t30_vs_t04)
583. DTE_prefilter <- merge(DTE_prefilter, classification[, c("isoform", "associated_gene", "structural_category")],
584. by.x = "pb_id", by.y = "isoform", all.x = T)
586. #Just take the transcripts that meet the DTU pre-filter:
587. prefilter <- aSwitchList_part2$isoformRepIF$isoform_id
589. DTE_prefilter <- DTE_prefilter %>% filter(pb_id %in% prefilter) #68,801
590. DTE <- DTE_prefilter %>% filter(padj < 0.05 &
591. (log2FoldChange < -1 | log2FoldChange > 1))
594. DTE$structural_category <- str_replace_all(DTE$structural_category, rep_str)
596. DTE_summ <- DTE %>%
597. group_by(condition, structural_category) %>%
598. summarise(count = n()) %>%
599. group_by(condition) %>%
600. mutate(perc = (count/sum(count))*100) %>%
601. ungroup()
603. DTE_summ$percentage = paste0(round(DTE_summ$perc, 0), "%")
606. DTE_summ$type <- "DTE"
607. DTU_summ$type <- "DTU"
608. DTU_summ <- rbind(DTU_summ, DTE_summ)
610. #FSM NIC NNC ISM
611. colours <- c("#009E73", "#D55E00", "#E69F00", "#0072B2")
615. DTU_summ$condition <- factor(DTU_summ$condition, levels = c("t00 vs t04", "t04 vs t30", "t00 vs t30"))
618. bars_by_tr <- ggplot(DTU_summ, aes(x = condition, y = count/(1000), fill = structural_category)) +
619. geom_bar(stat = "identity") +
620. geom_text(aes(label = percentage), position = position_stack(vjust = 0.5), size = 1.9) +
621. theme_classic(base_size = 8) +
622. scale_fill_manual(values = colours[c(1,4,2,3)]) +
623. theme(
624. legend.title = element_blank(),
625. legend.key.size = unit(0.35, 'cm'),
626. axis.line = element_line(linewidth = 0.6),
627. axis.ticks = element_line(linewidth = 0.6, colour = "black"),
628. axis.title.x = element_blank(),
629. axis.title.y = element_text(size = 8, colour = "black"),
630. axis.text = element_text(size = 6.4, colour = "black"),
631. strip.background = element_blank(),
632. strip.text = element_text(size = 8, face = "bold"),
633. axis.text.x=element_text(angle=-45, hjust = 0, vjust=1)) +
634. labs(y = expression(paste("Isoforms (x ", 10^{3}, ")"))) +
635. guides(y = guide_axis(minor.ticks = TRUE)) +
636. facet_wrap(~factor(type), scales = "free")
637. bars_by_tr
638. #
641. # FUNC. CONSEQ ENRICHMENT------------------------------------------------------
642. #Need to run the function fully to extract genes, etc.
643. switchAnalyzeRlist = aSwitchList_part2
644. consequencesToAnalyze = 'all'
645. alpha=0.05
646. dIFcutoff = 0.05
647. countGenes = TRUE
648. analysisOppositeConsequence=FALSE
649. plot=TRUE
650. localTheme = theme_bw(base_size = 12)
651. minEventsForPlotting = 10
652. returnResult=TRUE
653. returnSummary=TRUE
655. ### Consequences to analyze
656. acceptedTypes <- c(
657. # Transcript
658. 'tss',
659. 'tts',
660. 'last_exon',
661. 'isoform_length',
662. 'exon_number',
663. 'intron_structure',
664. 'intron_retention',
665. 'isoform_class_code',
666. # cpat
667. 'coding_potential',
668. # ORF
669. 'ORF_genomic',
670. 'ORF_length',
671. '5_utr_length',
672. '3_utr_length',
673. # seq similarity
674. 'isoform_seq_similarity',
675. 'ORF_seq_similarity',
676. '5_utr_seq_similarity',
677. '3_utr_seq_similarity',
678. # ORF
679. 'NMD_status',
680. # pfam
681. 'domains_identified',
682. 'genomic_domain_position',
683. 'domain_length',
684. 'domain_isotype',
686. # SignalIP
687. 'signal_peptide_identified',
689. # IDR
690. 'IDR_identified',
691. 'IDR_length',
692. 'IDR_type',
694. # sub cell
695. 'sub_cell_location',
696. 'sub_cell_shift_to_cell_membrane',
697. 'sub_cell_shift_to_cytoplasm',
698. 'sub_cell_shift_to_nucleus',
699. 'sub_cell_shift_to_Extracellular'
701. # topology
702. # 'isoform_topology'
703. )
705. consequencesAnalyzed <- unique(switchAnalyzeRlist$switchConsequence$featureCompared)
706. if ('all' %in% consequencesToAnalyze) {
707. consequencesToAnalyze <- consequencesAnalyzed
708. }
710. consequencesNotAnalyzed <- setdiff(consequencesToAnalyze, consequencesAnalyzed)
711. if (length(consequencesNotAnalyzed)) {
712. warning(
713. paste(
714. 'The following consequences appear not to have been analyzed and will therefor not be summarized:',
715. paste(consequencesNotAnalyzed, collapse = ', '),
716. sep = ' '
717. )
718. )
719. }
721. ### Extract non-location consequences
722. if(TRUE) {
723. localConseq <- switchAnalyzeRlist$switchConsequence[
724. which( !is.na(
725. switchAnalyzeRlist$switchConsequence$switchConsequence
726. ))
727. ,]
728. localConseq <- localConseq[which(
729. ! grepl('switch', localConseq$switchConsequence)
730. ),]
732. ### make list with levels
733. levelList <- list(
734. tss=c('Tss more upstream','Tss more downstream'),
735. tts=c('Tts more downstream','Tts more upstream'),
736. last_exon=c('Last exon more downstream','Last exon more upstream'),
737. isoform_length=c('Length gain','Length loss'),
738. isoform_seq_similarity=c('Length gain','Length loss'),
739. exon_number=c('Exon gain','Exon loss'),
740. intron_retention=c('Intron retention gain','Intron retention loss'),
741. ORF_length=c('ORF is longer','ORF is shorter'),
742. ORF=c('Complete ORF loss','Complete ORF gain'),
743. x5_utr_length=c('5UTR is longer','5UTR is shorter'),
744. x3_utr_length=c('3UTR is longer','3UTR is shorter'),
745. NMD_status=c('NMD sensitive','NMD insensitive'),
746. coding_potential=c('Transcript is coding','Transcript is Noncoding'),
747. domains_identified=c('Domain gain','Domain loss'),
748. domain_length=c('Domain length gain','Domain length loss'),
749. domain_isotype=c('Domain non-reference isotype gain','Domain non-reference isotype loss'),
750. IDR_identified = c('IDR gain','IDR loss'),
751. IDR_length = c('IDR length gain', 'IDR length loss'),
752. IDR_type = c('IDR w binding region gain', 'IDR w binding region loss'),
753. signal_peptide_identified=c('Signal peptide gain','Signal peptide loss'),
754. sub_cell_location = c('SubCell location gain','SubCell location loss'),
755. sub_cell_shift_to_cell_membrane = c('SubCell location memb gain','SubCell location memb loss'),
756. sub_cell_shift_to_cytoplasm = c('SubCell location cyto gain','SubCell location cyto loss'),
757. sub_cell_shift_to_nucleus = c('SubCell location nucl gain','SubCell location nucl loss'),
758. sub_cell_shift_to_Extracellular = c('SubCell location ext cell gain','SubCell location ext cell loss'),
759. isoform_topology = c('Topology complexity gain','Topology complexity loss'),
760. extracellular_region_count = c('Extracellular region gain', 'Extracellular region loss'),
761. intracellular_region_count = c('Intracellular region gain', 'Intracellular region loss'),
762. extracellular_region_length = c('Extracellular length gain','Extracellular length loss'),
763. intracellular_region_length = c('Intracellular length gain','Intracellular length loss')
764. )
765. levelListDf <- plyr::ldply(levelList, function(x) data.frame(feature=x, stringsAsFactors = FALSE))
767. ### Add consequence pairs
768. localConseq$conseqPair <- levelListDf$.id[match(localConseq$switchConsequence, levelListDf$feature)]
769. localConseq <- localConseq[which( !is.na(localConseq$conseqPair)),]
771. ### Subset to consequences analyzed
772. localConseq <- localConseq[which(
773. localConseq$featureCompared %in% consequencesToAnalyze
774. ),]
775. }
777. ### Subset to significant features
778. if(TRUE) {
779. ### Extract Sig iso
780. isoResTest <-
781. any(!is.na(
782. switchAnalyzeRlist$isoformFeatures$isoform_switch_q_value
783. ))
784. if (isoResTest) {
785. sigIso <- switchAnalyzeRlist$isoformFeatures[which(
786. switchAnalyzeRlist$isoformFeatures$isoform_switch_q_value < alpha &
787. abs(switchAnalyzeRlist$isoformFeatures$dIF) > dIFcutoff
788. ),
789. c('iso_ref', 'gene_ref')]
790. } else {
791. sigIso <- switchAnalyzeRlist$isoformFeatures[which(
792. switchAnalyzeRlist$isoformFeatures$gene_switch_q_value < alpha &
793. abs(switchAnalyzeRlist$isoformFeatures$dIF) > dIFcutoff
794. ),
795. c('iso_ref', 'gene_ref')]
796. }
798. if(isoResTest) {
799. localConseq <- localConseq[which(
800. localConseq$iso_ref_down %in% sigIso$iso_ref |
801. localConseq$iso_ref_up %in% sigIso$iso_ref
802. ),]
803. } else {
804. localConseq <- localConseq[which(
805. localConseq$gene_ref %in% sigIso$gene_ref
806. ),]
807. }
808. }
811. ### Summarize gain vs loss for each consequence in each condition
812. consequenceBalance <- plyr::ddply(
813. .data = localConseq,
814. .variables = c('condition_1','condition_2','conseqPair'),
815. #.inform = TRUE,
816. .fun = function(aDF) { # aDF <- localConseq[1:20,]
817. ### Add levels
818. if(analysisOppositeConsequence) {
819. localLvl <- rev(sort(
820. levelList[[ aDF$conseqPair[1] ]]
821. ))
822. } else {
823. localLvl <- sort(
824. levelList[[ aDF$conseqPair[1] ]]
825. )
826. }
827. aDF$switchConsequence <- factor(
828. aDF$switchConsequence,
829. levels=localLvl
830. )
832. ### Summarize category
833. if( countGenes ) {
834. df2 <- aDF[
835. which(!is.na(aDF$switchConsequence)),
836. c('gene_id','switchConsequence')
837. ]
838. localNumber <- plyr::ddply(df2, .drop = FALSE, .variables = 'switchConsequence', function(x) {
839. data.frame(
840. Freq = length(unique(x$gene_id))
841. )
842. })
843. colnames(localNumber)[1] <- 'Var1'
844. } else {
845. localNumber <- as.data.frame(table(aDF$switchConsequence))
846. }
848. if(nrow(localNumber) == 2) {
849. localTest <- suppressWarnings(
850. stats::binom.test(localNumber$Freq[1], sum(localNumber$Freq))
851. )
853. localRes <- data.frame(
854. feature=stringr::str_c(
855. localNumber$Var1[1],
856. ' (paired with ',
857. localNumber$Var1[2],
858. ')'
859. ),
860. propOfRelevantEvents=localTest$estimate,
861. stringsAsFactors = FALSE
862. )
864. localRes$propCiLo <- min(localTest$conf.int)
865. localRes$propCiHi <- max(localTest$conf.int)
866. localRes$propPval <- localTest$p.value
867. } else {
868. warning('Somthing strange happend - contact developer with reproducible example')
869. }
871. localRes$nUp <- localNumber$Freq[which( localNumber$Var1 == levels(localNumber$Var1)[1] )]
872. localRes$nDown <- localNumber$Freq[which( localNumber$Var1 == levels(localNumber$Var1)[2] )]
874. return(localRes)
875. }
876. )
878. consequenceBalance$propQval <- p.adjust(consequenceBalance$propPval, method = 'fdr')
879. consequenceBalance$Significant <- consequenceBalance$propQval < alpha
880. consequenceBalance$Significant <- factor(
881. consequenceBalance$Significant,
882. levels=c(FALSE,TRUE)
883. )
885. #Prep for plotting
886. consequenceBalance2 <- consequenceBalance[which(
887. (consequenceBalance$nUp + consequenceBalance$nDown) >= minEventsForPlotting
888. ),]
889. if(nrow(consequenceBalance2) == 0) {
890. stop('No features left for ploting after filtering with via "minEventsForPlotting" argument.')
891. }
893. consequenceBalance2$nTot <- consequenceBalance2$nDown + consequenceBalance2$nUp
895. ### Add comparison
896. consequenceBalance2$Comparison <- paste(
897. consequenceBalance2$condition_1,
898. 'vs',
899. consequenceBalance2$condition_2,
900. sep='\n'
901. )
903. ### Massage
904. consequenceBalance2$feature2 <- gsub(' \\(', '\n(', consequenceBalance2$feature)
906. consequenceBalance2$feature2 <- factor(
907. consequenceBalance2$feature2,
908. levels = rev(sort(unique(as.character(consequenceBalance2$feature2))))
909. )
911. # SPLICING PLOT -----------------------------------------------------------
912. extractSwitchPairs <- function(
913. switchAnalyzeRlist,
914. alpha = 0.05,
915. dIFcutoff = 0.05,
916. onlySigIsoforms = FALSE
917. ) {
918. ### Extract and massage data
919. if (TRUE) {
920. localData <- switchAnalyzeRlist$isoformFeatures[
921. which(
922. switchAnalyzeRlist$isoformFeatures$gene_switch_q_value < alpha &
923. abs(switchAnalyzeRlist$isoformFeatures$dIF) > dIFcutoff
924. ),
925. c(
926. 'iso_ref',
927. 'gene_ref',
928. 'isoform_switch_q_value',
929. 'gene_switch_q_value',
930. 'dIF'
931. )
932. ]
934. if (!nrow(localData)) {
935. stop('No genes were considered switching with the used cutoff values')
936. }
938. ### add switch direction
939. localData$switchDirection <- NA
940. localData$switchDirection[which(sign(localData$dIF) == 1)] <- 'up'
941. localData$switchDirection[which(sign(localData$dIF) == -1)] <- 'down'
943. ### Annotate significant features
944. isoResTest <-
945. any(!is.na(
946. switchAnalyzeRlist$isoformFeatures$isoform_switch_q_value
947. ))
948. if (isoResTest) {
949. localData$isoSig <-
950. localData$isoform_switch_q_value < alpha &
951. abs(localData$dIF) > dIFcutoff
952. } else {
953. localData$isoSig <-
954. localData$gene_switch_q_value < alpha &
955. abs(localData$dIF) > dIFcutoff
956. }
959. if(onlySigIsoforms) {
960. localData <- localData[which( localData$isoSig ),]
961. }
962. }
965. ### Create data sub-sets of interest
966. if(TRUE) {
967. sigUpData <- localData[which(
968. localData$isoSig & localData$switchDirection == 'up'
969. ),c('iso_ref','gene_ref')]
970. sigDnData <- localData[which(
971. localData$isoSig & localData$switchDirection == 'down'
972. ),c('iso_ref','gene_ref')]
974. colnames(sigUpData)[1] <- c('iso_ref_up')
975. colnames(sigDnData)[1] <- c('iso_ref_down')
978. if( ! onlySigIsoforms ) {
979. justUpData <- localData[which(
980. localData$switchDirection == 'up'
981. ),c('iso_ref','gene_ref')]
982. justDnData <- localData[which(
983. localData$switchDirection == 'down'
984. ),c('iso_ref','gene_ref')]
986. colnames(justUpData)[1] <- c('iso_ref_up')
987. colnames(justDnData)[1] <- c('iso_ref_down')
988. }
989. }
991. ### Join datasets to extract pairs
992. if(TRUE) {
993. if( onlySigIsoforms ) {
994. pairwiseIsoComparison <- dplyr::inner_join(
995. sigUpData,
996. sigDnData,
997. by= 'gene_ref',
998. multiple = "all"
999. )
1000. } else {
1001. ### Sig up and all down
1002. upPairs <- dplyr::inner_join(
1003. sigUpData,
1004. justDnData,
1005. by= 'gene_ref',
1006. multiple = "all"
1007. )
1008. ### Sig down and all up
1009. dnPairs <- dplyr::inner_join(
1010. justUpData,
1011. sigDnData,
1012. by= 'gene_ref',
1013. multiple = "all"
1014. )
1016. ### Combine
1017. pairwiseIsoComparison <- unique(
1018. rbind(
1019. upPairs,
1020. dnPairs
1021. )
1022. )
1023. pairwiseIsoComparison <- pairwiseIsoComparison[,c(
1024. 'gene_ref','iso_ref_up','iso_ref_down'
1025. )]
1027. ### Reorder
1028. pairwiseIsoComparison <- pairwiseIsoComparison[order(
1029. pairwiseIsoComparison$gene_ref,
1030. pairwiseIsoComparison$iso_ref_up,
1031. pairwiseIsoComparison$iso_ref_down
1032. ),]
1033. }
1034. }
1036. ### Add in additional data
1037. if(TRUE) {
1038. ### Add isoform names
1039. matchVectorUp <- match(
1040. pairwiseIsoComparison$iso_ref_up,
1041. switchAnalyzeRlist$isoformFeatures$iso_ref
1042. )
1043. matchVectorDn <- match(
1044. pairwiseIsoComparison$iso_ref_down,
1045. switchAnalyzeRlist$isoformFeatures$iso_ref
1046. )
1048. pairwiseIsoComparison$isoformUpregulated <-
1049. switchAnalyzeRlist$isoformFeatures$isoform_id[matchVectorUp]
1050. pairwiseIsoComparison$isoformDownregulated <-
1051. switchAnalyzeRlist$isoformFeatures$isoform_id[matchVectorDn]
1053. ### Gene infl
1054. pairwiseIsoComparison$gene_id <-
1055. switchAnalyzeRlist$isoformFeatures$gene_id[matchVectorUp]
1056. pairwiseIsoComparison$gene_name <-
1057. switchAnalyzeRlist$isoformFeatures$gene_name[matchVectorUp]
1059. ### Conditons
1060. pairwiseIsoComparison$condition_1 <-
1061. switchAnalyzeRlist$isoformFeatures$condition_1[matchVectorUp]
1062. pairwiseIsoComparison$condition_2 <-
1063. switchAnalyzeRlist$isoformFeatures$condition_2[matchVectorUp]
1065. }
1067. return( pairwiseIsoComparison )
1068. }
1070. #Use MIC IsoformSwitch object to extract the MICs
1071. aSwitchList_part2_MIC <- readRDS(file="./code/microexons/output/isoformswitch_MIC.rds")
1073. switchAnalyzeRlist = aSwitchList_part2_MIC
1074. splicingToAnalyze = 'all'
1075. alpha = 0.05
1076. dIFcutoff = 0.05
1077. onlySigIsoforms = T
1078. countGenes = T
1079. plot = TRUE
1080. localTheme = theme_bw(base_size = 14)
1081. minEventsForPlotting = 25 ##default = 10
1082. returnResult=TRUE
1083. returnSummary=TRUE
1085. ### Consequences to analyze
1086. acceptedTypes <- c("A3","A5","ATSS","ATTS","ES" ,"IR","MEE","MES", "MIC" )
1089. splicingAnalyzed <-
1090. intersect(
1091. acceptedTypes,
1092. colnames(switchAnalyzeRlist$AlternativeSplicingAnalysis)
1093. )
1094. if ('all' %in% splicingToAnalyze) {
1095. splicingToAnalyze <- splicingAnalyzed
1096. }
1098. splicingNotAnalyzed <-
1099. setdiff(splicingToAnalyze, splicingAnalyzed)
1100. if (length(splicingNotAnalyzed)) {
1101. warning(
1102. paste(
1103. 'The following consequences appear not to have been analyzed and will therefor not be summarized:',
1104. paste(splicingNotAnalyzed, collapse = ', '),
1105. sep = ' '
1106. )
1107. )
1108. }
1110. ### Get pairs
1111. if(TRUE) {
1112. pairwiseIsoComparison <- extractSwitchPairs(
1113. switchAnalyzeRlist,
1114. alpha = alpha,
1115. dIFcutoff = dIFcutoff,
1116. onlySigIsoforms = onlySigIsoforms
1117. )
1118. }
1120. ### Massage AS analysis
1121. if(TRUE) {
1122. localAS <- switchAnalyzeRlist$AlternativeSplicingAnalysis
1123. localAS <- localAS[which(
1124. localAS$isoform_id %in% pairwiseIsoComparison$isoformUpregulated |
1125. localAS$isoform_id %in% pairwiseIsoComparison$isoformDownregulated
1126. ),]
1128. ### Massage
1129. m1 <- reshape2::melt(localAS[,c(
1130. "isoform_id",
1131. "ES_genomic_start",
1132. "MEE_genomic_start",
1133. "MES_genomic_start",
1134. "IR_genomic_start",
1135. "A5_genomic_start",
1136. "A3_genomic_start",
1137. "ATSS_genomic_start",
1138. "ATTS_genomic_start",
1139. "MIC_genomic_start"
1140. )], id.vars = 'isoform_id')
1141. colnames(m1)[3] <- 'genomic_start'
1142. m1$AStype <- sapply(
1143. strsplit(as.character(m1$variable),'_'),
1144. function(x) x[1]
1145. )
1147. m2 <- reshape2::melt(localAS[,c(
1148. "isoform_id",
1149. "ES_genomic_end",
1150. "MEE_genomic_end",
1151. "MES_genomic_end",
1152. "IR_genomic_end",
1153. "A5_genomic_end",
1154. "A3_genomic_end",
1155. "ATSS_genomic_end",
1156. "ATTS_genomic_end",
1157. "MIC_genomic_end"
1158. )], id.vars = 'isoform_id')
1159. colnames(m2)[3] <- 'genomic_end'
1160. m2$AStype <- sapply(
1161. strsplit(as.character(m2$variable),'_'),
1162. function(x) x[1]
1163. )
1165. localAS <- dplyr::inner_join(
1166. m1[,c('isoform_id','AStype','genomic_start')],
1167. m2[,c('isoform_id','AStype','genomic_end')],
1168. by=c('isoform_id','AStype')
1169. )
1171. ### Add in NMD
1172. if('orfAnalysis' %in% names(switchAnalyzeRlist)) {
1173. localNMD <- data.frame(
1174. isoform_id=switchAnalyzeRlist$orfAnalysis$isoform_id,
1175. AStype='NMD',
1176. genomic_start=ifelse(switchAnalyzeRlist$orfAnalysis$PTC, '0,0', '0'),
1177. genomic_end=ifelse(switchAnalyzeRlist$orfAnalysis$PTC, '0,0', '0'),
1178. stringsAsFactors = FALSE
1179. )
1180. localAS <- rbind(localAS , localNMD)
1181. }
1182. }
1184. ### Add AS to pairs
1185. if(TRUE) {
1186. localConseq2 <- merge(
1187. pairwiseIsoComparison,
1188. localAS,
1189. by.x='isoformUpregulated',
1190. by.y='isoform_id'
1191. )
1193. localConseq3 <- merge(
1194. localConseq2,
1195. localAS,
1196. by.x=c('isoformDownregulated','AStype'),
1197. by.y=c('isoform_id','AStype'),
1198. suffixes = c("_up","_down")
1199. )
1202. ### Replace NAs so they can be compared (na just mean same as pre-transcript)
1203. localConseq3$genomic_start_up[which(
1204. is.na(localConseq3$genomic_start_up)
1205. )] <- 0
1206. localConseq3$genomic_end_up[which(
1207. is.na(localConseq3$genomic_end_up)
1208. )] <- 0
1209. localConseq3$genomic_start_down[which(
1210. is.na(localConseq3$genomic_start_down)
1211. )] <- 0
1212. localConseq3$genomic_end_down[which(
1213. is.na(localConseq3$genomic_end_down)
1214. )] <- 0
1216. ### Identify differences
1217. localConseq3$coordinatsDifferent <-
1218. localConseq3$genomic_start_up != localConseq3$genomic_start_down |
1219. localConseq3$genomic_end_up != localConseq3$genomic_end_down
1221. localConseq4 <- localConseq3[which(localConseq3$coordinatsDifferent),]
1223. if(nrow(localConseq4) == 0) {
1224. stop('No alternative splicing differences were found')
1225. }
1227. }
1229. ### Subset based on input
1230. if(TRUE) {
1231. localConseq4 <- localConseq4[which(
1232. localConseq4$AStype %in% splicingToAnalyze
1233. ),]
1234. if (!nrow(localConseq4)) {
1235. stop('No swithces with consequences were found')
1236. }
1238. }
1240. ### Analyze differences (both start and end coordinats)
1241. if(TRUE) {
1242. genomic_start_up <- strsplit(x = localConseq4$genomic_start_up , split = ';')
1243. genomic_start_down <- strsplit(x = localConseq4$genomic_start_down, split = ';')
1244. genomic_end_up <- strsplit(x = localConseq4$genomic_end_up , split = ';')
1245. genomic_end_down <- strsplit(x = localConseq4$genomic_end_down, split = ';')
1247. localConseq4$nrGain <- 0
1248. localConseq4$nrLoss <- 0
1249. for(i in seq_along(genomic_start_up)) {
1250. # combine start and end of exons
1251. localUp <- stringr::str_c(genomic_start_up[[i]] , '_', genomic_end_up [[i]])
1252. localDn <- stringr::str_c(genomic_start_down[[i]], '_', genomic_end_down[[i]])
1254. # remove coordinats from primary transcripts
1255. localUp <- localUp[which( localUp != '0_0')]
1256. localDn <- localDn[which( localDn != '0_0')]
1258. # calculate difference
1259. localConseq4$nrGain[i] <- sum(!localUp %in% localDn)
1260. localConseq4$nrLoss[i] <- sum(!localDn %in% localUp)
1261. }
1262. ### Modify ends (can only have 1)
1263. terminiIndex <- which( localConseq4$AStype %in% c('ATSS','ATTS') )
1264. localConseq4$nrGain[terminiIndex] <-
1265. 1 * sign(localConseq4$nrGain[terminiIndex])
1266. localConseq4$nrLoss[terminiIndex] <-
1267. 1 * sign(localConseq4$nrLoss[terminiIndex])
1269. localConseq4$nrDiff <- localConseq4$nrGain - localConseq4$nrLoss
1271. }
1273. ### Summarize gain vs loss for each AStype in each condition
1274. gainLossBalance <- plyr::ddply(
1275. .data = localConseq4,
1276. .variables = c('condition_1','condition_2','AStype'),
1277. .fun = function(aDF) { # aDF <- localConseq4[1:50,]
1278. if( countGenes ) {
1279. aDF2 <- plyr::ddply(aDF[,c('gene_ref','nrDiff')], .variables = 'gene_ref', function(aGene) {
1280. data.frame(
1281. nrDiff = sum(aGene$nrDiff)
1282. )
1283. })
1285. localRes <- data.frame(
1286. nUp = sum(aDF2$nrDiff > 0),
1287. nDown= sum(aDF2$nrDiff < 0)
1288. )
1289. } else {
1290. localRes <- data.frame(
1291. nUp=sum(aDF$nrDiff > 0),
1292. nDown=sum(aDF$nrDiff < 0)
1293. )
1294. }
1296. if( localRes$nUp > 0 | localRes$nDown > 0 ) {
1297. localTest <- suppressWarnings(
1298. stats::binom.test(localRes$nUp, localRes$nUp + localRes$nDown)
1299. )
1301. localRes$propUp <- localTest$estimate
1302. localRes$propUpCiLo <- min(localTest$conf.int)
1303. localRes$propUpCiHi <- max(localTest$conf.int)
1304. localRes$propUpPval <- localTest$p.value
1305. } else {
1306. localRes$propUp <- NA
1307. localRes$propUpCiLo <- NA
1308. localRes$propUpCiHi <- NA
1309. localRes$propUpPval <- NA
1310. }
1312. return(localRes)
1313. }
1314. )
1316. ### Massage for plotting
1317. if(TRUE) {
1318. gainLossBalance <- gainLossBalance[which(
1319. ! is.na(gainLossBalance$propUp)
1320. ),]
1321. gainLossBalance$propUpQval <- p.adjust(gainLossBalance$propUpPval, method = 'fdr')
1322. gainLossBalance$Significant <- gainLossBalance$propUpQval < alpha
1323. gainLossBalance$Significant <- factor(
1324. gainLossBalance$Significant,
1325. levels=c(FALSE,TRUE)
1326. )
1328. gainLossBalance$Comparison <- paste(
1329. gainLossBalance$condition_1,
1330. 'vs',
1331. gainLossBalance$condition_2,
1332. sep='\n'
1333. )
1335. ### Massage splicing type name
1336. if(TRUE) {
1337. gainLossBalance$AStype <- paste0(
1338. gainLossBalance$AStype, ' gain ',
1339. '(paired with ',gainLossBalance$AStype, ' loss)'
1340. )
1342. gainLossBalance$AStype <- gsub('MES gain', 'MES', gainLossBalance$AStype)
1343. gainLossBalance$AStype <- gsub('MES loss', 'MEI', gainLossBalance$AStype)
1345. gainLossBalance$AStype <- gsub('ES gain', 'ES', gainLossBalance$AStype)
1346. gainLossBalance$AStype <- gsub('ES loss', 'EI', gainLossBalance$AStype)
1348. #gainLossBalance$AStype <- gsub('IR gain', 'IR', gainLossBalance$AStype)
1349. #gainLossBalance$AStype <- gsub('IR loss', 'IS', gainLossBalance$AStype)
1350. }
1352. myOrder <- plyr::ddply(
1353. gainLossBalance,
1354. .variables = 'AStype',
1355. .fun = function(aDF) {
1356. mean(aDF$propUp)
1357. })
1358. myOrder <- myOrder$AStype[sort.list(myOrder$V1, decreasing = TRUE)]
1360. gainLossBalance$AStype <- factor(
1361. gainLossBalance$AStype,
1362. levels = myOrder
1363. )
1364. }
1365. #
1366. # Prep data for plotting
1367. gainLossBalance2 <- gainLossBalance[which(
1368. (gainLossBalance$nUp + gainLossBalance$nDown) >= minEventsForPlotting
1369. ),]
1371. gainLossBalance2$nTot <- gainLossBalance2$nUp + gainLossBalance2$nDown
1374. ####Combining the data to plot:
1375. opposite <- extractConsequenceEnrichment(
1376. aSwitchList_part2,
1377. alpha=0.05,
1378. dIFcutoff = 0.05,
1379. countGenes = TRUE,
1380. analysisOppositeConsequence=T,
1381. plot=F,
1382. returnResult=TRUE,
1383. returnSummary=TRUE)
1384. opposite <- opposite %>% filter(conseqPair == "NMD_status" |
1385. conseqPair == "domains_identified" |
1386. conseqPair == "IDR_identified")
1387. opposite$Conseq_type <- "Functional Consequence"
1388. opposite$nTot <- opposite$nUp + opposite$nDown
1389. opposite$Comparison <- paste0(opposite$condition_1, " vs ", opposite$condition_2)
1390. opposite <- opposite %>% relocate(c(nUp, nDown), .after = feature)
1391. opposite <- opposite %>% relocate(Comparison, .before = nTot)
1392. opposite <- opposite[,-c(3)]
1393. opposite <- opposite %>% relocate(c(Conseq_type), .after = nTot)
1396. combined_data <- consequenceBalance2[,-c(3,15)]
1397. combined_data$Conseq_type <- "Functional Consequence"
1398. combined_data <- combined_data %>% relocate(c(nUp, nDown), .after = feature)
1399. combined_data <- combined_data %>% relocate(Comparison, .before = nTot)
1400. combined_data <- combined_data %>% filter(feature != "NMD insensitive (paired with NMD sensitive)" &
1401. feature != "Domain gain (paired with Domain loss)" &
1402. feature != "IDR gain (paired with IDR loss)")
1405. colnames(combined_data) <- c("condition_1", "condition_2", "feature", "nUp", "nDown",
1406. "propUp" , "propUpCiLo", "propUpCiHi", "propUpPval",
1407. "propUpQval", "Significant", "Comparison" , "nTot", "Conseq_type")
1408. gainLossBalance2$Conseq_type <- "Alternative Splicing"
1409. colnames(gainLossBalance2) <- c("condition_1", "condition_2", "feature", "nUp", "nDown",
1410. "propUp" , "propUpCiLo", "propUpCiHi", "propUpPval",
1411. "propUpQval", "Significant", "Comparison" , "nTot", "Conseq_type")
1412. colnames(opposite) <- c("condition_1", "condition_2", "feature", "nUp", "nDown",
1413. "propUp" , "propUpCiLo", "propUpCiHi", "propUpPval",
1414. "propUpQval", "Significant", "Comparison" , "nTot", "Conseq_type")
1415. combined_data <- rbind(combined_data, gainLossBalance2)
1416. combined_data <- rbind(combined_data, opposite)
1418. combined_data$Comparison <- str_replace(combined_data$Comparison, "t00\nvs\nt04", "t00 vs t04")
1419. combined_data$Comparison <- str_replace(combined_data$Comparison, "t04\nvs\nt30", "t04 vs t30")
1420. combined_data$Comparison <- str_replace(combined_data$Comparison, "t00\nvs\nt30", "t00 vs t30")
1422. combined_data$feature <- as.factor(combined_data$feature)
1423. library(forcats)
1426. combined_data$feature <- str_replace(combined_data$feature, "3UTR is longer \\(paired with 3UTR is shorter\\)",
1427. "3' UTR is longer")
1428. combined_data$feature <- str_replace(combined_data$feature, "5UTR is longer \\(paired with 5UTR is shorter\\)",
1429. "5' UTR is longer")
1430. combined_data$feature <- str_replace(combined_data$feature, "Domain loss \\(paired with Domain gain\\)",
1431. "Domain loss")
1432. combined_data$feature <- str_replace(combined_data$feature, "Domain length gain \\(paired with Domain length loss\\)",
1433. "Domain length gain")
1434. combined_data$feature <- str_replace(combined_data$feature, "Domain non-reference isotype gain \\(paired with Domain non-reference isotype loss\\)",
1435. "Domain non-reference isotype gain")
1436. combined_data$feature <- str_replace(combined_data$feature, "Exon gain \\(paired with Exon loss\\)",
1437. "Exon number gain")
1438. combined_data$feature <- str_replace(combined_data$feature, "IDR loss \\(paired with IDR gain\\)",
1439. "IDR loss")
1440. combined_data$feature <- str_replace(combined_data$feature, "IDR length gain \\(paired with IDR length loss\\)",
1441. "IDR is longer")
1442. combined_data$feature <- str_replace(combined_data$feature, "IDR w binding region gain \\(paired with IDR w binding region loss\\)",
1443. "IDR with binding region gain")
1444. combined_data$feature <- str_replace(combined_data$feature, "Last exon more downstream \\(paired with Last exon more upstream\\)",
1445. "Last exon more downstream")
1446. combined_data$feature <- str_replace(combined_data$feature, "Length gain \\(paired with Length loss\\)",
1447. "Isoform is longer")
1448. combined_data$feature <- str_replace(combined_data$feature, "NMD sensitive \\(paired with NMD insensitive\\)",
1449. "NMD sensitive")
1450. combined_data$feature <- str_replace(combined_data$feature, "ORF is longer \\(paired with ORF is shorter\\)",
1451. "ORF is longer")
1452. combined_data$feature <- str_replace(combined_data$feature, "Tss more downstream \\(paired with Tss more upstream\\)",
1453. "TSS more downstream")
1454. combined_data$feature <- str_replace(combined_data$feature, "Tts more downstream \\(paired with Tts more upstream\\)",
1455. "TTS more downstream")
1456. combined_data$feature <- str_replace(combined_data$feature, "A3 gain \\(paired with A3 loss\\)",
1457. "3' splice site more downstream")
1458. combined_data$feature <- str_replace(combined_data$feature, "A5 gain \\(paired with A5 loss\\)",
1459. "5' splice site more upstream")
1460. combined_data$feature <- str_replace(combined_data$feature, "ATSS gain \\(paired with ATSS loss\\)",
1461. "ATSS gain")
1462. combined_data$feature <- str_replace(combined_data$feature, "ATTS gain \\(paired with ATTS loss\\)",
1463. "ATTS gain")
1464. combined_data$feature <- str_replace(combined_data$feature, "ES \\(paired with EI\\)",
1465. "Exon skipping")
1466. combined_data$feature <- str_replace(combined_data$feature, "MIC gain \\(paired with MIC loss\\)",
1467. "Microexon skipping")
1468. combined_data$feature <- str_replace(combined_data$feature, "IR gain \\(paired with IR loss\\)",
1469. "Intron retention")
1470. combined_data$feature <- str_replace(combined_data$feature, "MES \\(paired with MEI\\)",
1471. "Multiple exon skipping")
1473. combined_data$feature <- str_replace(combined_data$feature, "Signal peptide gain \\(paired with Signal peptide loss\\)",
1474. "Signal peptide gain")
1475. combined_data$feature <- str_replace(combined_data$feature, "SubCell location cyto gain \\(paired with SubCell location cyto loss\\)",
1476. "Shift to cytoplasm localization")
1477. combined_data$feature <- str_replace(combined_data$feature, "SubCell location ext cell gain \\(paired with SubCell location ext cell loss\\)",
1478. "Shift to extracellular localization")
1479. combined_data$feature <- str_replace(combined_data$feature, "SubCell location gain \\(paired with SubCell location loss\\)",
1480. "Subcellular localization gain")
1481. combined_data$feature <- str_replace(combined_data$feature, "SubCell location memb gain \\(paired with SubCell location memb loss\\)",
1482. "Shift to cell membrane localization")
1483. combined_data$feature <- str_replace(combined_data$feature, "SubCell location nucl gain \\(paired with SubCell location nucl loss\\)",
1484. "Shift to nuclear localization")
1486. combined_data$feature <- str_replace(combined_data$feature, "Topology complexity gain \\(paired with Topology complexity loss\\)",
1487. "Topology complexity gain")
1488. combined_data$feature <- str_replace(combined_data$feature, "Extracellular length gain \\(paired with Extracellular length loss\\)",
1489. "Extracellular region length gain")
1490. combined_data$feature <- str_replace(combined_data$feature, "Extracellular region gain \\(paired with Extracellular region loss\\)",
1491. "Extracellular region gain")
1492. combined_data$feature <- str_replace(combined_data$feature, "Intracellular length gain \\(paired with Intracellular length loss\\)",
1493. "Intracellular region length gain")
1494. combined_data$feature <- str_replace(combined_data$feature, "Intracellular region gain \\(paired with Intracellular region loss\\)",
1495. "Intracellular region gain")
1499. combined_data <- combined_data %>%
1500. mutate(feature = fct_relevel(feature,
1501. "3' splice site more downstream", "5' splice site more upstream",
1502. "ATSS gain", "ATTS gain", "Exon skipping", "Microexon skipping", "Intron retention", "Multiple exon skipping",
1503. # "Intracellular region gain",
1504. # "Extracellular region gain",
1505. # "Intracellular region length gain",
1506. # "Extracellular region length gain",
1507. # "Topology complexity gain",
1508. "Shift to extracellular localization",
1509. "Shift to cell membrane localization",
1510. "Shift to nuclear localization",
1511. "Shift to cytoplasm localization",
1512. "Subcellular localization gain",
1513. "Signal peptide gain",
1514. "IDR with binding region gain",
1515. "IDR is longer",
1516. "IDR loss",
1517. "Domain non-reference isotype gain",
1518. "Domain length gain",
1519. "Domain loss",
1520. "NMD sensitive",
1521. "ORF is longer",
1522. "Exon number gain",
1523. "Isoform is longer",
1524. "Last exon more downstream",
1525. "3' UTR is longer",
1526. "TTS more downstream",
1527. "5' UTR is longer",
1528. "TSS more downstream"))
1530. colnames(combined_data)[13] <- 'Genes'
1532. library(ggh4x)
1533. library(extrafont)
1534. library(scales)
1535. library(grid)
1537. font_import(paths = "C:/Windows/Fonts", pattern = "cambria", prompt = F) # This imports all the fonts from your system
1538. loadfonts(device = "win") # For Windows
1539. windowsFonts()
1541. plot_data <- combined_data %>% filter(#feature == "Multiple exon skipping" |
1542. feature == "Intron retention" |
1543. feature == "Exon skipping" |
1544. feature == "Microexon skipping" |
1545. # feature == "5' splice site more upstream" |
1546. # feature == "3' splice site more downstream" |
1547. # feature == "TSS more downstream" |
1548. feature == "5' UTR is longer" |
1549. # feature == "TTS more downstream" |
1550. feature == "3' UTR is longer" |
1551. # feature == "Isoform is longer" |
1552. # feature == "ORF is longer" |
1553. feature == "NMD sensitive" |
1554. feature == "Domain loss" |
1555. # feature == "Domain length gain" |
1556. feature == "IDR loss" #|
1557. # feature == "IDR is longer"
1558. )
1559. desired_order <- c("Exon skipping",
1560. "Microexon skipping",
1561. "Intron retention",
1562. "NMD sensitive",
1563. "Domain loss",
1564. "IDR loss",
1565. "5' UTR is longer",
1566. "3' UTR is longer")
1567. plot_data <- plot_data %>%
1568. mutate(feature = factor(feature, levels = rev(desired_order)))
1571. combined_plot_1 <- plot_data %>%
1572. ggplot(aes(y=feature, x=propUp, color=Significant)) +
1573. geom_errorbarh(aes(xmax = propUpCiLo, xmin=propUpCiHi), height = .3) +
1574. # geom_point(aes(size=Genes)) +
1575. geom_text(label = "\u25D6", aes(size=Genes), family = "Cambria") +
1576. facet_wrap(~ factor(Comparison, c("t00 vs t04", "t04 vs t30", "t00 vs t30")), scales = "free_x") +
1577. geom_vline(xintercept=0.5, linetype='dashed') +
1578. labs(title = "Isoform Switching Consequences",
1579. x="Fraction of Genes with Switches Primarily Resulting in Each Event") +
1580. theme_void(base_size = 8) +
1581. theme(plot.title = element_text(hjust = 0.5, size = 8, face = "bold"),
1582. strip.text = element_text(size = 8),
1583. axis.text = element_text(colour = "black", size = 6.4),
1584. axis.text.x=element_text(angle=-45, hjust = 0, vjust=1),
1585. legend.position = "top",
1586. legend.box = "horizontal",
1587. legend.box.just = "center",
1588. legend.text = element_text(size = 7),
1589. legend.key.spacing.y = unit(0, "pt"),
1590. legend.key.height = unit(-20, "lines"),
1591. axis.title.y = element_blank(),
1592. legend.margin = margin(
1593. t = 0,
1594. r = 0,
1595. b = 0,
1596. l = 0,
1597. unit = "cm"
1598. )) +
1599. scale_color_manual(name = paste0('FDR < ', alpha), values=c('black','red'), drop=FALSE) +
1600. scale_radius(limits=c(1, max(plot_data$Genes)),
1601. breaks = c(50, 500, 2000, 4000)) +
1602. guides(color = guide_legend(order=1, nrow = 2, byrow = T),
1603. size = guide_legend(order=2, override.aes = list(color = "black"), nrow = 2, byrow = T)) +
1604. coord_cartesian(xlim=c(0,1))
1606. combined_plot_1
1608. # Create a version of the plot for legend extraction
1609. combined_plot_symbols_only <- combined_plot_1 +
1610. theme(
1611. legend.text = element_text(colour = "transparent"), # invisible but still occupies space
1612. legend.title = element_text(colour = "transparent") # invisible but still occupies space
1613. )
1615. legend_1 <- ggpubr::get_legend(combined_plot_symbols_only)
1619. combined_plot_2 <- plot_data %>%
1620. ggplot(aes(y=feature, x=propUp, color=Significant)) +
1621. geom_errorbarh(aes(xmax = propUpCiLo, xmin=propUpCiHi), height = .3) +
1622. geom_text(label = "\u25D7", aes(size=Genes), family = "Cambria") +
1623. facet_wrap(~ factor(Comparison, c("t00 vs t04", "t04 vs t30", "t00 vs t30")), scales = "free_x") +
1624. geom_vline(xintercept=0.5, linetype='dashed') +
1625. labs(title = "Isoform Switching Consequences",
1626. x="Fraction of Genes with Switches Primarily Resulting in Each Event"
1627. ) +
1628. theme_void(base_size = 8) +
1629. theme(plot.title = element_text(hjust = 0.5, size = 8, face = "bold"),
1630. strip.text = element_text(size = 8),
1631. axis.text = element_text(colour = "black", size = 6.4),
1632. axis.text.x=element_text(angle=-45, hjust = 0, vjust=1),
1633. legend.position = "top",
1634. legend.box = "horizontal",
1635. legend.box.just = "center",
1636. legend.key.height = unit(-10, "lines"),
1637. legend.key.spacing.y = unit(0, "pt"),
1638. legend.text = element_text(size = 7),
1639. axis.title.y = element_blank(),
1640. legend.margin = margin(
1641. t = 0,
1642. r = 0,
1643. b = 0,
1644. l = 0,
1645. unit = "cm"
1646. )) +
1647. scale_color_manual(name = paste0('FDR < ', alpha), values=c('black','red'), drop=FALSE) +
1648. scale_radius(limits=c(1, max((plot_data$Genes))),
1649. breaks = c(50, 500, 2000, 4000)) +
1650. guides(color = guide_legend(order=1, nrow = 2, byrow = T),
1651. size = guide_legend(order=2, override.aes = list(color = "red"), nrow = 2, byrow = T)) +
1652. coord_cartesian(xlim=c(0,1))
1654. combined_plot_2
1657. # Extract the legend from the plot
1658. legend_2 <- ggpubr::get_legend(combined_plot_2)
1660. combined_plot <- plot_data %>%
1661. ggplot(aes(y=feature, x=propUp, color=Significant)) +
1662. geom_errorbarh(aes(xmax = propUpCiLo, xmin=propUpCiHi), height = .3) +
1663. # geom_point(aes(size=Genes)) +
1664. geom_text(label = "\u25CF", aes(size=Genes), family = "Cambria") +
1665. facet_wrap(~ factor(Comparison, c("t00 vs t04", "t04 vs t30", "t00 vs t30")), scales = "free_x") +
1666. geom_vline(xintercept=0.5, linetype='dashed') +
1667. labs(title = "Isoform Switching Consequences",
1668. x=str_wrap("Fraction of Genes with Switches Primarily Resulting in Each Event", 40)
1669. # y='Event/Consequence of Isoform Switch'
1670. ) +
1671. theme_bw(base_size = 8) +
1672. theme(plot.title = element_text(hjust = 0.5, vjust = 17, size = 8, face = "bold"),
1673. strip.text = element_text(size = 8),
1674. axis.text = element_text(colour = "black", size = 6.4),
1675. axis.text.x=element_text(angle=-45, hjust = 0, vjust=1),
1676. legend.position = "none",
1677. axis.title.y = element_blank() ) +
1678. scale_color_manual(name = paste0('FDR < ', alpha), values=c('black','red'), drop=FALSE) +
1679. scale_radius(limits=c(1, max((plot_data$Genes)))) +
1680. guides(color = guide_legend(order=1),
1681. size = guide_legend(order=2)) +
1682. coord_cartesian(xlim=c(0,1)) +
1683. scale_y_discrete(labels = label_wrap_gen(width = 17))
1685. combined_plot
1687. testp <- ggdraw() +
1688. draw_plot(combined_plot + theme(plot.margin = margin(t = 15, r = 6, b = 6, l = 10, unit = "pt"))
1689. , x = 0, y = 0, width = 1, height = 0.94) +
1690. draw_plot(legend_1, x = 0, y = 0.41, width = 1, height = 1) +
1691. draw_plot(legend_2, x = 0, y = 0.41, width = 1, height = 1)
1692. testp
1697. # PLOTTING SWITCH SUPP FIG ------------------------------------------------
1698. plot_data_full <- combined_data %>% filter(feature != "Shift to extracellular localization" &
1699. feature != "Shift to cell membrane localization" &
1700. feature != "Shift to nuclear localization" &
1701. feature != "Shift to cytoplasm localization" &
1702. feature != "Subcellular localization gain" &
1703. feature != "Signal peptide gain" &
1704. feature != "MEE gain (paired with MEE loss)" )
1707. combined_plot_1 <- plot_data_full %>%
1708. ggplot(aes(y=feature, x=propUp, color=Significant)) +
1709. geom_errorbarh(aes(xmax = propUpCiLo, xmin=propUpCiHi), height = .3) +
1710. # geom_point(aes(size=Genes)) +
1711. geom_text(label = "\u25D6", aes(size=Genes), family = "Cambria") +
1712. facet_grid(Conseq_type ~ factor(Comparison, c("t00 vs t04", "t04 vs t30", "t00 vs t30")), scales = "free_y") +
1713. force_panelsizes(c((plot_data %>% filter(Conseq_type == "Alternative Splicing") %>% tally() /3),
1714. (plot_data %>% filter(Conseq_type == "Functional Consequence") %>% tally() /3)),
1715. cols = 3,
1716. respect = NULL) +
1717. geom_vline(xintercept=0.5, linetype='dashed') +
1718. labs(title = "Isoform Switching Consequences",
1719. x="Fraction of Genes with Switches Primarily Resulting in Each Event") +
1720. theme_void(base_size = 8) +
1721. theme(plot.title = element_text(hjust = 0.5, size = 8, face = "bold"),
1722. strip.text = element_text(size = 8),
1723. axis.text = element_text(colour = "black", size = 6.4),
1724. axis.text.x=element_text(angle=-45, hjust = 0, vjust=1),
1725. legend.position = "top",
1726. legend.box = "horizontal",
1727. legend.box.just = "center",
1728. legend.text = element_text(size = 7),
1729. # legend.key.size = element_rect(size = 10),
1730. legend.key.spacing.y = unit(0, "pt"),
1731. legend.key.height = unit(-20, "lines"),
1732. axis.title.y = element_blank(),
1733. legend.margin = margin(
1734. t = 0,
1735. r = 0,
1736. b = 0,
1737. l = 0,
1738. unit = "cm"
1739. )) +
1740. scale_color_manual(name = paste0('FDR < ', alpha), values=c('black','red'), drop=FALSE) +
1741. scale_radius(limits=c(1, 6501),
1742. breaks = c(50, 500, 2500, 6500)) +
1743. guides(color = guide_legend(order=1, nrow = 2, byrow = T),
1744. size = guide_legend(order=2, override.aes = list(color = "black"), nrow = 2, byrow = T)) +
1745. coord_cartesian(xlim=c(0,1))
1747. combined_plot_1
1750. # Create a version of the plot for legend extraction
1751. combined_plot_symbols_only <- combined_plot_1 +
1752. theme(
1753. legend.text = element_text(colour = "transparent"), # invisible but still occupies space
1754. legend.title = element_text(colour = "transparent") # invisible but still occupies space
1755. )
1757. legend_1 <- ggpubr::get_legend(combined_plot_symbols_only)
1761. combined_plot_2 <- plot_data_full %>%
1762. ggplot(aes(y=feature, x=propUp, color=Significant)) +
1763. geom_errorbarh(aes(xmax = propUpCiLo, xmin=propUpCiHi), height = .3) +
1764. geom_text(label = "\u25D7", aes(size=Genes), family = "Cambria") +
1765. facet_grid(Conseq_type ~ factor(Comparison, c("t00 vs t04", "t04 vs t30", "t00 vs t30")), scales = "free_y") +
1766. force_panelsizes(c((plot_data %>% filter(Conseq_type == "Alternative Splicing") %>% tally() /3),
1767. (plot_data %>% filter(Conseq_type == "Functional Consequence") %>% tally() /3)), #7 AS, 15 func conseq
1768. cols = 3,
1769. respect = NULL) +
1770. geom_vline(xintercept=0.5, linetype='dashed') +
1771. labs(title = "Isoform Switching Consequences",
1772. x="Fraction of Genes with Switches Primarily Resulting in Each Event"
1773. # y='Event/Consequence of Isoform Switch'
1774. ) +
1775. theme_void(base_size = 8) +
1776. theme(plot.title = element_text(hjust = 0.5, size = 8, face = "bold"),
1777. strip.text = element_text(size = 8),
1778. axis.text = element_text(colour = "black", size = 6.4),
1779. axis.text.x=element_text(angle=-45, hjust = 0, vjust=1),
1780. legend.position = "top",
1781. legend.box = "horizontal",
1782. legend.box.just = "center",
1783. legend.key.height = unit(-10, "lines"),
1784. legend.key.spacing.y = unit(0, "pt"),
1785. legend.text = element_text(size = 7),
1786. axis.title.y = element_blank(),
1787. legend.margin = margin(
1788. t = 0,
1789. r = 0,
1790. b = 0,
1791. l = 0,
1792. unit = "cm"
1793. )) +
1794. scale_color_manual(name = paste0('FDR < ', alpha), values=c('black','red'), drop=FALSE) +
1795. scale_radius(limits=c(1, 6501),
1796. breaks = c(50, 500, 2500, 6500)) +
1797. guides(color = guide_legend(order=1, nrow = 2, byrow = T),
1798. size = guide_legend(order=2, override.aes = list(color = "red"), nrow = 2, byrow = T)) +
1799. coord_cartesian(xlim=c(0,1))
1801. combined_plot_2
1804. # Extract the legend from the plot
1805. legend_2 <- ggpubr::get_legend(combined_plot_2)
1807. combined_plot <- plot_data_full %>%
1808. ggplot(aes(y=feature, x=propUp, color=Significant)) +
1809. geom_errorbarh(aes(xmax = propUpCiLo, xmin=propUpCiHi), height = .3) +
1810. # geom_point(aes(size=Genes)) +
1811. geom_text(label = "\u25CF", aes(size=Genes), family = "Cambria") +
1812. facet_grid(Conseq_type ~ factor(Comparison, c("t00 vs t04", "t04 vs t30", "t00 vs t30")), scales = "free_y") +
1813. force_panelsizes(c((plot_data %>% filter(Conseq_type == "Alternative Splicing") %>% tally() /3),
1814. (plot_data %>% filter(Conseq_type == "Functional Consequence") %>% tally() /3)), #7 AS, 15 func conseq
1815. cols = 3,
1816. respect = NULL) +
1817. geom_vline(xintercept=0.5, linetype='dashed') +
1818. labs(title = "Isoform Switching Consequences",
1819. x=str_wrap("Fraction of Genes with Switches Primarily Resulting in Each Event", 40)
1820. # y='Event/Consequence of Isoform Switch'
1821. ) +
1822. theme_bw(base_size = 8) +
1823. theme(plot.title = element_text(hjust = 0.5, vjust = 17, size = 8, face = "bold"),
1824. strip.text = element_text(size = 8),
1825. axis.text = element_text(colour = "black", size = 6.4),
1826. axis.text.x=element_text(angle=-45, hjust = 0, vjust=1),
1827. legend.position = "none",
1828. axis.title.y = element_blank() ) +
1829. scale_color_manual(name = paste0('FDR < ', alpha), values=c('black','red'), drop=FALSE) +
1830. scale_radius(limits=c(1, 6501)) +
1831. guides(color = guide_legend(order=1),
1832. size = guide_legend(order=2)) +
1833. coord_cartesian(xlim=c(0,1)) +
1834. scale_y_discrete(labels = label_wrap_gen(width = 22))
1836. combined_plot
1838. #ggarrange(combined_plot ,legend_2)
1839. testp <- ggdraw() +
1840. draw_plot(combined_plot + theme(plot.margin = margin(t = 15, r = 6, b = 6, l = 10, unit = "pt"))
1841. , x = 0, y = 0, width = 1, height = 0.94) +
1842. draw_plot(legend_1, x = 0, y = 0.42, width = 1, height = 1) +
1843. draw_plot(legend_2, x = 0, y = 0.42, width = 1, height = 1)
1844. testp
1846. library(Cairo)
1848. CairoPDF(paste0("./figures/SUPP_switching.pdf"), height = 7, width = 6.5)
1849. print(testp)
1850. dev.off()
1852. svg(paste0("./figures/SUPP_switching.svg"), height = 7, width = 6.5)
1853. print(testp)
1854. dev.off()
1857. # AGO1 PLOT -----------------------------------------------
1858. #Need to manually run these IsoformSwitchAnalyzeR functions to run the modified function below
1859. source('./code/figures/Functions.R')
1861. switchAnalyzeRlist = aSwitchList_part2
1862. gene = "AGO1"
1863. isoform_id = NULL
1864. condition1 = "t00"
1865. condition2 = "t30"
1867. ### Advanced arguments
1868. IFcutoff = 0.05
1869. dIFcutoff = 0.05
1870. alphas = c(0.05, 0.001)
1871. rescaleTranscripts = TRUE
1872. plotTopology = FALSE
1873. reverseMinus = TRUE
1874. addErrorbars = TRUE
1875. logYaxis = FALSE
1876. localTheme = theme_bw(base_size = 6.4)
1877. additionalArguments = list(switchPlotTranscript, switchPlotGeneExp, switchPlotIsoExp, switchPlotIsoUsage)
1879. # check conditions
1880. if (TRUE) {
1881. ### Identify conditions
1882. allConditionPairs <-
1883. unique(switchAnalyzeRlist$isoformFeatures[,c(
1884. 'condition_1', 'condition_2'
1885. )])
1887. levelsToMatch <- as.vector(t(allConditionPairs))
1890. conditionsFound <-
1891. apply(allConditionPairs, 1, function(x) {
1892. x[1] == condition1 & x[2] == condition2
1893. })
1895. if (!any(conditionsFound)) {
1896. # the the conditions are not found try revers order
1897. conditionsFound <-
1898. apply(allConditionPairs, 1, function(x) {
1899. x[2] == condition1 & x[1] == condition2
1900. })
1902. # if now found change the order of conditions
1903. if (any(conditionsFound)) {
1904. temp <- condition1
1905. condition1 <- condition2
1906. condition2 <- temp
1907. } else {
1908. stop(
1909. paste(
1910. 'The wanted comparison: \'',
1911. condition1,
1912. ' vs ',
1913. condition2,
1914. '\' is not in the data - please revise accordingly'
1915. )
1916. )
1917. }
1918. }
1920. } else {
1921. condition1 <- allConditionPairs$condition_1
1922. condition2 <- allConditionPairs$condition_2
1923. }
1925. levelsToMatch <-
1926. levelsToMatch[which(
1927. levelsToMatch %in% c(condition1, condition2)
1928. )]
1931. ### Parse additional list argument
1932. if (TRUE) {
1933. ### Expression plot arguments
1934. if (TRUE) {
1935. ## Make list with default values - can be replaced by as.list(args(expressionPlots)) ?
1936. eArgList <- list(
1937. switchAnalyzeRlist = switchAnalyzeRlist,
1938. gene = gene,
1939. isoform_id = isoform_id,
1940. condition1 = condition1,
1941. condition2 = condition2,
1942. IFcutoff = IFcutoff,
1943. addErrorbars = addErrorbars,
1944. confidenceIntervalErrorbars = TRUE,
1945. confidenceInterval = 0.95,
1946. #alphas = c(0.05, 0.001),
1947. alphas = alphas,
1948. logYaxis = logYaxis,
1949. extendFactor = 0.05,
1950. localTheme = localTheme
1951. )
1953. ### Modify default arguments if nessesary
1954. # subset to only those not already covered by the arguments
1955. additionalArguments2 <-
1956. additionalArguments[which(
1957. names(additionalArguments) %in% c(
1958. 'confidenceIntervalErrorbars',
1959. 'confidenceInterval',
1960. 'alphas',
1961. 'extendFactor',
1962. 'rescaleRoot'
1963. )
1964. )]
1965. if (length(additionalArguments2) != 0) {
1966. newArgNames <- names(additionalArguments2)
1968. for (i in seq_along(additionalArguments2)) {
1969. eArgList[[newArgNames[i]]] <- additionalArguments2[[i]]
1970. }
1971. }
1973. }
1977. }
1979. ### Make subplots - these functions also check the input thoroughly
1981. # expression plots
1982. expressionPlots <- expressionAnalysisPlot(
1983. switchAnalyzeRlist = eArgList$switchAnalyzeRlist,
1984. gene = eArgList$gene,
1985. isoform_id = eArgList$isoform_id,
1986. condition1 = eArgList$condition1,
1987. condition2 = eArgList$condition2,
1988. IFcutoff = eArgList$IFcutoff,
1989. addErrorbars = eArgList$addErrorbars,
1990. confidenceIntervalErrorbars = eArgList$confidenceIntervalErrorbars,
1991. confidenceInterval = eArgList$confidenceInterval,
1992. alphas = eArgList$alphas,
1993. optimizeForCombinedPlot = TRUE,
1994. logYaxis = eArgList$logYaxis,
1995. extendFactor = eArgList$extendFactor,
1996. localTheme = eArgList$localTheme
1997. )
1998. #
2000. ### Change plot margin
2001. marginSize <- 0.1
2002. expressionPlots2 <-
2003. lapply(expressionPlots[1:2], function(x) {
2004. x +
2005. theme(
2006. plot.margin = margin(
2007. t = 6,
2008. r = 6,
2009. b = 6,
2010. l = 6,
2011. unit = "pt" ),
2012. legend.position="none"
2013. )
2014. })
2016. expressionPlots3 <-
2017. lapply(expressionPlots[3], function(x) {
2018. x +
2019. theme(
2020. margin(
2021. t = 6,
2022. r = 6,
2023. b = 6,
2024. l = 6,
2025. unit = "pt" ),
2026. legend.position="none"
2027. )
2028. })
2030. label_new <- c(
2031. "<span style='color:#009E73'><b>PB.1207.309",
2032. "<span style='color:#009E73'><b>PB.1207.317",
2033. "<span style='color:#D55E00'><b>PB.1207.336")
2035. #Run this file to bring in the edited AGO1 plotting functions
2036. source("./code/figures/ago_plotting.R")
2037. myPlot <- AGO1_plotting()
2039. library(ggtext)
2040. test <-
2041. myPlot +
2042. theme(plot.margin = margin(
2043. t = 6,
2044. r = 6,
2045. b = 6,
2046. l = 6,
2047. unit = "pt"),
2048. legend.margin = margin(t= -6,0,0,0,unit = "pt"),
2049. legend.position="bottom",
2050. legend.key.size = unit(0.25, 'cm'),
2051. plot.title = element_text(hjust = 0.5, face = "bold", size = 10),
2052. axis.text.y = element_markdown(size = 8),
2053. legend.title = element_blank(),
2054. legend.text = element_text(size = 8)) +
2055. scale_y_continuous(
2056. breaks = c(3,2,1), # Set the y-axis positions
2057. labels = label_new # Apply the custom labels
2058. ) +
2059. labs(title = "AGO1 - t00 vs t30") +
2060. guides(fill = guide_legend(nrow = 1))
2061. test
2062. #
2063. expressionPlots2$gene_expression$data$Analyis <- "DGE"
2064. expressionPlots3$isoform_usage$data$Analyis <- "DTU"
2065. expressionPlots2$isoform_expression$data$Analyis <- "DTE"
2068. AGO1 <- ggdraw() +
2069. draw_plot(test + theme(panel.grid.major = element_blank(),
2070. panel.grid.minor = element_blank(),
2071. plot.title = element_text(hjust = 0.5, size = 8, face = "bold")), x = 0, y = 0.6, width = 1, height = 0.4) +
2072. draw_plot(expressionPlots2$gene_expression +
2073. theme(axis.text.x=element_text(angle=-45, hjust = 0, vjust=1, size = 6.4, colour = "black"),
2074. axis.text.y=element_text(size = 6.4, colour = "black"),
2075. axis.title = element_text(size = 8),
2076. strip.text = element_text(size = 8) ),
2077. x = 0, y = 0, width = 0.3, height = 0.6) +
2078. draw_plot(expressionPlots3$isoform_usage +
2079. theme(plot.margin = margin(t = 6, r = 6, b = 6, l = 6, unit = "pt"),
2080. axis.text.x=element_text(angle=-45, hjust = 0, vjust=1, size = 6.4, colour = "black"),
2081. axis.text.y=element_text(size = 6.4, colour = "black"),
2082. axis.title = element_text(size = 8),
2083. strip.text = element_text(size = 8)) +
2084. facet_wrap(~ Analyis, labeller = label_wrap_gen(width = 10)) , x = 0.65, y = 0, width = 0.35, height = 0.6) +
2085. draw_plot(expressionPlots2$isoform_expression +
2086. theme(axis.text.x=element_text(angle=-45, hjust = 0, vjust=1, size = 6.4, colour = "black"),
2087. axis.text.y=element_text(size = 6.4, colour = "black"),
2088. axis.title = element_text(size = 8),
2089. strip.text = element_text(size = 8)) +
2090. facet_wrap(~ Analyis, labeller = label_wrap_gen(width = 10)), x = 0.3, y = 0, width = 0.35, height = 0.6)
2091. AGO1
2093. AGO1 <- AGO1 + theme(plot.margin = margin(t = 0, r = 16, b = 0, l = 6, unit = "pt"))
2097. # IR + NMD ----------------------------------------------------------------
2098. #IR transcripts:
2099. AS <- aSwitchList_part2$AlternativeSplicingAnalysis
2100. AS_IR <- AS %>% filter(IR > 0) #12224 transcripts with 1 or more retained introns
2102. AS_IR <- merge(AS_IR, classification[, c("isoform", "associated_gene")],
2103. by.x = "isoform_id", by.y = "isoform", all.x = T)
2104. CPM <- read.csv("./code/expression/output/transcript_CPM.csv")
2106. CPM$Avgt00 <- rowMeans(subset(CPM, select = c(1:3)))
2107. CPM$Avgt04 <- rowMeans(subset(CPM, select = c(4:9)))
2108. CPM$Avgt30 <- rowMeans(subset(CPM, select = c(10:15)))
2110. CPM$IR <- ifelse(CPM$pb_id %in% AS_IR$isoform_id, "IR", "no_IR")
2112. #Only keep isoforms in the AS from isoswitch object:
2113. CPM_AS <- CPM %>% filter(pb_id %in% aSwitchList_part2$AlternativeSplicingAnalysis$isoform_id)
2117. #PTC transcripts:
2118. PTC_tr <- aSwitchList_part2$isoformFeatures %>% filter(PTC == TRUE) %>% pull(isoform_id) %>% unique()
2120. CPM$PTC <- ifelse(CPM$pb_id %in% PTC_tr, "PTC", "no_PTC")
2121. CPM_PTC <- CPM %>% filter(pb_id %in% aSwitchList_part2$AlternativeSplicingAnalysis$isoform_id)
2124. #IR =
2125. IR_tr <- CPM %>% filter(IR == "IR") %>% pull(pb_id) %>% unique()
2127. #IR only:
2128. IR_only <- setdiff(IR_tr, PTC_tr) #9606
2129. #NMD only:
2130. PTC_only <- setdiff(PTC_tr, IR_tr)
2131. #Overlap:
2132. overlap <- intersect(IR_tr, PTC_tr)
2133. #Remaining:
2134. all <- unique(union(IR_tr, PTC_tr))
2135. remaining <- setdiff(CPM_AS$pb_id, all)
2137. transcript_groups <- list(
2138. IR_only = IR_only,
2139. PTC_only = PTC_only,
2140. overlap = overlap,
2141. remaining = remaining)
2143. # Convert the list into a long data.frame for joining
2144. group_df <- bind_rows(lapply(names(transcript_groups), function(name) {
2145. data.frame(pb_id = transcript_groups[[name]], group = name)
2146. }))
2148. # Join and label the group
2149. CPM_AS <- CPM_AS %>%
2150. left_join(group_df, by = "pb_id")
2153. cpm_long_AS <- CPM_AS %>%
2154. pivot_longer(c("Avgt00", "Avgt04", "Avgt30"),
2155. names_to = "Timepoint",
2156. values_to = "CPM")
2158. library(rstatix)
2159. library(scales)
2160. cpm_long_AS$log_cpm <- log2(cpm_long_AS$CPM + 1)
2162. stat.test.ir <- cpm_long_AS %>% group_by(Timepoint) %>% wilcox_test(log_cpm ~ group)
2163. stat.test.ir <- stat.test.ir %>% mutate(p.adj.signif = p.adjust(p, method='bonferroni')) #need to do sep, b/c of group_by
2164. stat.test.ir <- stat.test.ir %>% add_xy_position(x = "Timepoint")
2166. stat.test.ir$p.adj.signif <- sprintf("%.2e", stat.test.ir$p.adj.signif)
2167. stat.test.ir <- stat.test.ir %>% add_significance(p.col = "p.adj")
2170. stat.test.ir$y.position <- log2(stat.test.ir$y.position + 1)
2172. stat.test.ir <- stat.test.ir %>%
2173. mutate(group1 = recode(group1, "IR_only" = "IR only"),
2174. group2 = recode(group2, "IR_only" = "IR only"),
2175. group1 = recode(group1, "overlap" = "IR & PTC"),
2176. group2 = recode(group2, "overlap" = "IR & PTC"),
2177. group1 = recode(group1, "PTC_only" = "PTC only"),
2178. group2 = recode(group2, "PTC_only" = "PTC only"),
2179. group1 = recode(group1, "remaining" = "Others"),
2180. group2 = recode(group2, "remaining" = "Others"))
2182. cpm_long_AS <- cpm_long_AS %>%
2183. mutate(group = recode(group, "IR_only" = "IR only"),
2184. group = recode(group, "IR_only" = "IR only"),
2185. group = recode(group, "overlap" = "IR & PTC"),
2186. group = recode(group, "overlap" = "IR & PTC"),
2187. group = recode(group, "PTC_only" = "PTC only"),
2188. group = recode(group, "PTC_only" = "PTC only"),
2189. group = recode(group, "remaining" = "Others"),
2190. group = recode(group, "remaining" = "Others"))
2192. cpm_long_AS$group <- factor(cpm_long_AS$group, levels = c("IR only", "IR & PTC", "PTC only", "Others"))
2195. bp <- cpm_long_AS %>%
2196. ggplot(aes(x=factor(Timepoint), y=log_cpm)) +
2197. geom_boxplot(aes(fill = group, colour = group),
2198. # outlier.size = 0.0001,
2199. outlier.shape = NA,
2200. size = 0.2,
2201. alpha = 0.35,
2202. # size = 0.5,
2203. linewidth = 0.45,
2204. width = 0.7
2205. ) +
2206. ylim(0, 6.4) +
2207. labs(y = "log2(CPM + 1)") +
2208. stat_pvalue_manual(stat.test.ir,
2209. label = "p.adj.signif",
2210. size = 2.5 ,
2211. bracket.nudge.y = 1.2,
2212. step.group.by = "Timepoint",
2213. step.increase = 0.02,
2214. tip.length = 0.01,
2215. bracket.size = 0.2) +
2216. theme_classic(base_size = 8) +
2217. theme(legend.title = element_blank(),
2218. legend.text = element_text(size = 8),
2219. axis.title.x = element_blank(),
2220. legend.position = "right",
2221. axis.text = element_text(size = 8)) +
2222. scale_x_discrete(labels = c("Avgt00" = "t00", "Avgt04" = "t04", "Avgt30" = "t30"))
2223. bp
2224. #
2225. pdf(paste0("./figures/SUPP_IR_NMD_transcript.pdf"), height = 3.5, width = 4)
2226. print(bp)
2227. dev.off()
2229. svg(paste0("./figures/SUPP_IR_NMD_transcript.svg"), height = 3.5, width = 4)
2230. print(bp)
2231. dev.off()
2235. # APA ---------------------------------------------------------------------
2236. #K APA
2237. K_PPAU_t00_v_t30 <- read.csv("./code/APA/output/APA/PPAU_t00_v_t30.csv")
2238. cluster_t30_vs_t00_df <- read.csv("./code/APA/output/APA/Deseq2_cluster_lrmethod_t00_v_t30.csv")
2239. K_PPAU_t00_v_t30 <- merge(K_PPAU_t00_v_t30, cluster_t30_vs_t00_df[, c("X", "log2FoldChange", "padj")],
2240. by.x = "sequence_cluster", by.y = "X", all.x = T)
2241. K_PPAU_t00_v_t30 <- K_PPAU_t00_v_t30 %>% mutate(Analysis = "Long-read method",
2242. Comparison = "t00 vs t30")
2243. K_PPAU_t00_v_t30 <- K_PPAU_t00_v_t30[,c(1, 21, 22, 24:27)]
2244. colnames(K_PPAU_t00_v_t30)[1] <- "GeneName"
2247. K_PPAU_t00_v_t04 <- read.csv("./code/APA/output/APA/PPAU_t00_v_t04.csv")
2248. cluster_t04_vs_t00_df <- read.csv("./code/APA/output/APA/Deseq2_cluster_lrmethod_t00_v_t04.csv")
2249. K_PPAU_t00_v_t04 <- merge(K_PPAU_t00_v_t04, cluster_t04_vs_t00_df[, c("X", "log2FoldChange", "padj")],
2250. by.x = "sequence_cluster", by.y = "X", all.x = T)
2251. K_PPAU_t00_v_t04 <- K_PPAU_t00_v_t04 %>% mutate(Analysis = "Long-read method",
2252. Comparison = "t00 vs t04")
2253. K_PPAU_t00_v_t04 <- K_PPAU_t00_v_t04[,c(1, 21, 22, 24:27)]
2254. colnames(K_PPAU_t00_v_t04)[1] <- "GeneName"
2258. K_PPAU_t04_v_t30 <- read.csv("./code/APA/output/APA/PPAU_t04_v_t30.csv")
2259. cluster_t30_vs_t04_df <- read.csv("./code/APA/output/APA/Deseq2_cluster_lrmethod_t04_v_t30.csv")
2260. K_PPAU_t04_v_t30 <- merge(K_PPAU_t04_v_t30, cluster_t30_vs_t04_df[, c("X", "log2FoldChange", "padj")],
2261. by.x = "sequence_cluster", by.y = "X", all.x = T)
2262. K_PPAU_t04_v_t30 <- K_PPAU_t04_v_t30 %>% mutate(Analysis = "Long-read method",
2263. Comparison = "t04 vs t30")
2264. K_PPAU_t04_v_t30 <- K_PPAU_t04_v_t30[,c(1, 21, 22, 24:27)]
2265. colnames(K_PPAU_t04_v_t30)[1] <- "GeneName"
2270. #QAPA
2271. QAPA_PPAU_t00_v_30 <- read.csv("./code/APA/output/QAPA/QAPA_PPAU_t00_v_t30.csv")
2272. cluster_t00_vs_t30_df <- read.csv("./code/APA/output/QAPA/Deseq2_cluster_t00_v_t30.csv")
2273. QAPA_PPAU_t00_v_30 <- merge(QAPA_PPAU_t00_v_30, cluster_t00_vs_t30_df[, c("X", "log2FoldChange", "padj")],
2274. by.x = "GeneName", by.y = "X", all.x = T)
2275. QAPA_PPAU_t00_v_30 <- QAPA_PPAU_t00_v_30 %>% mutate(Analysis = "QAPA",
2276. Comparison = "t00 vs t30")
2277. QAPA_PPAU_t00_v_30 <- QAPA_PPAU_t00_v_30[,c(1, 22:27)]
2280. QAPA_PPAU_t00_v_04 <- read.csv("./code/APA/output/QAPA/QAPA_PPAU_t00_v_t04.csv")
2281. cluster_t00_vs_t04_df <- read.csv("./code/APA/output/QAPA/Deseq2_cluster_t00_v_t04.csv")
2282. QAPA_PPAU_t00_v_04 <- merge(QAPA_PPAU_t00_v_04, cluster_t00_vs_t04_df[, c("X", "log2FoldChange", "padj")],
2283. by.x = "GeneName", by.y = "X", all.x = T)
2284. QAPA_PPAU_t00_v_04 <- QAPA_PPAU_t00_v_04 %>% mutate(Analysis = "QAPA",
2285. Comparison = "t00 vs t04")
2286. QAPA_PPAU_t00_v_04 <- QAPA_PPAU_t00_v_04[,c(1, 22:27)]
2289. QAPA_PPAU_t04_v_30 <- read.csv("./code/APA/output/QAPA/QAPA_PPAU_t04_v_t30.csv")
2290. cluster_t04_vs_t30_df <- read.csv("./code/APA/output/QAPA/Deseq2_cluster_t04_v_t30.csv")
2291. QAPA_PPAU_t04_v_30 <- merge(QAPA_PPAU_t04_v_30, cluster_t04_vs_t30_df[, c("X", "log2FoldChange", "padj")],
2292. by.x = "GeneName", by.y = "X", all.x = T)
2293. QAPA_PPAU_t04_v_30 <- QAPA_PPAU_t04_v_30 %>% mutate(Analysis = "QAPA",
2294. Comparison = "t04 vs t30")
2295. QAPA_PPAU_t04_v_30 <- QAPA_PPAU_t04_v_30[,c(1, 22:27)]
2300. df_list <- list(K_PPAU_t00_v_t04, K_PPAU_t00_v_t30, K_PPAU_t04_v_t30,
2301. QAPA_PPAU_t00_v_04, QAPA_PPAU_t00_v_30, QAPA_PPAU_t04_v_30)
2302. df_all <- do.call(rbind, df_list)
2303. df_all <- df_all %>% tidyr::replace_na(list(padj = 1))
2306. library(ggrepel)
2307. library(tidyr)
2309. cluster_plot <- ggplot(df_all, aes(x=dPPAU, y=log2FoldChange)) +
2310. geom_point(alpha = 0.5, size = 1, color = ifelse((df_all$log2FoldChange > 1 & df_all$padj <0.05),
2311. "#FF958C",
2312. ifelse((df_all$log2FoldChange < -1 & df_all$padj <0.05), "#883677", "grey"))) +
2313. # geom_smooth(method = lm, se = FALSE) +
2314. stat_cor(method = "pearson", p.accuracy = 0.001, r.accuracy = 0.001, size = 3.2) +
2315. geom_vline(xintercept = c(-20, 20), linetype = "dashed") +
2316. geom_hline(yintercept = c(-1, 1), linetype = "dashed") +
2318. annotate("text", x = 28, y = -6, label = "Lengthening", hjust = 0, size = 2.8) +
2319. annotate("segment",
2320. x = 30, xend = 85,
2321. y = -7, yend = -7,
2322. arrow = arrow(length = unit(0.2, "cm")),
2323. colour = "red") +
2325. # Arrow + label for "shorter"
2326. annotate("text", x = -87, y = -6, label = "Shortening", hjust = 0, size = 2.8) +
2327. annotate("segment",
2328. x = -30, xend = -85,
2329. y = -7, yend = -7,
2330. arrow = arrow(length = unit(0.2, "cm")),
2331. colour = "blue") +
2333. theme_bw(base_size = 9) +
2334. theme(panel.spacing.x = unit(0.56, "cm")) +
2335. labs(y = "log2FC") +
2336. facet_grid(factor(Analysis) ~ factor(Comparison, c("t00 vs t04", "t04 vs t30", "t00 vs t30")) )
2337. cluster_plot
2339. pdf(paste0("./figures/SUPP_APA.pdf"), height = 5, width = 6.5)
2340. print(cluster_plot)
2341. dev.off()
2343. svg(paste0("./figures/SUPP_APA.svg"), height = 5, width = 6.5)
2344. print(cluster_plot)
2345. dev.off()
2349. # SFARI GENE ENRICHMENT ---------------------------------------------------
2350. # Function to compute Fisher's exact test for pairwise comparisons
2351. #AS term enrichment:
2352. #Run AS code above (at gene-level):
2353. localConseq4$condition <- paste0(localConseq4$condition_1, "vs", localConseq4$condition_2)
2354. localConseq4_summ <- localConseq4 %>% group_by(condition, AStype, gene_id) %>%
2355. summarise(Sum = sum(nrDiff)) %>%
2356. ungroup()
2358. AS <- localConseq4_summ
2360. #All genes (up or down)
2362. #To get all genes:
2363. gene_lists_t00_t04_AS <- AS %>%
2364. filter(condition == "t00vst04") %>%
2365. group_by(AStype) %>%
2366. filter(Sum != 0) %>%
2367. summarise(gene_ids = list(unique(gene_id)), .groups = "drop") %>%
2368. dplyr::select(AStype, gene_ids) %>%
2369. tibble::deframe()
2371. gene_lists_t00_t30_AS <- AS %>%
2372. filter(condition == "t00vst30") %>%
2373. group_by(AStype) %>%
2374. filter(Sum != 0) %>%
2375. summarise(gene_ids = list(unique(gene_id)), .groups = "drop") %>%
2376. dplyr::select(AStype, gene_ids) %>%
2377. tibble::deframe()
2379. gene_lists_t04_t30_AS <- AS %>%
2380. filter(condition == "t04vst30") %>%
2381. group_by(AStype) %>%
2382. filter(Sum != 0) %>%
2383. summarise(gene_ids = list(unique(gene_id)), .groups = "drop") %>%
2384. dplyr::select(AStype, gene_ids) %>%
2385. tibble::deframe()
2387. #Bg genes:
2388. bg.AS <- unique(aSwitchList_part2$isoformFeatures$gene_id) #10,529
2391. SFARI <- read.csv("./data/SFARI-Gene_genes_04-03-2025release_04-15-2025export.csv")
2394. #Testing for enrichment in the each functional consequence:
2395. SFARI.AS <- SFARI %>% filter(gene.symbol %in% bg.AS) %>% pull(gene.symbol) %>% unique()
2398. #Run tests for t00 vs t04:
2399. results_list_t00_t04 <- list()
2401. i = 1
2402. for (i in 1:9) {
2403. setA <- SFARI.AS
2404. setB <- gene_lists_t00_t04_AS[[i]]
2405. result <- compute_fisher_enrich(setA, setB, bg.AS)
2407. # Store the results in the list, and create the pair name
2408. pair_name <- paste0("SFARI vs ", names(gene_lists_t00_t04_AS[i]))
2409. results_list_t00_t04[[pair_name]] <- result
2410. }
2412. results_df_t00_t04 <- do.call(rbind, results_list_t00_t04)
2413. results_df_t00_t04 <- data.frame(pair = names(results_list_t00_t04), results_df_t00_t04)
2414. results_df_t00_t04$p_adjusted <- p.adjust(results_df_t00_t04$p_value, method = "bonferroni")
2416. #Run tests for t00 vs t30:
2417. results_list_t00_t30 <- list()
2419. i = 1
2420. for (i in 1:9) {
2421. setA <- SFARI.AS
2422. setB <- gene_lists_t00_t30_AS[[i]]
2423. result <- compute_fisher_enrich(setA, setB, bg.AS)
2425. # Store the results in the list, and create the pair name
2426. pair_name <- paste0("SFARI vs ", names(gene_lists_t00_t30_AS[i]))
2427. results_list_t00_t30[[pair_name]] <- result
2428. }
2430. results_df_t00_t30 <- do.call(rbind, results_list_t00_t30)
2431. results_df_t00_t30 <- data.frame(pair = names(results_list_t00_t30), results_df_t00_t30)
2432. results_df_t00_t30$p_adjusted <- p.adjust(results_df_t00_t30$p_value, method = "bonferroni")
2434. #Run tests for t00 vs t04:
2435. results_list_t04_t30 <- list()
2437. i = 1
2438. for (i in 1:9) {
2439. setA <- SFARI.AS
2440. setB <- gene_lists_t04_t30_AS[[i]]
2441. result <- compute_fisher_enrich(setA, setB, bg.AS)
2443. # Store the results in the list, and create the pair name
2444. pair_name <- paste0("SFARI vs ", names(gene_lists_t04_t30_AS[i]))
2445. results_list_t04_t30[[pair_name]] <- result
2446. }
2448. results_df_t04_t30 <- do.call(rbind, results_list_t04_t30)
2449. results_df_t04_t30 <- data.frame(pair = names(results_list_t04_t30), results_df_t04_t30)
2450. results_df_t04_t30$p_adjusted <- p.adjust(results_df_t04_t30$p_value, method = "bonferroni")
2452. #Combined data:
2453. results_df_t00_t04$comparison <- "t00 vs t04"
2454. results_df_t00_t30$comparison <- "t00 vs t30"
2455. results_df_t04_t30$comparison <- "t04 vs t30"
2457. results_AS <- rbind(rbind(results_df_t00_t04, results_df_t00_t30), results_df_t04_t30)
2458. results_AS <- results_AS %>%
2459. mutate(
2460. consequence = str_remove(pair, "^SFARI vs\\s+")
2461. )
2462. results_AS
2464. ######Func conseq term enrichment:
2465. bg.func <- bg.AS
2466. #SFARI <- read.csv("./data/SFARI-Gene_genes_04-03-2025release_04-15-2025export.csv")
2469. #Testing for enrichment in the each functional consequence:
2470. SFARI.func <- SFARI.AS
2472. #Genes in each conseq:
2473. Conseq_genes_count <- extractConsequenceEnrichment(
2474. aSwitchList_part2,
2475. alpha=0.05,
2476. dIFcutoff = 0.05,
2477. countGenes = T,
2478. analysisOppositeConsequence=F,
2479. plot=F,
2480. returnResult=TRUE,
2481. returnSummary=F)
2482. #Conseq_genes_count <- Conseq_genes_count %>% filter((condition_1 == "t00" & condition_2 == "t30"))
2483. Conseq_genes_count$comparison <- paste0(Conseq_genes_count$condition_1, " vs ", Conseq_genes_count$condition_2)
2484. Conseq_genes_count <- Conseq_genes_count %>% filter(featureCompared == "3_utr_length" |
2485. featureCompared == "5_utr_length" |
2486. featureCompared == "domain_length" |
2487. featureCompared == "domains_identified" |
2488. featureCompared == "IDR_identified" |
2489. featureCompared == "IDR_length" |
2490. featureCompared == "isoform_length" |
2491. featureCompared == "NMD_status" |
2492. featureCompared == "ORF_length"
2493. # featureCompared == "tss" |
2494. # featureCompared == "tts"
2495. )
2497. Conseq_genes_count <- Conseq_genes_count %>% group_by(comparison, gene_id) %>%
2498. mutate(nrDiff = ifelse((switchConsequence == "3UTR is longer"|
2499. switchConsequence == "5UTR is longer"|
2500. switchConsequence == "Domain gain"|
2501. switchConsequence == "Domain length gain"|
2502. switchConsequence == "IDR gain"|
2503. switchConsequence == "IDR length gain"|
2504. switchConsequence == "Length gain"|
2505. switchConsequence == "NMD sensitive"| ##This is the opposite.
2506. switchConsequence == "ORF is longer"
2507. # switchConsequence == "Tss more upstream"|
2508. # switchConsequence == "Tts more downstream"
2509. ), 1, -1)) %>% ungroup()
2510. info <- Conseq_genes_count %>% dplyr::select(switchConsequence, featureCompared) %>% unique()
2512. Conseq_genes_summ <- Conseq_genes_count %>% group_by(comparison, switchConsequence, gene_id) %>%
2513. dplyr::summarize(sum = sum(nrDiff)) %>%
2514. ungroup()
2515. Conseq_genes_summ <- merge(Conseq_genes_summ, info[,c('switchConsequence', 'featureCompared')])
2518. #To get all genes:
2519. gene_lists_t00_t04 <- Conseq_genes_summ %>%
2520. filter(comparison == "t00 vs t04") %>%
2521. group_by(featureCompared) %>%
2522. filter(sum != 0) %>%
2523. summarise(gene_ids = list(unique(gene_id)), .groups = "drop") %>%
2524. dplyr::select(featureCompared, gene_ids) %>%
2525. tibble::deframe()
2527. gene_lists_t00_t30 <- Conseq_genes_summ %>%
2528. filter(comparison == "t00 vs t30") %>%
2529. group_by(featureCompared) %>%
2530. filter(sum != 0) %>%
2531. summarise(gene_ids = list(unique(gene_id)), .groups = "drop") %>%
2532. dplyr::select(featureCompared, gene_ids) %>%
2533. tibble::deframe()
2535. gene_lists_t04_t30 <- Conseq_genes_summ %>%
2536. filter(comparison == "t04 vs t30") %>%
2537. group_by(featureCompared) %>%
2538. filter(sum != 0) %>%
2539. summarise(gene_ids = list(unique(gene_id)), .groups = "drop") %>%
2540. dplyr::select(featureCompared, gene_ids) %>%
2541. tibble::deframe()
2543. #Run tests for t00 vs t04:
2544. results_list_t00_t04 <- list()
2546. i = 1
2547. for (i in 1:9) {
2548. setA <- SFARI.func
2549. setB <- gene_lists_t00_t04[[i]]
2550. result <- compute_fisher_enrich(setA, setB, bg.func)
2552. # Store the results in the list, and create the pair name
2553. pair_name <- paste0("SFARI vs ", names(gene_lists_t00_t04[i]))
2554. results_list_t00_t04[[pair_name]] <- result
2555. }
2557. results_df_t00_t04 <- do.call(rbind, results_list_t00_t04)
2558. results_df_t00_t04 <- data.frame(pair = names(results_list_t00_t04), results_df_t00_t04)
2559. results_df_t00_t04$p_adjusted <- p.adjust(results_df_t00_t04$p_value, method = "bonferroni")
2561. #Run tests for t00 vs t30:
2562. results_list_t00_t30 <- list()
2564. i = 1
2565. for (i in 1:9) {
2566. setA <- SFARI.func
2567. setB <- gene_lists_t00_t30[[i]]
2568. result <- compute_fisher_enrich(setA, setB, bg.func)
2570. # Store the results in the list, and create the pair name
2571. pair_name <- paste0("SFARI vs ", names(gene_lists_t00_t30[i]))
2572. results_list_t00_t30[[pair_name]] <- result
2573. }
2575. results_df_t00_t30 <- do.call(rbind, results_list_t00_t30)
2576. results_df_t00_t30 <- data.frame(pair = names(results_list_t00_t30), results_df_t00_t30)
2577. results_df_t00_t30$p_adjusted <- p.adjust(results_df_t00_t30$p_value, method = "bonferroni")
2579. #Run tests for t00 vs t04:
2580. results_list_t04_t30 <- list()
2582. i = 1
2583. for (i in 1:9) {
2584. setA <- SFARI.func
2585. setB <- gene_lists_t04_t30[[i]]
2586. result <- compute_fisher_enrich(setA, setB, bg.func)
2588. # Store the results in the list, and create the pair name
2589. pair_name <- paste0("SFARI vs ", names(gene_lists_t04_t30[i]))
2590. results_list_t04_t30[[pair_name]] <- result
2591. }
2593. results_df_t04_t30 <- do.call(rbind, results_list_t04_t30)
2594. results_df_t04_t30 <- data.frame(pair = names(results_list_t04_t30), results_df_t04_t30)
2595. results_df_t04_t30$p_adjusted <- p.adjust(results_df_t04_t30$p_value, method = "bonferroni")
2597. #Combined data:
2598. results_df_t00_t04$comparison <- "t00 vs t04"
2599. results_df_t00_t30$comparison <- "t00 vs t30"
2600. results_df_t04_t30$comparison <- "t04 vs t30"
2602. results_func_conseq <- rbind(rbind(results_df_t00_t04, results_df_t00_t30), results_df_t04_t30)
2603. results_func_conseq <- results_func_conseq %>%
2604. mutate(
2605. consequence = str_remove(pair, "^SFARI vs\\s+")
2606. )
2608. results_AS$conseq_type <- "Alternative Splicing"
2609. results_func_conseq$conseq_type <- "Functional Consequence"
2611. ###Combined results + prep:
2612. combined_enrich <- rbind(results_AS, results_func_conseq)
2614. combined_enrich$consequence <- str_replace(combined_enrich$consequence, "MES",
2615. "Multiple exon splicing")
2616. combined_enrich$consequence <- str_replace(combined_enrich$consequence, "3_utr_length",
2617. "3' UTR length")
2618. combined_enrich$consequence <- str_replace(combined_enrich$consequence, "5_utr_length",
2619. "5' UTR length")
2620. combined_enrich$consequence <- str_replace(combined_enrich$consequence, "domain_length",
2621. "Domain length")
2622. combined_enrich$consequence <- str_replace(combined_enrich$consequence, "domains_identified",
2623. "Domain gain/loss")
2624. combined_enrich$consequence <- str_replace(combined_enrich$consequence, "IDR_identified",
2625. "IDR gain/loss")
2626. combined_enrich$consequence <- str_replace(combined_enrich$consequence, "IDR_length",
2627. "IDR length")
2628. combined_enrich$consequence <- str_replace(combined_enrich$consequence, "isoform_length",
2629. "Isoform length")
2630. combined_enrich$consequence <- str_replace(combined_enrich$consequence, "NMD_status",
2631. "NMD status")
2632. combined_enrich$consequence <- str_replace(combined_enrich$consequence, "ORF_length",
2633. "ORF length")
2635. combined_enrich$consequence <- str_replace(combined_enrich$consequence, "A3",
2636. "Alternative 3' splice site")
2637. combined_enrich$consequence <- str_replace(combined_enrich$consequence, "A5",
2638. "Alternative 5' splice site")
2639. combined_enrich$consequence <- str_replace(combined_enrich$consequence, "ES",
2640. "Exon splicing")
2641. combined_enrich$consequence <- str_replace(combined_enrich$consequence, "IR",
2642. "Intron splicing")
2644. combined_enrich$consequence <- str_replace(combined_enrich$consequence, "MIC",
2645. "Microexon splicing")
2649. combined_enrich <- combined_enrich %>%
2650. mutate(consequence = fct_relevel(consequence,
2651. "Alternative 3' splice site",
2652. "Alternative 5' splice site",
2653. "Exon splicing",
2654. "Microexon splicing",
2655. "Intron splicing",
2656. "Multiple exon splicing",
2657. "IDR length",
2658. "IDR gain/loss",
2659. "Domain length",
2660. "Domain gain/loss",
2661. "NMD status",
2662. "ORF length",
2663. "Isoform length",
2664. "3' UTR length",
2665. "5' UTR length"))
2667. combined_enrich$sig <- ifelse(combined_enrich$p_adjusted <= 0.05, TRUE, FALSE)
2668. combined_enrich$comparison <- factor(combined_enrich$comparison, levels = c("t00 vs t04", "t04 vs t30", "t00 vs t30"))
2671. ##FIg:
2672. combined_enrich_filter <- combined_enrich %>% filter(consequence != "ATTS" &
2673. consequence != "ATSS" &
2674. consequence != "MEE" &
2675. consequence != "Isoform length" &
2676. consequence != "Domain length" &
2677. consequence != "IDR length" &
2678. consequence != "Alternative 3' splice site" &
2679. consequence != "Alternative 5' splice site" &
2680. consequence != "Multiple exon splicing" &
2681. consequence != "ORF length")
2683. desired_order <- c("Exon splicing",
2684. "Microexon splicing",
2685. "Intron splicing",
2686. "NMD status",
2687. "Domain gain/loss",
2688. "IDR gain/loss",
2689. "5' UTR length",
2690. "3' UTR length")
2691. combined_enrich_filter <- combined_enrich_filter %>%
2692. mutate(consequence = factor(consequence, levels = rev(desired_order)))
2693. #
2695. ##Plot:
2696. ASD_bar <- combined_enrich_filter %>%
2697. ggplot(aes(x = odds_ratio.odds.ratio,
2698. y = consequence,
2699. fill = overlap)) +
2700. geom_col(width = 0.6) +
2701. facet_wrap(~ factor(comparison, c("t00 vs t04", "t04 vs t30", "t00 vs t30")),
2702. scales = "free_x") +
2703. scale_fill_gradient(name = "Gene Count", low = "lightblue", high = "darkblue") +
2704. scale_x_continuous(limits = c(0, 2.85),
2705. expand = expansion(mult = c(0, .1))) +
2706. labs(title = "SFARI Gene Enrichment",
2707. x = "Odds Ratio",
2708. y = NULL) +
2709. geom_vline(xintercept = 1, linetype = "dashed", colour = "black") +
2710. geom_text(data = subset(combined_enrich_filter, sig == TRUE),
2711. aes(x = odds_ratio.odds.ratio,
2712. y = consequence,
2713. label = "*"),
2714. hjust = -0.4,
2715. vjust = 0.75,
2716. colour = "black",
2717. size = 3.5) +
2718. theme_bw(base_size = 8) +
2719. theme(
2720. plot.title = element_text(hjust = 0.5,
2721. size = 8, face = "bold"),
2722. strip.text = element_text(size = 8),
2723. axis.text = element_text(colour = "black", size = 6.4),
2724. axis.title.y = element_blank(),
2725. legend.position = "none"
2726. )
2727. ASD_bar
2729. legend_ASD_bar <- ggpubr::get_legend(ASD_bar +
2730. theme( legend.position = "top",
2731. legend.box = "horizontal",
2733. legend.text = element_text(size = 7),
2734. legend.key.size = unit(0.25, 'cm'),
2736. legend.margin = margin(
2737. t = 0,
2738. r = 0,
2739. b = 0,
2740. l = 0,
2741. unit = "cm"
2742. ) )+
2743. guides(fill = guide_colourbar(title.vjust = 0.8,
2744. barwidth = 3.5))
2745. )
2747. combined_ASD_bar <- ggdraw() +
2748. draw_plot(ASD_bar + theme(plot.title = element_text(vjust = 16),
2749. axis.title.x = element_text(vjust = -6))
2750. , x = 0, y = 0, width = 1, height = 0.934) +
2751. draw_plot(legend_ASD_bar, x = 0, y = 0.447, width = 1, height = 1)
2752. combined_ASD_bar
2755. #
2757. # GO TERMS - SWITCHING GENES ----------------------------------------------
2758. library(gprofiler2)
2759. run_GO_enrichment <- function(gene_list, bg, searches) {
2760. go_result <- gost(gene_list,
2761. organism = "hsapiens",
2762. significant = T,
2763. exclude_iea = F,
2764. user_threshold = 0.05,
2765. correction_method = "fdr",
2766. custom_bg = bg,
2767. sources = searches)
2768. return(go_result)
2769. }
2771. #Run AS code above (at gene-level):
2772. localConseq4$condition <- paste0(localConseq4$condition_1, "vs", localConseq4$condition_2)
2773. localConseq4_summ <- localConseq4 %>% group_by(condition, AStype, gene_id) %>%
2774. summarise(Sum = sum(nrDiff)) %>%
2775. ungroup()
2777. ###T00 vs T30
2778. AS_t00_t30 <- localConseq4_summ %>% filter(condition == "t00vst30")
2780. #All genes (up or down)
2781. genes_ATTS <- AS_t00_t30 %>% filter(AStype == "ATTS" & Sum != 0) %>% pull(gene_id)
2782. genes_MIC <- AS_t00_t30 %>% filter(AStype == "MIC" & Sum != 0) %>% pull(gene_id)
2783. genes_A3 <- AS_t00_t30 %>% filter(AStype == "A3" & Sum != 0) %>% pull(gene_id)
2784. genes_A5 <- AS_t00_t30 %>% filter(AStype == "A5" & Sum != 0) %>% pull(gene_id)
2785. genes_ATSS <- AS_t00_t30 %>% filter(AStype == "ATSS" & Sum != 0) %>% pull(gene_id)
2786. genes_ES <- AS_t00_t30 %>% filter(AStype == "ES" & Sum != 0) %>% pull(gene_id)
2787. genes_MES <- AS_t00_t30 %>% filter(AStype == "MES" & Sum != 0) %>% pull(gene_id)
2788. genes_IR <- AS_t00_t30 %>% filter(AStype == "IR" & Sum != 0) %>% pull(gene_id)
2790. #Bg genes:
2791. all_AS_genes <- unique(aSwitchList_part2$isoformFeatures$gene_id) #10529
2794. #GO terms for all AS results
2795. searches <- c("GO:BP", "GO:MF", "GO:CC", "MIRNA")
2796. gene_lists <- list(genes_ATTS = genes_ATTS,
2797. genes_MIC = genes_MIC,
2798. genes_A3 = genes_A3,
2799. genes_A5 = genes_A5,
2800. genes_ATSS = genes_ATSS,
2801. genes_ES = genes_ES,
2802. genes_MES = genes_MES,
2803. genes_IR = genes_IR)
2805. # Initialize an empty list to store results
2806. enrichment_results <- list()
2807. i = 1
2808. # Loop through each gene list and perform GO enrichment
2809. for (i in seq_along(gene_lists)) {
2810. enrichment_results[[names(gene_lists)[i]]] <- run_GO_enrichment(gene_lists[[i]],
2811. bg = all_AS_genes,
2812. searches = searches)
2813. }
2814. #
2816. go_results_IR <- enrichment_results$genes_IR$result
2817. go_results_IR$query <- "IR GO Terms - t00 vs t30"
2818. # RP_IR <- gene_lists[["genes_IR"]][grepl("^RP", gene_lists[["genes_IR"]])]
2819. go_results_IR <- go_results_IR %>%
2820. group_by(source) %>%
2821. arrange(p_value) %>%
2822. dplyr::slice(1:5) %>%
2823. ungroup()
2824. go_results_IR$term_name_wrapped <- str_wrap(go_results_IR$term_name, width = 20)
2825. # Create a dot plot for each gene list and GO term
2826. IR_go <- ggplot(go_results_IR, aes(x = reorder(term_name_wrapped, -p_value), y = -log10(p_value), fill = intersection_size)) +
2827. geom_col(width = c(0.6*1, 0.6*1.1, 0.6*1.1, 0.6*1.1, 0.6*1.1, 0.6*1.1, 0.6*1)) +
2828. scale_y_continuous(expand = expansion(mult = c(0, .1))) +
2829. labs(title = "IR GO Terms - t00 vs t30",
2830. y = expression(paste(-log[10], "(adjusted ", italic("P"), ")")) ) +
2831. theme_bw(base_size = 8) +
2832. coord_flip() +
2833. scale_fill_gradient(name = "Gene Count", low = "lightblue", high = "darkblue") +
2834. theme(plot.title = element_text(hjust = 0.5, size = 8, colour = "black", face = "bold"),
2835. axis.title.y = element_blank(),
2836. axis.text = element_text(size = 6.4, colour = "black"),
2837. axis.title = element_text(size = 8, colour = "black"),
2838. legend.position = "bottom",
2839. strip.text = element_text(size = 6.4, colour = "black"),
2840. # plot.margin = margin(t = 6, b = 0, l = 16, r = 6, unit = "pt"),
2841. legend.margin = margin(t = -5, b = 0, l = -25, r = 0, unit = "pt"),
2842. legend.key.size = unit(0.25, 'cm')
2843. ) +
2844. facet_wrap(~source,
2845. scales = "free_y",
2846. dir = "v",
2847. strip.position = "right",
2848. ncol = 1) + # Facet by source if applicable
2849. force_panelsizes(row = c(1.25, 4.5, 1.25))
2850. IR_go
2855. # EXPORT FIG 4 ------------------------------------------------------------
2856. library(Cairo)
2858. Figure_4 <- ggdraw() +
2859. draw_plot(bars_by_tr, x = 0, y = 0.78, width = 0.5, height = 0.22) +
2860. draw_plot(final_plot, x = 0.5, y = 0.78, width = 0.5, height = 0.22) +
2861. draw_plot(testp, x = 0, y = 0.35, width = 0.5, height = 0.43) +
2862. draw_plot(combined_ASD_bar, x = 0.5, y = 0.38, width = 0.5, height = 0.3815) +
2863. draw_plot(IR_go, x = 0.01, y = 0, width = 0.35, height = 0.35) +
2864. draw_plot(AGO1, x = 0.35, y = 0, width = 0.65, height = 0.35) +
2866. draw_plot_label(label = c("A", "B", "C", "D", "E", "F"),
2867. size = 12,
2868. x = c(0, 0.5, 0, 0.5, 0, 0.355),
2869. y = c(1, 1, 0.78, 0.78, 0.35, 0.35))
2870. Figure_4
2871. #
2873. CairoPDF(paste0("./figures/Fig_4.pdf"), height = 8, width = 6.5)
2874. print(Figure_4)
2875. dev.off()
2878. svg(paste0("./figures/Fig_4.svg"), height = 8, width = 6.5)
2879. print(Figure_4)
2880. dev.off()
2886. # SUPP. FIGS & TABLES ------------------------------------------------------
2889. # GO TERMS - SWITCHING GENES ----------------------------------------------
2890. library(gprofiler2)
2891. run_GO_enrichment <- function(gene_list, bg, searches) {
2892. # Convert gene symbols to Entrez IDs (if necessary)
2893. # gene_ids <- bitr(gene_list, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = OrgDb)
2894. go_result <- gost(gene_list,
2895. organism = "hsapiens",
2896. significant = T,
2897. exclude_iea = F,
2898. user_threshold = 0.05,
2899. correction_method = "fdr",
2900. custom_bg = bg,
2901. sources = searches)
2902. return(go_result)
2903. }
2906. #Run AS code above (at gene-level):
2907. localConseq4$condition <- paste0(localConseq4$condition_1, "vs", localConseq4$condition_2)
2908. localConseq4_summ <- localConseq4 %>% group_by(condition, AStype, gene_id) %>%
2909. summarise(Sum = sum(nrDiff)) %>%
2910. ungroup()
2912. ###T00 vs T30
2913. AS_t00_t30 <- localConseq4_summ %>% filter(condition == "t00vst30")
2915. #All genes (up or down)
2916. genes_ATTS <- AS_t00_t30 %>% filter(AStype == "ATTS" & Sum != 0) %>% pull(gene_id)
2917. genes_MIC <- AS_t00_t30 %>% filter(AStype == "MIC" & Sum != 0) %>% pull(gene_id)
2918. genes_A3 <- AS_t00_t30 %>% filter(AStype == "A3" & Sum != 0) %>% pull(gene_id)
2919. genes_A5 <- AS_t00_t30 %>% filter(AStype == "A5" & Sum != 0) %>% pull(gene_id)
2920. genes_ATSS <- AS_t00_t30 %>% filter(AStype == "ATSS" & Sum != 0) %>% pull(gene_id)
2921. genes_ES <- AS_t00_t30 %>% filter(AStype == "ES" & Sum != 0) %>% pull(gene_id)
2922. genes_MES <- AS_t00_t30 %>% filter(AStype == "MES" & Sum != 0) %>% pull(gene_id)
2923. genes_IR <- AS_t00_t30 %>% filter(AStype == "IR" & Sum != 0) %>% pull(gene_id)
2925. #Bg genes:
2926. all_AS_genes <- unique(aSwitchList_part2$isoformFeatures$gene_id) #10529
2929. #GO terms for all AS results
2930. searches <- c("GO:BP", "GO:MF", "GO:CC", "MIRNA")
2931. gene_lists <- list(genes_ATTS = genes_ATTS,
2932. genes_MIC = genes_MIC,
2933. genes_A3 = genes_A3,
2934. genes_A5 = genes_A5,
2935. genes_ATSS = genes_ATSS,
2936. genes_ES = genes_ES,
2937. genes_MES = genes_MES,
2938. genes_IR = genes_IR)
2940. # Initialize an empty list to store results
2941. enrichment_results <- list()
2942. i = 1
2943. # Loop through each gene list and perform GO enrichment
2944. for (i in seq_along(gene_lists)) {
2945. enrichment_results[[names(gene_lists)[i]]] <- run_GO_enrichment(gene_lists[[i]],
2946. bg = all_AS_genes,
2947. searches = searches)
2948. }
2949. #
2951. ###T00 vs T04
2952. AS_t00_t04 <- localConseq4_summ %>% filter(condition == "t00vst04")
2954. #All genes (up or down)
2955. genes_ATTS <- AS_t00_t04 %>% filter(AStype == "ATTS" & Sum != 0) %>% pull(gene_id)
2956. genes_MIC <- AS_t00_t04 %>% filter(AStype == "MIC" & Sum != 0) %>% pull(gene_id)
2957. genes_A3 <- AS_t00_t04 %>% filter(AStype == "A3" & Sum != 0) %>% pull(gene_id)
2958. genes_A5 <- AS_t00_t04 %>% filter(AStype == "A5" & Sum != 0) %>% pull(gene_id)
2959. genes_ATSS <- AS_t00_t04 %>% filter(AStype == "ATSS" & Sum != 0) %>% pull(gene_id)
2960. genes_ES <- AS_t00_t04 %>% filter(AStype == "ES" & Sum != 0) %>% pull(gene_id)
2961. genes_MES <- AS_t00_t04 %>% filter(AStype == "MES" & Sum != 0) %>% pull(gene_id)
2962. genes_IR <- AS_t00_t04 %>% filter(AStype == "IR" & Sum != 0) %>% pull(gene_id)
2964. #Bg genes:
2965. all_AS_genes
2967. #GO terms for all AS results
2968. searches <- c("GO:BP", "GO:MF", "GO:CC", "MIRNA")
2969. gene_lists <- list(genes_ATTS = genes_ATTS,
2970. genes_MIC = genes_MIC,
2971. genes_A3 = genes_A3,
2972. genes_A5 = genes_A5,
2973. genes_ATSS = genes_ATSS,
2974. genes_ES = genes_ES,
2975. genes_MES = genes_MES,
2976. genes_IR = genes_IR)
2978. # Initialize an empty list to store results
2979. enrichment_results <- list()
2980. i = 1
2981. # Loop through each gene list and perform GO enrichment
2982. for (i in seq_along(gene_lists)) {
2983. enrichment_results[[names(gene_lists)[i]]] <- run_GO_enrichment(gene_lists[[i]],
2984. bg = all_AS_genes,
2985. searches = searches)
2986. }
2989. ###T04 vs T30
2990. AS_t04_t30 <- localConseq4_summ %>% filter(condition == "t04vst30")
2992. #All genes (up or down)
2993. genes_ATTS <- AS_t04_t30 %>% filter(AStype == "ATTS" & Sum != 0) %>% pull(gene_id)
2994. genes_MIC <- AS_t04_t30 %>% filter(AStype == "MIC" & Sum != 0) %>% pull(gene_id)
2995. genes_A3 <- AS_t04_t30 %>% filter(AStype == "A3" & Sum != 0) %>% pull(gene_id)
2996. genes_A5 <- AS_t04_t30 %>% filter(AStype == "A5" & Sum != 0) %>% pull(gene_id)
2997. genes_ATSS <- AS_t04_t30 %>% filter(AStype == "ATSS" & Sum != 0) %>% pull(gene_id)
2998. genes_ES <- AS_t04_t30 %>% filter(AStype == "ES" & Sum != 0) %>% pull(gene_id)
2999. genes_MES <- AS_t04_t30 %>% filter(AStype == "MES" & Sum != 0) %>% pull(gene_id)
3000. genes_IR <- AS_t04_t30 %>% filter(AStype == "IR" & Sum != 0) %>% pull(gene_id)
3002. #Bg genes:
3003. all_AS_genes
3005. #GO terms for all AS results
3006. searches <- c("GO:BP", "GO:MF", "GO:CC", "MIRNA")
3007. gene_lists <- list(genes_ATTS = genes_ATTS,
3008. genes_MIC = genes_MIC,
3009. genes_A3 = genes_A3,
3010. genes_A5 = genes_A5,
3011. genes_ATSS = genes_ATSS,
3012. genes_ES = genes_ES,
3013. genes_MES = genes_MES,
3014. genes_IR = genes_IR)
3016. # Initialize an empty list to store results
3017. enrichment_results <- list()
3018. i = 1
3019. # Loop through each gene list and perform GO enrichment
3020. for (i in seq_along(gene_lists)) {
3021. enrichment_results[[names(gene_lists)[i]]] <- run_GO_enrichment(gene_lists[[i]],
3022. bg = all_AS_genes,
3023. searches = searches)
3024. }
3028. # GO TERMS - FUNC CONSEQ --------------------------------------------------
3029. ###Do for func. conseq:
3030. Conseq_genes <- extractConsequenceEnrichment(
3031. aSwitchList_part2,
3032. alpha=0.05,
3033. dIFcutoff = 0.05,
3034. countGenes = T,
3035. analysisOppositeConsequence=F,
3036. plot=F,
3037. returnResult=TRUE,
3038. returnSummary=T)
3041. ##T00 vs T30:
3042. Conseq_genes_count <- extractConsequenceEnrichment(
3043. aSwitchList_part2,
3044. alpha=0.05,
3045. dIFcutoff = 0.05,
3046. countGenes = T,
3047. analysisOppositeConsequence=F,
3048. plot=F,
3049. returnResult=TRUE,
3050. returnSummary=F)
3051. Conseq_genes_count <- Conseq_genes_count %>% filter((condition_1 == "t00" & condition_2 == "t30"))
3053. Conseq_genes_count$nrDiff <- ifelse((Conseq_genes_count$switchConsequence == "3UTR is longer"|
3054. Conseq_genes_count$switchConsequence == "5UTR is longer"|
3055. Conseq_genes_count$switchConsequence == "Domain gain"|
3056. Conseq_genes_count$switchConsequence == "Domain length gain"|
3057. Conseq_genes_count$switchConsequence == "IDR gain"|
3058. Conseq_genes_count$switchConsequence == "IDR length gain"|
3059. Conseq_genes_count$switchConsequence == "Length gain"|
3060. Conseq_genes_count$switchConsequence == "NMD sensitive"| ##This is the opposite.
3061. Conseq_genes_count$switchConsequence == "ORF is longer"|
3062. Conseq_genes_count$switchConsequence == "Tss more upstream"|
3063. Conseq_genes_count$switchConsequence == "Tts more downstream"), 1, -1)
3065. Conseq_genes_summ <- Conseq_genes_count %>%
3066. group_by(switchConsequence, gene_id) %>%
3067. summarize(sum = sum(nrDiff)) %>%
3068. ungroup()
3069. summary_genes <- Conseq_genes_summ %>%
3070. group_by(switchConsequence) %>%
3071. filter(sum != 0)
3072. #
3074. #I can just take the genes that match the ones above (b/c all others are 0)
3075. genes_3UTR <- Conseq_genes_summ %>% filter(switchConsequence == "3UTR is longer" |
3076. switchConsequence == "3UTR is shorter") %>% pull(gene_id)
3077. genes_5UTR <- Conseq_genes_summ %>% filter(switchConsequence == "5UTR is longer" |
3078. switchConsequence == "5UTR is shorter") %>% pull(gene_id)
3079. genes_domain_number <- Conseq_genes_summ %>% filter(switchConsequence == "Domain gain" |
3080. switchConsequence == "Domain loss") %>% pull(gene_id)
3081. genes_domain_length <- Conseq_genes_summ %>% filter(switchConsequence == "Domain length gain" |
3082. switchConsequence == "Domain length loss") %>% pull(gene_id)
3083. genes_IDR_number <- Conseq_genes_summ %>% filter(switchConsequence == "IDR gain" |
3084. switchConsequence == "IDR loss") %>% pull(gene_id)
3085. genes_IDR_length <- Conseq_genes_summ %>% filter(switchConsequence == "IDR length gain" |
3086. switchConsequence == "IDR length loss") %>% pull(gene_id)
3087. genes_length <- Conseq_genes_summ %>% filter(switchConsequence == "Length gain" |
3088. switchConsequence == "Length loss") %>% pull(gene_id)
3089. genes_NMD <- Conseq_genes_summ %>% filter(switchConsequence == "NMD sensitive" |
3090. switchConsequence == "NMD insensitive") %>% pull(gene_id)
3091. genes_ORF_length <- Conseq_genes_summ %>% filter(switchConsequence == "ORF is longer" |
3092. switchConsequence == "ORF is shorter") %>% pull(gene_id)
3093. genes_TSS_change <- Conseq_genes_summ %>% filter(switchConsequence == "Tss more upstream" |
3094. switchConsequence == "Tss more downstream") %>% pull(gene_id)
3095. genes_TTS_change <- Conseq_genes_summ %>% filter(switchConsequence == "Tts more downstream" |
3096. switchConsequence == "Tts more upstream") %>% pull(gene_id)
3099. # Your gene lists in a named list
3100. gene_lists <- list(
3101. genes_3UTR = genes_3UTR,
3102. genes_5UTR = genes_5UTR,
3103. genes_domain_number = genes_domain_number,
3104. genes_domain_length = genes_domain_length,
3105. genes_IDR_number = genes_IDR_number,
3106. genes_IDR_length = genes_IDR_length,
3107. genes_length = genes_length,
3108. genes_NMD = genes_NMD,
3109. genes_ORF_length = genes_ORF_length,
3110. genes_TSS_change = genes_TSS_change,
3111. genes_TTS_change = genes_TTS_change)
3114. #GO terms for all AS results
3115. searches <- c("GO:BP", "GO:MF", "GO:CC", "MIRNA")
3116. bg_genes <- unique(aSwitchList_part2$isoformFeatures$gene_id)
3117. # Initialize an empty list to store results
3118. enrichment_results <- list()
3119. i = 1
3120. # Loop through each gene list and perform GO enrichment
3121. for (i in seq_along(gene_lists)) {
3122. enrichment_results[[names(gene_lists)[i]]] <- run_GO_enrichment(gene_lists[[i]],
3123. bg = bg_genes,
3124. searches = searches)
3125. }
3126. #go_results_IDR_number <- enrichment_results$genes_IDR_number$result
3127. t00_t30 <- enrichment_results
3131. ##T00 vs T04:
3132. Conseq_genes_count <- extractConsequenceEnrichment(
3133. aSwitchList_part2,
3134. alpha=0.05,
3135. dIFcutoff = 0.05,
3136. countGenes = T,
3137. analysisOppositeConsequence=F,
3138. plot=F,
3139. returnResult=TRUE,
3140. returnSummary=F)
3141. Conseq_genes_count <- Conseq_genes_count %>% filter((condition_1 == "t00" & condition_2 == "t04"))
3143. Conseq_genes_count$nrDiff <- ifelse((Conseq_genes_count$switchConsequence == "3UTR is longer"|
3144. Conseq_genes_count$switchConsequence == "5UTR is longer"|
3145. Conseq_genes_count$switchConsequence == "Domain gain"|
3146. Conseq_genes_count$switchConsequence == "Domain length gain"|
3147. Conseq_genes_count$switchConsequence == "IDR gain"|
3148. Conseq_genes_count$switchConsequence == "IDR length gain"|
3149. Conseq_genes_count$switchConsequence == "Length gain"|
3150. Conseq_genes_count$switchConsequence == "NMD sensitive"| ##This is the opposite.
3151. Conseq_genes_count$switchConsequence == "ORF is longer"|
3152. Conseq_genes_count$switchConsequence == "Tss more upstream"|
3153. Conseq_genes_count$switchConsequence == "Tts more downstream"), 1, -1)
3155. Conseq_genes_summ <- Conseq_genes_count %>%
3156. group_by(switchConsequence, gene_id) %>%
3157. summarize(sum = sum(nrDiff)) %>%
3158. ungroup()
3159. summary_genes <- Conseq_genes_summ %>%
3160. group_by(switchConsequence) %>%
3161. filter(sum != 0)
3163. #I can just take the genes that match the ones above (b/c all others are 0)
3164. genes_3UTR <- Conseq_genes_summ %>% filter(switchConsequence == "3UTR is longer" |
3165. switchConsequence == "3UTR is shorter") %>% pull(gene_id)
3166. genes_5UTR <- Conseq_genes_summ %>% filter(switchConsequence == "5UTR is longer" |
3167. switchConsequence == "5UTR is shorter") %>% pull(gene_id)
3168. genes_domain_number <- Conseq_genes_summ %>% filter(switchConsequence == "Domain gain" |
3169. switchConsequence == "Domain loss") %>% pull(gene_id)
3170. genes_domain_length <- Conseq_genes_summ %>% filter(switchConsequence == "Domain length gain" |
3171. switchConsequence == "Domain length loss") %>% pull(gene_id)
3172. genes_IDR_number <- Conseq_genes_summ %>% filter(switchConsequence == "IDR gain" |
3173. switchConsequence == "IDR loss") %>% pull(gene_id)
3174. genes_IDR_length <- Conseq_genes_summ %>% filter(switchConsequence == "IDR length gain" |
3175. switchConsequence == "IDR length loss") %>% pull(gene_id)
3176. genes_length <- Conseq_genes_summ %>% filter(switchConsequence == "Length gain" |
3177. switchConsequence == "Length loss") %>% pull(gene_id)
3178. genes_NMD <- Conseq_genes_summ %>% filter(switchConsequence == "NMD sensitive" |
3179. switchConsequence == "NMD insensitive") %>% pull(gene_id)
3180. genes_ORF_length <- Conseq_genes_summ %>% filter(switchConsequence == "ORF is longer" |
3181. switchConsequence == "ORF is shorter") %>% pull(gene_id)
3182. genes_TSS_change <- Conseq_genes_summ %>% filter(switchConsequence == "Tss more upstream" |
3183. switchConsequence == "Tss more downstream") %>% pull(gene_id)
3184. genes_TTS_change <- Conseq_genes_summ %>% filter(switchConsequence == "Tts more downstream" |
3185. switchConsequence == "Tts more upstream") %>% pull(gene_id)
3188. # Your gene lists in a named list
3189. gene_lists <- list(
3190. genes_3UTR = genes_3UTR,
3191. genes_5UTR = genes_5UTR,
3192. genes_domain_number = genes_domain_number,
3193. genes_domain_length = genes_domain_length,
3194. genes_IDR_number = genes_IDR_number,
3195. genes_IDR_length = genes_IDR_length,
3196. genes_length = genes_length,
3197. genes_NMD = genes_NMD,
3198. genes_ORF_length = genes_ORF_length,
3199. genes_TSS_change = genes_TSS_change,
3200. genes_TTS_change = genes_TTS_change)
3203. #GO terms for all AS results
3204. searches <- c("GO:BP", "GO:MF", "GO:CC", "MIRNA")
3205. bg_genes
3206. # Initialize an empty list to store results
3207. enrichment_results <- list()
3208. i = 1
3209. # Loop through each gene list and perform GO enrichment
3210. for (i in seq_along(gene_lists)) {
3211. enrichment_results[[names(gene_lists)[i]]] <- run_GO_enrichment(gene_lists[[i]],
3212. bg = bg_genes,
3213. searches = searches)
3214. }
3215. t00_t04 <- enrichment_results
3219. ##T04 vs T30:
3220. Conseq_genes_count <- extractConsequenceEnrichment(
3221. aSwitchList_part2,
3222. alpha=0.05,
3223. dIFcutoff = 0.05,
3224. countGenes = T,
3225. analysisOppositeConsequence=F,
3226. plot=F,
3227. returnResult=TRUE,
3228. returnSummary=F)
3229. Conseq_genes_count <- Conseq_genes_count %>% filter((condition_1 == "t04" & condition_2 == "t30"))
3231. Conseq_genes_count$nrDiff <- ifelse((Conseq_genes_count$switchConsequence == "3UTR is longer"|
3232. Conseq_genes_count$switchConsequence == "5UTR is longer"|
3233. Conseq_genes_count$switchConsequence == "Domain gain"|
3234. Conseq_genes_count$switchConsequence == "Domain length gain"|
3235. Conseq_genes_count$switchConsequence == "IDR gain"|
3236. Conseq_genes_count$switchConsequence == "IDR length gain"|
3237. Conseq_genes_count$switchConsequence == "Length gain"|
3238. Conseq_genes_count$switchConsequence == "NMD sensitive"| ##This is the opposite.
3239. Conseq_genes_count$switchConsequence == "ORF is longer"|
3240. Conseq_genes_count$switchConsequence == "Tss more upstream"|
3241. Conseq_genes_count$switchConsequence == "Tts more downstream"), 1, -1)
3243. Conseq_genes_summ <- Conseq_genes_count %>%
3244. group_by(switchConsequence, gene_id) %>%
3245. summarize(sum = sum(nrDiff)) %>%
3246. ungroup()
3247. summary_genes <- Conseq_genes_summ %>%
3248. group_by(switchConsequence) %>%
3249. filter(sum != 0)
3252. #I can just take the genes that match the ones above (b/c all others are 0)
3253. genes_3UTR <- Conseq_genes_summ %>% filter(switchConsequence == "3UTR is longer" |
3254. switchConsequence == "3UTR is shorter") %>% pull(gene_id)
3255. genes_5UTR <- Conseq_genes_summ %>% filter(switchConsequence == "5UTR is longer" |
3256. switchConsequence == "5UTR is shorter") %>% pull(gene_id)
3257. genes_domain_number <- Conseq_genes_summ %>% filter(switchConsequence == "Domain gain" |
3258. switchConsequence == "Domain loss") %>% pull(gene_id)
3259. genes_domain_length <- Conseq_genes_summ %>% filter(switchConsequence == "Domain length gain" |
3260. switchConsequence == "Domain length loss") %>% pull(gene_id)
3261. genes_IDR_number <- Conseq_genes_summ %>% filter(switchConsequence == "IDR gain" |
3262. switchConsequence == "IDR loss") %>% pull(gene_id)
3263. genes_IDR_length <- Conseq_genes_summ %>% filter(switchConsequence == "IDR length gain" |
3264. switchConsequence == "IDR length loss") %>% pull(gene_id)
3265. genes_length <- Conseq_genes_summ %>% filter(switchConsequence == "Length gain" |
3266. switchConsequence == "Length loss") %>% pull(gene_id)
3267. genes_NMD <- Conseq_genes_summ %>% filter(switchConsequence == "NMD sensitive" |
3268. switchConsequence == "NMD insensitive") %>% pull(gene_id)
3269. genes_ORF_length <- Conseq_genes_summ %>% filter(switchConsequence == "ORF is longer" |
3270. switchConsequence == "ORF is shorter") %>% pull(gene_id)
3271. genes_TSS_change <- Conseq_genes_summ %>% filter(switchConsequence == "Tss more upstream" |
3272. switchConsequence == "Tss more downstream") %>% pull(gene_id)
3273. genes_TTS_change <- Conseq_genes_summ %>% filter(switchConsequence == "Tts more downstream" |
3274. switchConsequence == "Tts more upstream") %>% pull(gene_id)
3277. # Your gene lists in a named list
3278. gene_lists <- list(
3279. genes_3UTR = genes_3UTR,
3280. genes_5UTR = genes_5UTR,
3281. genes_domain_number = genes_domain_number,
3282. genes_domain_length = genes_domain_length,
3283. genes_IDR_number = genes_IDR_number,
3284. genes_IDR_length = genes_IDR_length,
3285. genes_length = genes_length,
3286. genes_NMD = genes_NMD,
3287. genes_ORF_length = genes_ORF_length,
3288. genes_TSS_change = genes_TSS_change,
3289. genes_TTS_change = genes_TTS_change)
3292. #GO terms for all AS results
3293. searches <- c("GO:BP", "GO:MF", "GO:CC", "MIRNA")
3294. bg_genes
3295. # Initialize an empty list to store results
3296. enrichment_results <- list()
3297. i = 1
3298. # Loop through each gene list and perform GO enrichment
3299. for (i in seq_along(gene_lists)) {
3300. enrichment_results[[names(gene_lists)[i]]] <- run_GO_enrichment(gene_lists[[i]],
3301. bg = bg_genes,
3302. searches = searches)
3303. }
3304. t04_t30 <- enrichment_results

Figures.R at commit 080e8a9, under MIT · at the source