---
title: "Pathogen-agnostic Sero-analysis Tutorial"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Pathogen-agnostic Sero-analysis Tutorial}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
options(rmarkdown.html_vignette.check_title = FALSE)
setup <- function() {
  needed <- c("knitr", "rmarkdown", "tidyverse", "kableExtra")
  
  lapply(needed, function(pkg) {
    if (requireNamespace(pkg, quietly = TRUE)) {
      library(pkg, character.only = TRUE)
    }
  })
}

setup()
library(SeroTrackR)
```

If you would like to run the `{SeroTrackR}` R package data processing pipeline for work not relating to the fit-for-purpose tool leveraging eight *Plasmodium vivax* serological exposure markers, then this is the tutorial for you! 

## Essential

Please ensure that you have read and prepared your raw Luminex files and plate layout files as per the instructions in the [setup](setup.html) vignette. 

## Step 1: Read in Raw Data 

For this tutorial, we will be using the built-in files within `{SeroTrackR}` that are extrapolated from *Plasmodium vivax* samples. 

```{r setup 1}
library(SeroTrackR)
library(tidyverse)

your_raw_data <- c(
  system.file("extdata", "example_MAGPIX_plate1.csv", package = "SeroTrackR"),
  system.file("extdata", "example_MAGPIX_plate2.csv", package = "SeroTrackR"), 
  system.file("extdata", "example_MAGPIX_plate3.csv", package = "SeroTrackR")
)
your_plate_layout <- system.file("extdata", "example_platelayout_1.xlsx", package = "SeroTrackR")
```

To run your OWN data, follow the code below and uncomment (i.e., remove the hashtags): 

```{r, eval=FALSE}
your_raw_data <- c(
  "PATH/TO/YOUR/FILE/plate1.csv",
  "PATH/TO/YOUR/FILE/plate2.csv",
  "PATH/TO/YOUR/FILE/plate3.csv"
)
your_plate_layout <- "PATH/TO/YOUR/FILE/plate_layout.xlsx"
```

You can read any Luminex serology data file using the `readSeroData()` function. 

The Luminex platform can be specified using the `platform` argument, with either `bioplex`, `magpix` or `intelliflex`. 

For the xPONENT software-based exported files (MAGPIX or INTELLIFLEX), the version of the software should be specified as either **"4.2"** or **"4.3"**, with the version **"4.2"** as the default.

```{r}
sero_data <- readSeroData(
  raw_data = your_raw_data, 
  platform = "magpix",        # default 
  version = "4.2"             # default 
)
```

This exports a list of data frames: 

1. `data_raw`: The raw, unedited files: A check that the correct file/s were added.
2. `results`: The processed serological data files: Antigen names compatible with the R package, clear columns for each antigen and their MFI results.
3. `counts`: The bead counts detected from the Luminex instrument per antigen, per sample, for each plate. 
4. `blanks`: The MFI results for the blank wells. 
5. `stds`: The MFI results for the standard curve samples.
6. `run`: The run information from the Luminex machine. 

You can access each data frame, for example the `results` data frame, like this:

```{r}
sero_data$results
```

Next, you need to read in the plate layout! This file will allow us to relabel the raw luminex data from "unknowns" to our study Sample ID's. 

```{r}
plate_list <- readPlateLayout(
  plate_layout = your_plate_layout, 
  sero_data  = sero_data
)
```

Great! Now our serological data has been processed and our plate layouts are added into R-compatible data frames ready for analysis!

## Step 2: Quality control checks

Quality control checks can also be run on any serological data file that has been processed using the `readSeroData()` and `readPlateLayout()` functions. The `runQC()` function will: 

1. Process bead count data (assigns a warning if bead counts are less than 15 or above 500 per antigen).
2. Assign Sample ID's to Luminex outputs.

```{r}
qc_results  <- runQC(
  sero_data = sero_data,       # load in your serological data variable 
  plate_list = plate_list      # load your plate list variable 
)
```

This exports a list of data frames: 

1. `processCounts_output`: The bead counts per antigen, per sample, for each plate. Intermediary file for downstream processing. Warning = 0 for pass and = 1 for fail. 
2. `getCounts_output`: The processed bead counts file for each sample, summarising whether a sample needs to be repeated. Intermediary file for downstream processing.
3. `sampleid_output`: Aligning the Luminex anonymised data (SampleID) to the plate layout ID's (Sample). Intermediary file for downstream processing.
4. `getAntigenCounts_output`: The processed bead counts file per antigen per sample. Used for `plotBeadCounts()`. 
5. `getCountsQC_output`: The bead counts and QC results per antigen. **Key file**. 

You can access each data frame, for example the `getCountsQC_output` data frame, like this:

```{r, eval=FALSE}
qc_results$getCountsQC_output
```

```{r, echo=FALSE}
qc_results$getCountsQC_output %>% head()
```

You can also visualise QC results using the `sero_data` and `qc_results` variables. 

### Standard Curve Plot

The `plotStds()` function plots the standard curve, generated from the antibody data from the standards you indicated in your plate layout (e.g. S1-S10) and Median Fluorescent Intensity (MFI) units are displayed in log10-adjusted scale. In general, the standard curves should look relativel linear (only when the y-axis is on log-adjusted scale). 

For your standard curve, you will need to specify the `location = "PNG"`.

*Note*: The argument was originally written to distinguish between PvSEM assays that used different pooled positive control plasma from highly immune adults, which the `location = "ETH"` argument requires a different processing pipeline - here this is irrelevant! 

In this plot, there are grey dots in the background. This only appears because we *are* using the PvSEMs in the fit-for-purpose tool and it signifies the standard curve range observed at WEHI. In your dataset, you will not have these grey dots. 

```{r}
plotStds(sero_data, location = "PNG", experiment_name = "experiment1")
```

For more information on how to interrogate and visualise standard curves, see the Explore Standard Curves page. 

### Bead Counts QC Plot

The `plotCounts()` function provides a summary of the bead counts for each plate well are displayed, with blue indicating there are sufficient beads (≥15) or red when there are not enough. If any of the wells are red, they should be double-checked manually and re-run on a new plate if required.

```{r}
plotCounts(qc_results, experiment_name = "experiment1")
```

```{r}
getRepeats(qc_results, plate_list)
```

The `getRepeats()` function will inform you whether there are "No repeats necessary" or provide a list of samples to be re-run. In the example data, the beads in plate 2 wells A1 and A2 will need to be repeated. 

For more information on how to visualise these data, see the Bead Counts page. 

### Blanks QC Plot

The Median Fluorescent Intensity (MFI) units for each antigen is displayed for your blank samples. In general, each blank sample should have ≤50 MFI for each antigen, if they are higher they should be cross-checked manually.

```{r}
plotBlanks(sero_data, experiment_name = "experiment1")
```

In the example data, blank samples recorded higher MFI values for LF005 on plate 1 and should be checked to confirm this is expected from the assay.

## Step 3: MFI to RAU conversion 

The automated data processing allows you to convert your Median Fluorescent Intensity (MFI) data into Relative Antibody Units (RAU) by fitting a 5-parameter logistic function to the standard curve on a per-antigen level using the `MFItoRAU()` function suite. 

**Key input variables**:
 
- `sero_data`, `plate_list`, `qc_results`: Previously defined variables from the `readSeroData()`, `readPlateLayout()` and `runQC()` functions. 
- `std_point`: Standard Point Curve: 5 = 5-point curve, 10 = 10-point curve. 

```{r, warning=FALSE, message=FALSE}
mfi_to_rau_output <- MFItoRAU(
  sero_data = sero_data, 
  plate_list = plate_list,
  qc_results = qc_results, 
  std_point = 10
)
```

The `mfi_to_rau_output` variable contains the processed MFI and RAU data for each antigen per sample. There are three data frames output from this variable: 

1. ``[[1]]``: Contains correct SampleID's, well number and letters, MFI and RAU values per antigen and QC results. 
1. ``[[2]]``: Contains correct SampleID's, MFI and RAU values per antigen and QC results. **Cleaner data frame!**
1. ``[[3]]``: Raw MFI to RAU conversion data per antigen. 

```{r, eval=FALSE}
mfi_to_rau_output[[2]]
```

```{r, echo=FALSE}
mfi_to_rau_output[[2]] %>% head()
```

The results from this log-log conversion should look relatively linear for each antigen and can be plotted using the `plotModel()` function suite. 

```{r}
plotModel(mfi_to_rau_output, sero_data)
```

For more information on how to visualise MFI and RAU data, see the MFI and RAU plots page. 

## Alternative Pipeline Function 

A separate function called `runPvSeroPipeline()` was introduced to simplify writing all of these functions. For the pathogen-agnostic case, please ensure that `classify = "No"` so that the PvSEM classification algorithm is not applied. 

```{r runPvSeroPipeline-without-classification, eval=FALSE}
no_classification_final_analysis <- runPvSeroPipeline(
  raw_data = your_raw_data, 
  plate_layout = your_plate_layout, 
  platform = "magpix", 
  location = "ETH", 
  experiment_name = "experiment1", 
  classify = "No", ########################## key if you do NOT want any classification performed i.e., you do not have PvSEM antigens 
  algorithm_type = "antibody_model", 
  sens_spec = "balanced"
)
```

You can then access the datasets and visualisations using the following functions: 

```{r, eval=FALSE}
# MFI and RAU data
no_classification_final_analysis[[1]] 

# QC plots
#### Standard Curve Plot
no_classification_final_analysis[[2]]

#### Bead Counts QC Plot
no_classification_final_analysis[[3]] # Plot
no_classification_final_analysis[[4]] # Samples to repeat 

#### Blanks QC Plot
no_classification_final_analysis[[5]]

#### Model Output Plot
no_classification_final_analysis[[6]]
```

## Final Code

```{r, eval=FALSE}
library(SeroTrackR)
library(tidyverse)

your_raw_data <- c(
  system.file("extdata", "example_MAGPIX_plate1.csv", package = "SeroTrackR"),
  system.file("extdata", "example_MAGPIX_plate2.csv", package = "SeroTrackR"), 
  system.file("extdata", "example_MAGPIX_plate3.csv", package = "SeroTrackR")
)
your_plate_layout <- system.file("extdata", "example_platelayout_1.xlsx", package = "SeroTrackR")

# read serology data
sero_data <- readSeroData(raw_data = your_raw_data, platform = "magpix", version = "4.2")

# read plate data
plate_list <- readPlateLayout(plate_layout = your_plate_layout, sero_data  = sero_data)

# run quality control checks 
qc_results  <- runQC(sero_data = sero_data, plate_list = plate_list)

# visualise quality control
plotStds(sero_data, location = "PNG", experiment_name = "experiment1")

# bead counts plot
plotCounts(qc_results, experiment_name = "experiment1")

# get repeats
getRepeats(qc_results, plate_list)

# blank well plot
plotBlanks(sero_data, experiment_name = "experiment1")

# perform MFI to RAU conversion 
mfi_to_rau_output <- MFItoRAU(
  sero_data = sero_data, 
  plate_list = plate_list,
  qc_results = qc_results, 
  std_point = 10
)

# visualise MFI to RAU conversion 
plotModel(mfi_to_rau_output, sero_data)
```