_comparison.Rmd

### {{ case }} vs {{ ctrl}}

```{r}
#| label: groupDiff{{random}}

# This function will run both QC steps (krsa_filter_lowPeps, krsa_filter_nonLinear) and krsa_filter_ref_pep
pep_passed_qc <- krsa_quick_filter(data = data_pw_max, data2 = data_modeled$scaled,
    signal_threshold = 5, r2_threshold = 0.8,
    groups = c("{{case}}", "{{ctrl}}"))

# This function calculates log2 fold change values between the defined groups
# The byChip argument lets you calculates the log2 fold change the results within each chip
differential_phosphorylated_peptides <- krsa_group_diff(data_modeled$scaled, c("{{case}}", "{{ctrl}}"), pep_passed_qc, byChip = T)

# save LFC table
write_csv(differential_phosphorylated_peptides, "results/dpp_{{case}}_{{ctrl}}-{{chip}}.csv")

# Extract top peptides based on the LFC cutoff using average of LFCs across chips
significant_peptides <-
  krsa_get_diff(differential_phosphorylated_peptides,
                totalMeanLFC ,
                c(0.15, 0.2, 0.3, 0.4)) %>% list("meanLFC" = .)

# Extract top peptides based on the LFC cutoff using per chip LFCs
sigPepsPerChip <- krsa_get_diff_byChip(differential_phosphorylated_peptides, LFC, c(0.15, 0.2, 0.3, 0.4))

# Combine the peptides hits in one list
sigPeps_total <- list(significant_peptides, sigPepsPerChip) %>% unlist(recursive = F) %>%  unlist(recursive = F)


```

#### Heatmap

After applying the *Filtering Parameters* for this group comparison, only *`r length(significant_peptides$meanLFC[["0.2"]])`*/141 peptides carried forward in the analysis (i.e. *`r length(significant_peptides$meanLFC[["0.2"]])` hits*). Below are some figures to visualize the differences between these samples for considering these *hits*.

```{r fig.height=8, fig.width=6}
#| label: heatmapInd{{random}}
#| fig.width: 6
#| fig.height: 8
#| fig-cap: "Random Caption"

# generates a heatmap using the selected groups and peptides
krsa_heatmap(data_modeled$normalized, significant_peptides$meanLFC$`0.2`, groups = c("{{case}}", "{{ctrl}}"), scale = "row")

```

#### Violin Plot

Below, the violin plot visualizes the distribution of selected peptides for the analysis.

```{r}
#| label: violinIndPlot{{random}}
#| fig-width: 6
#| fig-height: 6
#| fig-cap: "Violin plot of two groups"

# generates a violin plot using the selected groups and peptides
krsa_violin_plot(data_modeled$scaled, significant_peptides$meanLFC$`0.2`, "Barcode", groups = c("{{case}}", "{{ctrl}}"))

```

#### Waterfall Plot

This waterfall represents the log2 fold changes between the two groups at each peptide.

```{r}
#| label: waterfall{{random}}
#| fig-height: 8
#| fig-width: 6
#| fig-cap: "Waterfall Plot to show the distribution of change in peptide phosphorylation"

# generates a waterfall of the log2 fold change values for the selected peptide (top peptides)
krsa_waterfall(differential_phosphorylated_peptides, lfc_thr = 0.2, byChip = F)
```

#### Upstream Kinase Analysis

The lab carefully curated and mapped the kinases that can act and phosphorylate each peptide present on the chip. This was achieved by using multiple sources including GPS 3.0, Kinexus Phosphonet, PhosphoELM and PhosphoSite Plus. Based on that association between peptides and kinases, a random sampling analysis is performed for these hits. The basic idea of *KRSA* is: For each iteration (*2000* iterations performed in this analysis), the same number of hits are randomly selected from the total 141/or 193 peptides present on the chip. Predicted kinases are then mapped to this sample list of peptides and number of kinases are determined. The kinase count from the actual hits and random sampling is then compared to determine the significance.

```{r}
#| label: krsa{{random}}

# STK chip

if (params$chip_type == "STK") {
    chipCov <- KRSA_coverage_STK_PamChip_87102_v2
    KRSA_file <- KRSA_Mapping_STK_PamChip_87102_v1
} else if (params$chip_type == "PTK") {
    chipCov <- KRSA_coverage_PTK_PamChip_86402_v1
    KRSA_file <- KRSA_Mapping_PTK_PamChip_86402_v1
}


# For parallel computing, load the furrr package:
# opens multiple R sessions to run faster
plan(multisession)

# Run the KRSA function across the different sets of peptides using the furrr package for parallel computing
mutiple_krsa_outputs <-
  future_map(sigPeps_total, krsa, .options = furrr_options(seed = T))
saveRDS(mutiple_krsa_outputs, file = "multiple_krsa_output_{{case}}_{{ctrl}}_{{chip}}.RDS")

plan(sequential)

# Tidy output
df <- data.frame(matrix(unlist(mutiple_krsa_outputs), ncol = max(lengths(mutiple_krsa_outputs)), byrow = TRUE))
df <- setNames(do.call(rbind.data.frame, mutiple_krsa_outputs), names(mutiple_krsa_outputs$meanLFC.0.2))

df <- df %>% rownames_to_column("method") %>% select(Kinase, Z, method) %>% 
  mutate(method = str_extract(method, "\\w+\\.\\w+\\.\\w+")) %>% 
  mutate(method = gsub("(^\\w+)[\\.]", "\\1>", method)) %>% 
  mutate_if(is.numeric, round, 2)

df %>% 
  pivot_wider(names_from = method, values_from = Z) -> df2
  
# Creates an average Z score table using the across chip analysis
df %>% 
  filter(grepl("mean", method)) %>% 
  select(Kinase, Z, method) %>% group_by(Kinase) %>% mutate(AvgZ = mean(Z)) -> AvgZTable

# Creates an average Z score table using the within chip analysis
df %>% 
  filter(!grepl("mean", method)) %>% 
  select(Kinase, Z, method) %>% group_by(Kinase) %>% mutate(AvgZ = mean(Z)) -> AvgZTable2

AvgZTable2 |> 
  select(Kinase, AvgZ) |> 
  unique() |> 
  mutate(Direction = if_else(AvgZ >= 0, "Up", "Down")) |> 
  group_by(Direction) |> 
  nest() |> 
  mutate(filtered = map(data, \(x) {x |> mutate(AbsZ = abs(AvgZ)) |> arrange(desc(AbsZ)) |> slice_head(n = 10) |> select(-AbsZ)})) |> 
  select(-data) |> 
  unnest(filtered) |> 
  ungroup() |> 
  select(-Direction) |> 
  arrange(desc(AvgZ)) |> 
  knitr::kable()


# save file
AvgZTable2 %>% write_csv("results/withinChip_KRSA_Table_{{case}}_{{ctrl}}.csv")

# Extract top kinases based on abs(Z) score
kinases_hits <- AvgZTable2 |> 
  select(Kinase, AvgZ) |> 
  unique() |> 
  ungroup() |> 
  slice_max(AvgZ, n = 10) |> 
  pull(Kinase)

#krsa_top_hits(AvgZTable2, 1.75)
#krsa_top_hits(AvgZTable2, 1.5)

# Show the number of peptides per each set in atable
krsa_show_peptides(sigPeps_total) |> 
  knitr::kable()

```

#### Z Scores Plot

We will plot the individual and averaged Z scores using both the across and within chip analyses.

```{r}
#| label: zscoresPlot{{random}}
#| fig-height: 6
#| fig-width: 6
#| fig-cap: "Waterfall plot of the Z Scores of each kinase family"

# Generates Z scores waterfall plots
krsa_zscores_plot(AvgZTable2)

pdf(
  "combined_krsa_zscores.pdf",
  width = 8.5,
  height = 11,
  bg = "white"
)
AvgZTable2 |> ungroup() |> filter(abs(AvgZ) >= 2) |> krsa_zscores_plot()
dev.off()

```

#### Reverse KRSA Plot

We will use the reverse KRSA plot function, to plot the log2 fold change values for all peptides mapped to kinase hits. This will help us examine the activity of the kinase

```{r}
#| label: revKRSAPlot{{random}}
#| fig-height: 6
#| fig-width: 6
#| fig-cap: "Kinase Activity summary for each kinase family based on peptide phsophorylation"

# plot the reverse KRSA figure for top kinases to determine their activity levels
krsa_reverse_krsa_plot(chipCov, differential_phosphorylated_peptides, kinases_hits, 0.2, byChip = FALSE)
```

#### Coverage Plot

To view the coverage of kinases across the full list of peptides on the chip, we will use the coverage plot function

```{r}
#| label: covPlot{{random}}
#| fig-height: 6
#| fig-width: 6
#| fig-cap: "Percentge of peptides each kinase family phosphorylates"

# generates a kinase coverage plot 
krsa_coverage_plot(chipCov, AvgZTable2, chipType)

```

#### Ball Model Network

We will view the ball model network function, to generate a model representing the protein-protein interactions between kinases

```{r}
#| label: netPlot{{random}}
#| fig-width: 8
#| fig-height: 8

# Plot the network ball model
krsa_ball_model(kinases_hits, AvgZTable2, 10, 2.5, 4.8)


```

\newpage