Burkett-Run-All-STK.Rmd

---
title: "Mouse Empathy Dataset - Combined - STK"
subtitle: "Burkett/Jen"
date: "`r format(Sys.time(), '%d %B, %Y')`"
author: "Cognitive Disorders Research Lab"
format: 
  pdf:
    documentclass: "scrreport"
    fig-pos: 'htbp'
    classoptions: ["oneside", "letterpaper"]
    number-sections: true
    include-in-header: 
      text: |
        \usepackage{float}
editor: source
toc: true
params:
  chip_type: "STK" # Change the Chip Type to PTK if needed
  signal_file_r1: "kinome_data/STK/Run1_KRSAsOutput_Median_SigmBg_230714112932.txt" # Change the file name to the one you have
  saturation_file_r1: "kinome_data/STK/Run1_KRSAsOutput_Signal_Saturation_230714112932.txt" # Change the file name to the one you have
  signal_file_r2: "kinome_data/STK/Run2_KRSAsoutput_Median_SigmBg_230714113751.txt" # Change the file name to the one you have
  saturation_file_r2: "kinome_data/STK/Run2_KRSAsoutput_Signal_Saturation_230714113752.txt" # Change the file name to the one you have
  signal_file_r3: "kinome_data/STK/Run3_KRSAsOutput_Median_SigmBg_230714115014.txt" # Change the file name to the one you have
  saturation_file_r3: "kinome_data/STK/Run3_KRSAsOutput_Signal_Saturation_230714115014.txt" # Change the file name to the one you have
  threshold: 2 # Change the threshold if necessary
---

```{r}
#| label: "setup"
#| include: FALSE
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  cache = TRUE,
  echo = FALSE,
  message = FALSE,
  warning = FALSE
)

library(KRSA)
library(knitr)
library(tidyverse)

library(gt)
library(furrr)

run_prefix <- "run-all-STK"
```

# Introduction

## Background

The Pamstation12 instrument provides a profiling of kinase activity of cell or tissue samples. The device is loaded with either serine/threonine or tyrosine microarray chips. Each chip has 4 wells so four samples can be loaded on a single chip, and the Pamstation12 can accommodate 3 chips per run. The microarray represents 144 (STK chip) or 196 (PTK chip) reporter peptides that can be phosphorylated by serine/threonine or tyrosine kinases. The device measures the degree of the phosphorylation in real time by detecting fluorescently labeled antibodies at different exposure times. The list of peptides present in each microarray can be viewed here: [STK chip](https://pamgene.com/wp-content/uploads/2020/09/STK-144-PamChip-87102.pdf), [PTK chip](https://pamgene.com/wp-content/uploads/2020/09/PTK-196-PamChip-86402.pdf)

\newpage

![Pamgene Kinase Activity Platform](images/pamgene_workflow.png){fig-align="center"}

![Pamgene Platform Detection](images/pamgene_detectionFig.png){fig-align="center"}

# Run Design

Designing the placement of the samples on the chips and arrays is important to consider due to the variability across different chips and batches. During the run some wells are subject to fail and their data cannot be analyzed and shown below as red.

```{r}
#| label: "design"
#| fig-width: 4
#| fig-height: 4
#| eval: false

tt <- tibble(
  num = rep(c(3.5,2.5,1.5,0.5),3),
  chips = rep(c(2,5,8),each=4),
  colr = c(rep(1,12)),
  nms = rep(LETTERS[1:4], 3)
)
rundes <- ggplot() +
  geom_rect(aes(xmin=1, xmax=3, ymin=0,ymax=4), fill="gray30") +
  geom_rect(aes(xmin=4, xmax=6, ymin=0,ymax=4), fill="gray30") +
  geom_rect(aes(xmin=7, xmax=9, ymin=0,ymax=4), fill="gray30") +
  xlim(0,10) + 
  theme_void() 

rundes + geom_point(data = tt, aes(chips, num,color=factor(colr)), size = 13, show.legend = F) +
  geom_text(data = tt, aes(chips,num,label = nms), size = 1.6) +
  scale_color_manual(values=c("green3", "red3"))
```

# Results

## Image Analysis

The first step of analyzing the run is to convert the images taken by the PamStation of each array at different exposure times to numerical values This is done by the Bionavigator software developed by Pamgene. The software recognizes the grid of the array with the aid of the searching algorithm (Pamgrid) to correctly identify each spot on the array. The numbers produced by this software represent the median value of the foreground pixels minus the median value of the background pixels to produce the median signal minus background (Median_SigmBg).

## Reading Data

The first step will be reading the crosstab view bionavigator files (Median_SigmBg and Signal_Saturation) and defining the PamChip type (STK or PTK). The raw data is read and then transformed to be in tidy format for an easier analysis, modeling, and visualizing.

```{r}
#| label: "read-data"

chipType <- params$chip_type

# Read crosstab view bionavigator files and tidy data
raw_data_r1 <- krsa_read(
  signal_file = params$signal_file_r1,
  signal_saturation = params$saturation_file_r1
)

raw_data_r2 <- krsa_read(
  signal_file = params$signal_file_r2,
  signal_saturation = params$saturation_file_r2
)

raw_data_r3 <- krsa_read(
  signal_file = params$signal_file_r3,
  signal_saturation = params$saturation_file_r3
)

raw_data <- bind_rows(raw_data_r1, raw_data_r2, raw_data_r3)

```

## QC Initial Steps and Groups Assignments

We will perform a couple of quality control steps to deal with negative values in the data and adjust based on signal saturation (optional). Next, we will define a new column to represent the grouping. And then, we will extract end point signal values

```{r}
#| label: "quality-control"

# qc processing
qc_data <- krsa_qc_steps(raw_data) %>% 
  mutate(
  Group = `Comment 1`,
  SampleName = paste0(SampleName, "_" , Barcode)
)

groups <- unique(qc_data$Group)

# extract end level signal values @ all max exposure time (100ms)
krsa_extractEndPointMaxExp(qc_data, chipType) -> data_pw_max

# extract end level signal values end level @ all exposure times
krsa_extractEndPoint(qc_data, chipType) -> data_pw

```

## QC Steps and Model Fitting

We will filter out peptides with low signals. In order to combine the values from different exposure times into a single value, a simple linear regression model of the *Median_SigmBg* as a function of exposure time is fitted. The slope of of the model fit and $R^2$ are then used for quality control and samples comparison. The slope is also scaled by multiplying by 100 and log2 transformed (*Slope_Transformed*). We then filter out peptides with poor linear fit and references peptides.

```{r}
#| label: "model"

# Filter out peptides with low signals
high_peptides <- krsa_filter_lowPeps(data_pw_max, 5)

# Fit the linear model
# This will produce a list of data frames:
# scaled: the Slope_Transformed values (see above for more info on Slope_Transformed)
# normalized: the Slope_Transformed values but normalized by Chip/Barcode
# Grouped: The mean of Slope_Transformed values for all samples within a group
data_modeled <- krsa_scaleModel(data_pw, high_peptides)
  
# Filter out peptides weak linear fit
linear_peptides <- krsa_filter_nonLinear(data_modeled$scaled, 0.8)

# Filter out reference peptides
clean_peptides <- krsa_filter_ref_pep(linear_peptides)

saveRDS(data_modeled, str_glue("datastore/{run_prefix}-data_modeled-{params$chip_type}.RDS"))
```

## Global Signal Intensity

For a global signal intensity across all samples/groups, few figures can be plotted based on the *Slope_Transformed* values.

### Global CV Plots

We will plot the coefficient of variation on both the normal and normalized fits. This will help us to identify groups with high variation that could be explained by sample outliers.

```{r}
#| label: "global-cv-plot"
#| fig-width: 6
#| fig-height: 6
#| fig-cap: "Coefficient of Variation plotted for each peptide across all 4 groups"

# Plot a CV figure using the modeled normalized data
krsa_cv_plot(data_modeled$normalized, clean_peptides)
```

### Global Violin Plots

We will plot violin figures to examine global signal differences between groups/samples.

```{r}
#| label: "global-violin-plot-by-group"
#| fig-width: 6
#| fig-height: 6
#| fig-cap: "Violin Plots for signal intensity Distribution Across Groups for all replicates"

# Plot a violin figure and facet by the (Group) variable
krsa_violin_plot(data_modeled$scaled, clean_peptides, "Group")


```

```{r}
#| label: "global-violin-plot-by-chip"
#| fig-width: 6
#| fig-height: 6
#| fig-cap: "Violin Plots for signal intensity Distribution Across Chips for all replicates"

# Plot a violin figure and facet by the (Barcode) variable
krsa_violin_plot(data_modeled$scaled, clean_peptides, "Barcode")


```

### Global Heatmaps

The heatmap represent all the peptides present on the chip except the positive/internal controls and peptides that failed to pass QC. The heatmaps are scaled by row to highlight the peptide signal differences across the samples. A hierarchical unsupervised clustering is applied both on the peptides and the samples to group potentially similar signatures.

```{r}
#| label: "global-signal-heatmap-individual"
#| fig-cap: "Row and chip normalized intensity values for the selected peptides"
#| fig-height: 6
#| fig-width: 6

# Generates a heatmap using the modeled normalized data
krsa_heatmap(data_modeled$normalized, clean_peptides, scale = "row")

```

```{r}
#| label: "global-signal-heatmap-grouped"
#| fig-cap: "Row and group normalized intensity values for the selected peptides"
#| fig-height: 6
#| fig-width: 6
#| cache: true

# Generates a heatmap using the modeled grouped data
krsa_heatmap_grouped(data_modeled$grouped, clean_peptides, scale = "row")

pdf("combined_heatmap.pdf", width = 8.5, height = 11, bg = "white")
krsa_heatmap_grouped(data_modeled$grouped, clean_peptides, scale = "row")
dev.off()

```

## Group Comparison

To compare between samples, a two-group comparison is performed. In this case, there are three group comparisons:

- B vs A

The *Slope_Transformed* ratio between each group, paired by chip, is calculated to the fold change. Based on the fold change, peptides that pass a certain fold change threshold are considered significant hits. Also, quality control steps applied in each comparison to filter out peptides that do not reach specific criteria:

- The *Median_SigmBg* at max exposure *100ms* must be above a certain value\
- $R^2$ of the linear model fit must be above a threshold value

These *Filtering Parameters* (fold change threshold, QC criteria) can be modified to adjust the stringency of the analysis. The *Filtering Parameters* that are used for this analysis:

- The *Median_SigmBg* at max exposure *100ms* must be equal or above 5\
- $R^2$ of the linear model fit must be above or equal 0.8\
- Log fold change (LFC) cutoffs at (0.2,0.3,0.4)

```{r}
#| label: "group-comp-list"

# You can define the comparison groups here.
# The group names are in the `groups` variable and you can use subscripting to select the groups
# For example, to compare the first and fourth groups, you can use `groups[1]` and `groups[4]`
# Then copy the next block as many times as there are comparisons and change the `COMP1` to `COMP2`, `COMP3`, etc.
# Also make sure that the label is unique for each block

# Define Groups to be compared
comparisons <- list(COMP1 = c(groups[1], groups[2]))
```

```{r}
#| results: "asis"
#| echo: FALSE
#| include: TRUE
#| label: "childA"
#| fig-width: 6
#| fig-height: 6
#| fig-cap: "Heatmap of pairwise group comparison"

case <- comparisons$COMP1[1]
ctrl <- comparisons$COMP1[2]
random <- "A"
threshold <- params$threshold

src <- knitr::knit_expand(file = "_comparison.Rmd", case = case, ctrl = ctrl, threshold = threshold, chip = params$chip_type)

res <- knitr::knit_child(text = src, quiet = TRUE)

cat(res, sep = "\n")
```

# Session Info

```{r session_info}
devtools::session_info()
```