CompExonSpace/PCA_GA(FigS9).Rmd at main · Novartis/CompExonSpace · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95

```{r}
#load packages
library(dplyr)
library(ggplot2)
library(ggrepel)
library(stringr)

```

```{r}

# function to load pca data (one combined file per experiment)
run_ga_pca <- function(exp_id,
                       data_dir = Sys.getenv("MYDATA"),
                       top_n = 1000) {

  message("Running GA‑filtered PCA for ", exp_id)

  # Load merged PSI file for experiment 1, 2, 3, ..
  filepath <- file.path(data_dir, "psi_persample", paste0(exp_id, ".csv"))
  df <- read.csv(filepath, stringsAsFactors = FALSE, check.names = FALSE)
  rownames(df) <- make.unique(df$exonID)

  # Filter for Exons with GA-ending
  df <- df |> dplyr::select(-exonID) |> dplyr::filter(Exon2 == "GA")

  experiment_cols <- grep(paste0("^", exp_id, "_"), colnames(df), value = TRUE)
  psidata <- df[, experiment_cols, drop = FALSE]

  # Variance + top exons
  variances <- apply(psidata, 1, var, na.rm = TRUE)
  variances <- variances[!is.na(variances)]

  top <- names(sort(variances, decreasing = TRUE))[1:min(top_n, length(variances))]
  psidata_top <- psidata[top, , drop = FALSE]

  # PCA analysis
  psidata_top_t <- t(psidata_top)
  pca_res <- prcomp(psidata_top_t, scale. = TRUE)

  # get the sample info
  sample_info <- data.frame(
    Sample = rownames(psidata_top_t),
    Condition = gsub(paste0("^", exp_id, "_(.*)_PSI_R[0-9]+$"),
                     "\\1", rownames(psidata_top_t)),
    Replicate = gsub(".*_R([0-9]+)$", "R\\1", rownames(psidata_top_t))
  )

  pca_df <- data.frame(
    PC1 = pca_res$x[, 1],
    PC2 = pca_res$x[, 2],
    Condition = sample_info$Condition,
    Replicate = sample_info$Replicate
  )

  # Plot everything
  p <- ggplot(pca_df, aes(x = PC1, y = PC2, color = Condition)) +
    geom_point(size = 2, alpha = 0.8) +
    geom_text_repel(aes(label = Replicate), size = 4, show.legend = FALSE) +
    theme_minimal() +
    labs(
      title = paste0("PCA of GA‑exons, ", exp_id, " (", top_n, " most variable exons)"),
      x = "PC1",
      y = "PC2"
    ) +
    theme(plot.title = element_text(hjust = 0.5))

  print(p)
  invisible(p)
}


```

```{r}
#run analysis and plot per experiment
run_ga_pca("Exp1")
```

```{r}
run_ga_pca("Exp2")
```

```{r}
run_ga_pca("Exp3")
```

```{r}
run_ga_pca("Exp4")
```

```{r}
run_ga_pca("Exp5")
```