---
title: "Comparing STR profiles with STRprofilerR"
author: "Jared Andrews"
date: "`r Sys.Date()`"
output:
  BiocStyle::html_document:
    toc_float: true
vignette: >
  %\VignetteIndexEntry{Comparing STR profiles with STRprofilerR}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
```

# Introduction

Cell lines drift, get mislabelled, and cross-contaminate each other. Short
tandem repeat (STR) profiling is the standard way to catch it: genotype a panel
of highly polymorphic repeat markers, and compare the result against the
material the model was supposed to come from.

**STRprofilerR** compares STR profiles to each other or to a
reference database, scores the similarity three ways, flags samples that look
mixed, and can search the Cellosaurus knowledge base for a match. It is an R
port of the [STRprofiler](https://github.com/j-andrews7/STRprofiler) Python
package.

**STRprofilerR is intended for research purposes only.** A low score is a
prompt to investigate, not a verdict.

# Installation

You can install **STRprofilerR** from Bioconductor with:

```{r install, eval = FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("STRprofilerR")

# You can install the development version from Github
pak::pak("j-andrews7/STRprofilerR")
```

# Reading profiles

STR files typically come in one of two shapes (at least from our facility/platform), and `readSTRProfiles()` reads both. A **long**
file has one row per sample and one column per marker:

```{r read-long}
library(STRprofilerR)
longFile <- system.file("extdata", "ExampleSTR_long.csv", package = "STRprofilerR")
head(read.csv(longFile, check.names = FALSE))
```

```{r read-long-2}
profiles <- readSTRProfiles(longFile, sampleCol = "Sample Name")
profiles
```

A **wide** file has one row per sample/marker pair, with the alleles spread
across numbered columns. Columns whose names contain `Allele` are collected;
size and height columns are ignored. Point `markerCol` at the column containing marker names and
the layout is detected automatically.

Files of either shape, and of any supported format, can be read together. The
marker set is the union across files, and samples missing a marker are recorded
as untyped rather than as having no alleles. Untyped markers are skipped during scoring rather than counted as
mismatches.

```{r read-both}
xlsxFile <- system.file("extdata", "ExampleSTR.xlsx", package = "STRprofilerR")

if (requireNamespace("readxl", quietly = TRUE)) {
    both <- readSTRProfiles(c(longFile, xlsxFile), sampleCol = "Sample Name")
    print(both)
}
```

## Alleles are stored pre-split

Inside a `STRProfiles` object each marker is a
`r Biocpkg("IRanges")` `CharacterList`, so the alleles for a sample are already
a character vector rather than a string waiting to be split:

```{r alleles}
alleles(profiles)[["marker1"]]
```

Allele calls are cleaned on the way in: duplicates collapse, numeric alleles
sort ahead of string ones, and a trailing `.0` is dropped while a genuine
decimal such as `9.3` is kept.

```{r clean}
cleanAlleles("10.0,10,13,13.0,14,14 ")
```

## Marker names

Instruments and databases spell markers inconsistently. `markerAliases()` holds
the harmonisation table, applied by default:

```{r markers}
harmonizeMarkers(c("Penta D", "Penta_E", "vWA"))
```

Extend it with `extraAliases` when your instrument produces something the
default table does not cover.

Every marker is also classified at ingest, which is how amelogenin is found
without you having to name the column:

```{r marker-data}
markerData(profiles)
```

## Metadata

Database files often carry non-marker columns. Name them in `metadataCols` and
they end up in `sampleData()`, ignored for scoring purposes:

```{r metadata}
dbFile <- system.file("extdata", "main_database.csv", package = "STRprofilerR")
reference <- readSTRProfiles(dbFile)

dim(reference)
head(sampleData(reference), 3)
```

If using an existing database with extra sample metadata, be sure to set this so those columns aren't used in scoring.

This file is the database the Shiny application uses by default: profiles of
patient-derived models from The Jackson Laboratory PDX program and the NCI
Patient-Derived Models Repository. It's also provided as the
`pdxSTRProfiles` dataset used [below](#workflow).

# Scoring

Let $M$ be the markers where **both** profiles were typed, $s$ the number of
alleles they share across those markers, and $q$ and $r$ the number of alleles
each carries across those same markers:

$$\textrm{Tanabe} = 100 \times \frac{2s}{q + r}
\qquad
\textrm{Masters (query)} = 100 \times \frac{s}{q}
\qquad
\textrm{Masters (reference)} = 100 \times \frac{s}{r}$$

Tanabe is the Sørensen–Dice coefficient and is symmetric: swap the two profiles
and the score is unchanged. The Masters scores are not, which can be useful. Masters (query) asks how much of the query the reference accounts for,
which is useful when a query may be contaminated. A query containing everything the reference has *plus* extra
alleles scores 100 on Masters (reference) while Tanabe drops.

Amelogenin is a sex marker rather than a polymorphic STR, so it is left out
unless you ask for it with `useAmel = TRUE`.

```{r score}
queries <- readSTRProfiles(longFile, sampleCol = "Sample Name")
refs <- readSTRProfiles(
    system.file("extdata", "ExampleSTR_database.csv", package = "STRprofilerR"),
    sampleCol = "Sample Name"
)

scoreProfiles(queries, refs)
```

Scoring is vectorised as a handful of sparse matrix products rather than a loop
over pairs, so comparing a batch against a database of thousands is very quick. `chunkSize` controls how many queries are held in memory at once,
and `minScore` discards weak pairs as they are produced:

```{r score-filter}
scoreProfiles(queries, refs, minScore = 90)
```

# Comparing end to end

`compareProfiles()` runs scoring, flags mixing, and summarises.

```{r compare}
cmp <- compareProfiles(queries, refs)
cmp
```

The summary carries one row per query - the best two hits, and everything
passing each threshold:

```{r summary}
as.data.frame(summary(cmp))[, c("Sample", "mixed", "topHit", "nextBest")]
```

Leave `reference` out to compare a batch against itself, dropping
self-comparisons:

```{r compare-self}
compareProfiles(refs)
```

## Mixing

A diploid genome gives at most two alleles per autosomal marker. While models can get pretty messed up and may occasionally harbor a third allele, consistently
seeing three or more in many markers suggests more than one genome in the tube.

To help identify such cases, STRprofilerR's `flagMixedSamples()` counts markers with **more than two** alleles and flags a
sample when that count **exceeds** `threeAlleleThreshold` (3, by default).

```{r mixing}
batch <- readSTRProfiles(
    system.file("extdata", "Example_Batch_File.csv", package = "STRprofilerR")
)

flagMixedSamples(batch)
```

It is worth checking flagged samples more closely.

# Writing results

`writeSTRResults()` produces the file set the Python command line tool does,
prefixed `STRprofilerR` rather than `strprofiler`:

```{r write}
outDir <- file.path(tempdir(), "STRprofilerR-vignette")
written <- writeSTRResults(cmp, outDir)

basename(written)
```

That is a summary table, one table per query listing every reference it was
compared against with that reference's alleles, a self-contained HTML view of
the summary, and a log recording the parameters, inputs, and versions used.

# Querying Cellosaurus

When you have no local reference for a line,
[CLASTR](https://www.cellosaurus.org/str-search/) searches the human profiles in
Cellosaurus. Check your marker names against its controlled vocabulary first, as CLASTR silently ignores markers it does not recognise:

```{r validate}
validateClastrMarkers(markers(refs))
```

The query itself needs network access, so it is not run here:

```{r clastr, eval = FALSE}
profile <- c(
    Amelogenin = "X", CSF1PO = "11,12", D2S1338 = "19,23", D3S1358 = "15,17",
    D5S818 = "11,12", D7S820 = "10", D8S1179 = "10", D13S317 = "11,12",
    D16S539 = "11,12", D18S51 = "13", D19S433 = "14", D21S11 = "29,30",
    FGA = "20,22", PentaD = "11,13", PentaE = "14,16", TH01 = "6,9",
    TPOX = "8,9", vWA = "17,19"
)

hits <- clastrQuery(profile, scoreFilter = 90)
as.data.frame(hits)[, c("accession", "name", "score", "problem")]
```

A cell line with a known contamination history carries a note in `problem`. A result that CLASTR returns as two
profiles appears as two rows, with `(Best)` and `(Worst)` appended to the
accession.

`clastrBatchQuery()` writes CLASTR's own multi-sheet workbook instead, when you
want its presentation rather than the data.

# The command line application

The CLI is built on [Rapp](https://github.com/r-lib/Rapp). Install the launcher
once:

```{r cli-install, eval = FALSE}
install.packages("Rapp")
Rapp::install_pkg_cli_apps("STRprofilerR")
```

```bash
strprofilerr compare --database refs.csv -o ./results STR1.xlsx STR2.csv
strprofilerr clastr --score-filter 90 -o ./results batch.csv
strprofilerr app --database refs.csv
strprofilerr compare --help
```

Flags mirror the Python CLI, and both spellings are accepted, so
`--sample-col` and `--sample_col` both work and existing invocations keep
running.

# The Shiny application

For colleagues who would rather not write code, `STRprofilerApp()` serves the
same comparisons in a web browser. It is a port of the `strprofiler` Python
package's application and needs the `shiny`, `bslib`, and `DT` packages.

```{r app, eval = FALSE}
STRprofilerApp()

# Against your own database instead of the bundled one.
STRprofilerApp(database = "our_database.csv")
```

The application has four tabs:

- **Single Query**: type a profile in marker by marker and score it against
  the database or search Cellosaurus with it. Alleles differing from the query
  are highlighted in the report.
- **Batch Query**: upload a file of profiles and compare it against the
  database, within itself, or against Cellosaurus.
- **Database File Management**: swap in a custom database for the session.
- **Usage Guide**: how to read the reports, and what to cite.

From a shell, `strprofilerr app --database our_database.csv` does the same. To
host the application for others, deploy an `app.R` that calls
`STRprofilerR::STRprofilerApp()` with your database. One ships with the package,
at `system.file("app", "app.R", package = "STRprofilerR")`.

# STRprofilerR in a larger workflow {#workflow}

STR authentication is rarely an end in itself. It is a quality control step that
decides which models go on into experiments, and its results belong with the
data those experiments produce.

## Bioconductor data structures

`STRProfiles` is built from `r Biocpkg("S4Vectors")` and `r Biocpkg("IRanges")`
classes rather than bespoke ones. `sampleData()`, `markerData()`, `scores()`,
and `summary()` all return `DataFrame`s, and each marker's alleles are a
`CharacterList`, so the usual Bioconductor idioms for subsetting, merging, and
list-wise operations apply. Counting the alleles each sample carries at a
marker, for instance:

```{r bioc-classes}
lengths(alleles(batch)[["D5S818"]])
```

## The bundled PDX database

`pdxSTRProfiles` holds 1,258 STR profiles from The Jackson Laboratory PDX program and the NCI
Patient-Derived Models Repository, ready to use as a reference:

```{r pdx}
data(pdxSTRProfiles)
pdxSTRProfiles

pdxCmp <- compareProfiles(batch, pdxSTRProfiles)
as.data.frame(summary(pdxCmp))[, c("Sample", "mixed", "topHit", "nextBest")]
```

## Generating profiles from raw data

`STRProfiles()` provides the ability to build profiles
a `data.frame` already in memory, such as the result of a LIMS query, an
instrument export after some reshaping, or STR calls made from sequencing data
by another tool. It needs one row per sample and one column per marker, and
need not cover the full panel:

```{r from-df}
calls <- data.frame(
    Sample = "PDX_042",
    Amelogenin = "X", CSF1PO = "11", D13S317 = "9,14", D16S539 = "12",
    D18S51 = "17", D21S11 = "29,31", D3S1358 = "15", D5S818 = "11",
    D7S820 = "9,11", D8S1179 = "13,15", FGA = "21,22", TH01 = "9.3",
    TPOX = "8", vWA = "14"
)
newModel <- STRProfiles(calls)

scoreProfiles(newModel, pdxSTRProfiles, minScore = 80)
```

## Carrying authentication info with an experiment

Generally, authentication is frequently performed prior to larger experiments, and these profiles may be worth carrying along with said assay data. Recording each sample's authentication result in the
`colData` of its `r Biocpkg("SummarizedExperiment")` keeps the quality control
with the data, so it survives subsetting and goes wherever the object goes. The
same applies to the `colData` of a `SingleCellExperiment` or
`MultiAssayExperiment`.

In this example, the experiment's sample sheet records which model each sample should be,
and `Sample_C` has been swapped:

```{r se, eval = requireNamespace("SummarizedExperiment", quietly = TRUE), message = FALSE}
library(SummarizedExperiment)

# A stand-in for an experiment on the three models in the batch.
counts <- matrix(
    rpois(30, lambda = 50),
    nrow = 10,
    dimnames = list(paste0("gene", 1:10), rownames(batch))
)

se <- SummarizedExperiment(
    assays = list(counts = counts),
    colData = DataFrame(
        model = c("J000077451", "J000077960", "J000077451"),
        row.names = rownames(batch)
    )
)

# Each sample's best match in the reference.
pairs <- scores(pdxCmp)
pairs <- pairs[order(pairs$query, -pairs$tanabeScore), ]
best <- pairs[!duplicated(pairs$query), ]
best <- best[match(colnames(se), best$query), ]

se$strMatch <- best$reference
se$strTanabe <- best$tanabeScore
se$strMixed <- summary(pdxCmp)[colnames(se), "mixed"]
se$authenticated <- se$strMatch == se$model & se$strTanabe >= 80 & !se$strMixed
colData(se)

# Keep only authenticated samples for analysis.
se[, se$authenticated]
```

## Maintaining a reference database

Once a new model is authenticated, its profile could become a reference for the
next batch. `c()` combines `STRProfiles` objects, taking the union of their
markers, and `writeSTRProfiles()` writes a file that `readSTRProfiles()`, the
command line application's `--database`, and `STRprofilerApp()` can all read:

```{r grow-db}
inHouse <- c(pdxSTRProfiles, newModel)
inHouse

writeSTRProfiles(inHouse, file.path(tempdir(), "in_house_database.csv"))
```

This makes updating a reference database straightforward.

# Session info

```{r session}
sessionInfo()
```