# Validation of DNA Interpretation Systems {#dnaval}

#### Sam Tyner{-}

### Acknowledgements {-}

This work would not have been possible without the excellent documentation of the [`strvalidator`](https://sites.google.com/site/forensicapps/strvalidator) package [@strval]. Thank you to the package's author, [Oskar Hansson, Ph.D](https://www.linkedin.com/in/goto-oskarhansson/), who has authored many, many supporting documents, tutorials, etc. for his `strvalidator` package. Thank you, Oskar! 

```{r setup03, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, cache=TRUE, fig.align = "center", message = FALSE, warning = FALSE, tidy = TRUE)
```

```{r dnaimage, echo = FALSE, out.width="100%"}
knitr::include_graphics("img/dna.png")
```

## Introduction

The earliest documented use of [DNA profiling](#def:dnaprof) in the legal system was an immigration dispute in the United Kingdom [@butler05]. A young man of Ghanaian descent with family in the UK was believed to have forged his Ghanaian passport and had an expired British passport^[See the interview (https://www.bbc.co.uk/sounds/play/b008fcdz) at 19:23.]. DNA profiling techniques developed by Sir Alec Jeffreys were used to prove that the he was indeed his mother's son, and thus he did have a right to immigrate to the UK. The technique was subsequently used for many other parentage cases, and soon after, DNA profiling was used for the first time to convict someone of murder in 1986 [@butler09].

When DNA profiling began, an individual's blood sample was taken to create their DNA profile. Now, DNA can be taken by a cheek swab, and minute traces of touch DNA can tie a perpetrator to the scene of the crime. This is thanks to the [polymerase chain reaction](#def:pcr) (PCR), a method of copying a DNA sample over and over again to amplify the genetic signal for profile extraction. Once a DNA sample is amplified with PCR, different DNA markers can be analyzed to make an identification. The standard for forensic DNA typing is to use [short tandem repeats](#def:str) (STRs) as the DNA marker. Other markers, [single nucleotide polymorphisms](#def:snp) (SNPs) and the mitochondrial genome ([mtDNA](#def:mtdna)), have different uses. SNPs can be used to identify ancestry or visible traits of a human, while mtDNA is used in cases where DNA is highly degraded [@dnareview]. Because STR is the standard, we dedicate the rest of this chapter to its methodology. 

### Procedure for DNA Analysis using STRs

In order to understand the STR methodology, we first need to understand what is being analysed. We present the comparison of genetic and printed information from @butler09 in Table \@ref(tab:comparetable). When forensic scientists analyze a DNA sample, they are looking for repeated "words" or DNA sequences in different "paragraphs," or loci. The locus information is stored in the [chromosome](#def:chromo), which is the "page" the genetic information is on. Your *chromosomes* exist in the *nucleus* of every *cell* in your *body*, just like a *page* is within a *chapter* in a *book* in a *library*.  STR markers are a set of loci or genes. At each locus, the number of times a tetranucleotide sequence (e.g. AAGC) repeats is counted [@butler05]. This count indicates the [allele](#def:allele), or gene variation, at that particular locus. 

```{r comparetable, results='asis', echo = FALSE}
x <- c("Library", "Book", "Chapter", "Page", "Paragraph", "Word", "Letter")
y <- c("Body", "Cell", "Nucleus", "Chromosome", "[Locus](#def:locus) or gene", "Short DNA sequence", "DNA [nucleotide](#def:nucleotide)")
knitr::kable(data.frame(x,y), col.names = c("Printed Information", "Genetic Information"), caption = "Recreation of Table 2.1 from @butler09")
```

In forensic DNA profiling, a particular set loci are examined for comparison. The number of loci depends on the equipment and method used to analyzed the DNA sample, but can be as high as 27 for the particular method we discuss here [@loci27]. As of January 1, 2017, there are 20 core loci in CODIS, the Combined DNA Index System, which is the FBI's national program for DNA [databases](#def:database) and software used in criminal justice systems across the United States [@codis]. These sets of loci were chosen because of their high variability in the population. To find the alleles at each loci, the DNA sample is amplified using PCR, and then run through [capillary electrophoresis](#def:ce) (CE). The result of CE is the DNA profile, with the alleles on each locus indicated by different colored peaks from a chemical dyeing process. <!--For interpreting the results of the chemical analysis, the loci need to be described in terms of their stutter ratio, heterozygote peak balance, inter-locus balance, and stochastic threshold. --> 

The amplification process introduces random change known as *slippage*, which creates [*stutter peaks*](#def:stutter) in the observed DNA profile that are different than the true allele peaks [@butler09]. In addition, different labs may use different machines and materials in forensic analysis resulting in different measurements for the same DNA sample. Thus, the validation of methods and materials is a crucial step. According to the Scientific Working Group on DNA Analysis Methods (SWGDAM), "validation is a process by which a procedure is evaluated to determine its efficacy and reliability for forensic casework and/or database analysis" [@swgdamval].  Validation helps minimize error in forensic DNA analysis and helps keep results consistent across laboratories and materials.

The process of validation for forensic DNA methodology is expensive, time consuming, and unstandardized, and the R package `strvalidator` was created to help solve these problems in forensic DNA analysis [@strval]. The `strvalidator` package makes the validation process faster by automating data analysis with respect to "heterozygote balance, stutter ratio, inter-locus balance, and the stochastic threshold" [@strval]. In the remainder of this chapter, we introduce the type of data to import for use of this package, the primary functions of the package, and show an example of each of the four aforementioned validation steps in R.  


## Data 

The `strvalidator` package takes files exported from the [GeneMapper®](https://www.thermofisher.com/order/catalog/product/4475073) software, or a similar expert system that exports tab-delimited text files, as inputs. The data exported from these software programs typically come in a very wide format, and on import it needs to be transformed into a long format more appropriate for data analysis. In Figure \@ref(fig:widelonggif), we visualize the process of data being transformed from wide to long format and back. In wide format, variable values are column names, while in the long format these column names become part of the data.

```{r widelonggif, out.width = "80%", fig.align="center", fig.cap="Animation heuristic showing the transformation from long form to wide form data and back. Code for GIF from [Omni Analytics Group](https://omnianalytics.io/2018/08/30/animating-the-data-transformation-process-with-gganimate/)", echo = FALSE}
# code to recreate the gif is in the file "code/dna-animation-scratch.R"
knitr::include_graphics("img/dna-anim.gif")
```

The `strvalidator` package contains import methods to make sure that the data imported from other software is in the right form for validation analysis. The GeneMapper® software creates one column for each possible allele observed at a locus and their corresponding sizes, heights, and data points. Once the data have been trimmed and slimmed, they look something like this: 

```{r echo=FALSE}
library(strvalidator)
data(set4)
myDNAdata <- set4
```
```{r}
head(myDNAdata)
```

where

- `Sample.Name` is the name of the sample being analyzed 
- `Marker` is the locus in the DNA analysis kit 
- `Dye` is the dye channel for that locus
- `Allele` is the allele (# of sequence repeats) at that locus 
- `Height` is the observed peak height after amplification in RFUs (RFU = Relative Fluorescence Unit)


## R Package

The package `strvalidator` has a [graphical user interface](#def:gui) (GUI) to perform analyses so that no coding knowledge is necessary to run these analyses. The author of the package, Oskar Hansson, has written an extensive tutorial^[*STR-validator* 2.0.0 tutorial, published 8 August 2017 (https://drive.google.com/file/d/0B2v6NDpFIgvDel9wSlFfb3hvekk/view).] on the GUI. As this book is focused on open science, we do not use the GUI because it does not output the underlying code used for the point-and-click analyses. Instead, we use the code that powers the GUI directly. This code is called the "basic layer" of the package by @strval.

The data are read into R via the `import()` function. This function combines the processes of trimming and slimming the data. Trimming selects the columns of interest for your analysis (e.g. `Sample.Name`, `Allele`, `Height`), while slimming converts the data from wide format to long format, as shown in Figure \@ref(fig:widelonggif). 

After the data has be loaded, there are four main families of functions in the `strvalidator` package that are used for analysis. 

- `add*()`: Add to the DNA data. For example, use `addMarker()` to add locus information or `addSize()` to add the fragment size in base pair (bp) for each allele.
- `calculate*()`: Compute properties of the DNA data. For example, use `calculateHb()` to compute heterozygous balance for the data  or `calculateLb()` to compute the inter-locus balance (profile balance) of the data. 
- `remove*()` :  Remove artifacts from the data with `removeArtefact()` and remove spikes from the data with `removeSpike()`. 
- `table*()` : Summarize results from one of the `calculate*()` analyses. For example, `tableStutter()` summarizes the results from `calculateStutter()`. 

For complete definitions and explanations of all functions available in `strvalidator`, please see the [`strvalidator` manual](https://drive.google.com/file/d/0B2v6NDpFIgvDMzlydllKaHVBYW8/view).  There are many other capabilities of `strvalidator` that do not cover in this chapter for the sake of brevity. 

## Drawing Conclusions

There is no one tidy way to conclude a DNA validation analysis, which may be done for new machines, new kits, or any internal validation required. The `strvalidator` package's primary purpose is to *import* large validation data sets and *analyze* the results of the validation experiment according to different metrics [@riman16]. A more complete description of the necessary validation studies is found in @swgdamval, and full step-by-step tutorials can be found in @riman16 and @strvalweb. 

For validation analysis with respect to heterozygote balance, stutter ratio, inter-locus balance, and stochastic threshold, there are recommended guidelines to follow. 

### Stutter ratio 

The `*Stutter()` functions in `strvalidator` can analyze ratios of different types of stutter such as the backward stutter, the forward stutter, and the allowed overlap (none, stutter, or allele), as shown in Figure \@ref(fig:stutterfig). Each of @hill2011, @westen2012, @brookes2012, and @tvedebrink2012 show greater stutter with more repeats, and these results are similar to those in @strval. In addition they found that some loci, such as TH01, experience less stutter on average than others. 

### Heterozygote balance

For guidelines specific to the PowerPlex® ESX 17 and ESI 17 systems featured in @strval, refer to @hill2011. Generally speaking, per @gill97, the heterozygote balance should be no less less than 60%. 

###  Inter-locus balance

Per @strval, there are two methods in `strvalidator` to compute inter-locus balance.

1. As the proportion of the total peak height of a profile 
2. Relative to the highest peak total within a single locus in the profile, with the option to compute this value for each dye channel. 
    
Ideally, the loci would be perfectly balanced and the total peak height in each locus would be equal to $\frac{1}{n}$ where $n$ is the number of loci in the kit [@strval]. 


### Stochastic threshold

The *stochastic threshold* (ST) or *interpretation threshold* is the "point above which there is a low probability that the second allele in a truly heterozygous sample has not dropped out" [@butler09]. The ST is used to assess dropout risk in `strvalidator`. Another important threshold in DNA interpretation is the analytical threshold (AT), which is a peak height (for example, 50 RFUs) above which peaks "are considered an analytical signal and thus recorded by the data analysis software" [@butler09].  @strval refer to analytical threshold (AT) as the limit of detection threshold (LDT). Peaks above the AT are considered *signal*, and any peaks below the AT are considered *noise*. The ST is the RFU value above which it is reasonable to assume that, at a given locus, allelic dropout of a sister allele has not occurred. .^[ThermoFisher has a great five minute video explainer of the two different threshold types (https://youtu.be/lbCQmMbfIoM).] Peaks that appear to be homozygous but have heights above the AT and below the ST may not be true homozygotes and may have experienced stochastic effects, such as allele dropout or elevated stutter. Usually, these stochastic events only happen for very small amounts of DNA that have been amplified.

In `strvalidator`, dropout is scored according to the user-provided LDT value and the reference data provided. The risk of dropout is then modeled using a logistic regression of the calculated dropout score on the allele heights. Then for an acceptable level of dropout risk, say 1%, the stochastic threshold is computed according to the logistic regression model. Thus, the ST is the peak height at which the probability of dropout is less than or equal to 1%.  


## Case Study

We do a simple case study using eight repeated samples from the same individual that are included in the `strvalidator` package. 

<!--
There are four pieces in validation: 
1. stutter ratio 
2. heterozygote balance (intra-locus balance)
3. inter-locus balance 
4. stochastic threshold. 
-->

### Get the data

We'll use the package data `set1`, which is data from the genotyping of eight replicate measurements of a  positive control sample, one replicate of a negative control sample, and the ladder used in analysis. The PowerPlex® ESX 17 System from the Promega Corporation^[Technical manual: https://www.promega.com/-/media/files/resources/protocols/technical-manuals/101/powerplex-esx-17-system-protocol.pdf?la=en] was used on these samples for amplification of 17 loci recommended for analysis by the European Network of Forensic Science Institutes (ENFSI) and the European DNA Profiling Group (EDNAP), the European equivalent of SWGDAM. The known reference sample used is the `ref1` data in the `strvalidator` package. 

First, we load the data, then slim it for analysis. Then, we use `generateEPG()` to visualize an electropherogram-like plot of the data. This function, like the other plotting functions in `strvalidator`, is built on the `ggplot2` package [@R-ggplot2]. We also use the [`dplyr`](https://dplyr.tidyverse.org/) package throughout for data manipulation tasks [@R-dplyr].   

```{r set1, fig.show='hide'}
library(strvalidator)
library(dplyr)
library(ggplot2)
data(set1)
head(set1)
# slim and trim the data
set1.slim <- slim(set1, fix = c("Sample.Name", "Marker", "Dye"),
                  stack = c("Allele", "Height"), keep.na = FALSE)
dim(set1)
dim(set1.slim)
head(set1.slim)
p <- set1.slim %>% 
  filter(Sample.Name != "Ladder") %>% 
  generateEPG(kit = "ESX17")
```
```{r epg, fig.cap="Electropherogram-like `ggplot2` plot of the mean of all 8 samples in `set1`"}
p +  ggtitle("Mean peak heights for 8 samples from PC shown")
```

Next, get the reference sample data. 

```{r ref1, fig.show='hide'}
data(ref1)
head(ref1)
ref1.slim <- slim(ref1, fix = c("Sample.Name", "Marker"), 
                  stack = "Allele", keep.na = FALSE)
head(ref1.slim)
p <- generateEPG(ref1.slim, kit = "ESX17") + ggtitle("True profile for sample PC") 
```

```{r ref1epg, fig.cap="The reference profile electrogpherogram, `ref1`."}
p
```

```{r oldplot,echo=FALSE, eval=FALSE}
alevels <- as.character(sort(as.numeric(unique(set1.slim$Allele)[-c(1:3)])))
plotsamples <- set1.slim %>% filter(Sample.Name != "Ladder", Allele != "OL", Marker != "AMEL") %>% sortMarker(kit = "ESX17")
plotsamples$Allele <- factor(plotsamples$Allele, levels = alevels, ordered = T)
plotsamples$Height <- as.numeric(plotsamples$Height)
plotref <- ref1.slim %>% sortMarker(kit = "ESX17") %>% filter(Marker != "AMEL")
plotref$Allele <- factor(plotref$Allele, levels = alevels, ordered = T)
ggplot() + 
  geom_line(data = plotsamples, aes(x = Allele, y = Height, group = Sample.Name, 
                                    color = Sample.Name)) + 
  geom_segment(data = plotref, aes(x = Allele, xend = Allele, y = -Inf, yend = Inf), alpha = .7) + 
  scale_color_brewer(palette = "Set1")+
  guides(color = guide_legend(nrow = 1))+
  facet_wrap(~Marker, ncol = 4, scales= "free", drop = T) + 
  theme(legend.position = "top", axis.text.x = element_text(size = rel(.8), angle = 270, hjust = 0, vjust = .5)) +
  labs(x = "Allele (True allele shown in black)", y = "Peak Height (RFUs)", color = "Sample") 
plotsamples <- set1.slim %>% filter(Sample.Name != "Ladder", Allele != "OL", Marker == "AMEL") %>% sortMarker(kit = "ESX17")
plotref <- ref1.slim %>% sortMarker(kit = "ESX17") %>% filter(Marker == "AMEL")


```



### Check the stutter ratio 

```{r stutterfig, echo=FALSE, fig.align="center", out.width="50%", fig.cap="Figure 2 from @strval. The analysis range, 2 back stutters and 1 forward stutter is shown at 3 levels of overlap."}
knitr::include_graphics("img/strval-fig2.jpg")
```

Stutter peaks are byproducts of the DNA amplification process, and their presence muddles data interpretation [@strval]. Stutter is caused by strand slippage in PCR [@butler09]. This slippage causes small peaks to appear next to true peaks, and a threshold is needed to determine if a peak is caused by slippage or if it could be a mixture sample with a minor contributor. We calculate the stutter for the eight replicates in `set1` using one back stutter, no forward stutter and no overlap. We compare these values to the 95$^{th}$ percentiles in Table 3 of @strval. See Figure \@ref(fig:stutterfig) for an example of stutter.   


```{r stutter1, fig.cap="Stutter ratios by allele for each of the eight samples in the `set1` data, computed for one back stutter, zero forward stutter, and no overlap. Note that SR increases with allele length (e.g. D10S1248; D2S1338; D12S391). Horizontal dotted lines represent the 95th percentile of stutter ratio values from the study done in @strval."}
# make sure the right samples are being analyzed 
checkSubset(data = set1.slim, ref = ref1.slim)
# supply the false stutter and true stutter values for your data. these are from the GUI. 
stutter_false_val <- c(-1.9, -1.8, -1.7, -0.9, -0.8, -0.7, 0.9, 0.8, 0.7)
stutter_replace_val <- c(-1.3, -1.2, -1.1, -0.3, -0.2, -0.1, 0.3, 0.2, 0.1)
# calculate the stutter values 
set1_stutter <- calculateStutter(set1.slim, ref1.slim, back = 1, forward = 0, interference = 0, 
                 replace.val = stutter_false_val, by.val = stutter_replace_val)

stutterplot <- addColor(set1_stutter, kit = "ESX17") %>% 
  sortMarker(kit = "ESX17", add.missing.levels = FALSE)
marks <- levels(stutterplot$Marker)[-1]
stutterplot$Marker<- factor(as.character(stutterplot$Marker), levels = marks)
compare_dat <- data.frame(Marker = ref1$Marker[-1], 
                          perc95 = (c(11.9, 4.6, 10.9, 10.7,12.1, 12.0, 11.1, 10.4,16.0,
                                     11.4, 9.1, 10.1, 8.3, 14.4, 10.1, 12.8))/100)
compare_dat <- filter(compare_dat, Marker %in% stutterplot$Marker)

ggplot() + 
  geom_point(data = stutterplot, position = position_jitter(width = .1), 
             aes(x = Allele, y = Ratio, color = as.factor(Type)), alpha = .7) + 
  geom_hline(data = compare_dat, aes(yintercept = perc95), linetype = "dotted") + 
  facet_wrap(~Marker, ncol=4, scales = "free_x", drop = FALSE) + 
  labs(x = "True Allele", y= "Stutter Ratio", color = "Type")
```

Figure \@ref(fig:stutter1) shows the ratio of stutter for each of the eight control samples in `set1`. The horizontal dotted lines show the 95$^{th}$ percentile of the stutter ratio values computed in the same way from 220 samples in @strval. There are a few stutter values above the dotted line,  but overall the values correspond to what we expect to happen in a sample with only one contributor. Unusual values are shown in Table \@ref(tab:stuthigh).

```{r stutterbig, echo = FALSE}
stutterplot %>% left_join(compare_dat) %>% mutate(large = Ratio > perc95) %>% 
  filter(large) %>% 
  select(Sample.Name:Type, `95th perc.` = perc95) %>% 
  knitr::kable(caption = "Stutter peaks larger than the 95$^{th}$ percentile of peak values for the study in @strval.", label = "stuthigh", digits = 3)
```

### Check heterozygote balance (intra-locus balance)

Computing the heterozygote peak balance (Hb) is most important for analyzing samples with two or more contributors. We calculate Hb values for the eight repeated samples in `set1` below using Equation 3 from @strval to compute the ratio. 

```{r hb, fig.cap="Hb ratio values for the eight samples in `set1`. Most ratios are above the 0.6 threshold."}
# checkSubset(data = set3, ref = ref3)
set1_hb <- calculateHb(data = set1.slim, ref = ref1.slim, hb = 3, kit = "ESX17", 
            sex.rm = TRUE, qs.rm = TRUE, ignore.case = TRUE)
hbplot <- addColor(set1_hb, kit = "ESX17") %>% 
  sortMarker(kit = "ESX17", add.missing.levels = FALSE)
hbplot$Marker <- factor(as.character(hbplot$Marker), levels = marks)


ggplot(data = hbplot) + 
  geom_point(aes(x = MPH, y = Hb, color = Dye), position = position_jitter(width = .1)) + 
  geom_hline(yintercept = .6, linetype = "dotted") + 
  facet_wrap(~Marker, nrow = 4, scales = "free_x", drop = FALSE) + 
  scale_color_manual(values = c("blue", "green", "black", "red")) + 
  labs(x = "Mean Peak Height (RFU)", y = "Ratio", color = "Dye") + 
  guides(color = guide_legend(nrow=1)) + 
  theme(axis.text.x = element_text(size = rel(.8)),  legend.position = 'top')
```

Figure \@ref(fig:hb) shows the Hb values for the eight samples in `set1`. The balance ratio is typically no less than 0.6 according to @gill97, but there are a few exceptions to this rule in the `set1` sample, shown in Table \@ref(tab:smhb)

```{r which6, echo = FALSE}
set1_hb %>% filter(Hb < .6) %>% knitr::kable(caption = "Observations in the `set1` data which have Hb value less than 0.6.", label = "smhb", digits = 3)
```

<!--Another way to look at the heterozygote balance is with a box plot shown in Figure \@ref(fig:hbbox). 

```{r hbbox}
ggplot(data = hbplot, aes(x = Marker, y = Hb, color = Dye)) + 
  geom_boxplot(alpha = .7) + 
  geom_hline(yintercept = .6, alpha = .6, linetype = "dotted") + 
  scale_color_manual(values = c("blue", "green", "black", "red")) + 
  labs(x = "Locus", y = "Ratio") + 
  theme(axis.text.x = element_text(size = rel(.8)),  legend.position = "top")

```
-->

### Check inter-locus balance

Inter-locus balance (Lb) is a measure of peak balances across loci [@strval]. The total height of the peaks in all loci should be spread evenly across each individual locus in a sample. In the `set1` data, 17 loci are measured, thus each individual locus balance should be about $\frac{1}{17}^{th}$ of the total height of all peaks in RFUs. 

```{r lb, fig.cap="Inter-locus balance for the eight PC samples. At each locus, the value should be about 1/17. The peak heights should ideally be similar in each locus."}
set1_lb <- calculateLb(data = set1.slim, ref = ref1.slim, kit = "ESX17", option = "prop", 
            by.dye = FALSE, ol.rm = TRUE, sex.rm = FALSE, qs.rm = TRUE,  
            ignore.case = TRUE, na = 0)
set1_height <- calculateHeight(data = set1.slim, ref = ref1.slim, kit = "ESX17", 
                               sex.rm = FALSE, qs.rm = TRUE,na.replace = 0)
set1_lb %>% left_join(
  set1_height %>% 
    select(Sample.Name:Marker, Dye, TPH, H, Expected, Proportion) %>% distinct(), 
  by = c("Sample.Name", "Marker", "Dye", "TPPH" = "TPH")
) -> set1_lb


set1_lb <- sortMarker(set1_lb, kit = "ESX17", add.missing.levels = TRUE)

ggplot(set1_lb) + 
  geom_boxplot(aes(x = Marker, y = Lb, color = Dye), alpha = .7) + 
  scale_color_manual(values = c("blue", "green", "black", "red")) + 
  geom_hline(yintercept = 1/17, linetype = "dotted") +
  theme(legend.position = "top", axis.text.x = element_text(size = rel(.8), angle = 270, hjust = 0, vjust = .5)) + 
  labs(y = "Lb (proportional method)")
```

<!--
```{r, eval = FALSE}

ggplot(set1_lb) + 
  geom_line(aes(x = Marker, y = Lb, group = Sample.Name, color = Sample.Name), alpha = .7) + 
  scale_color_brewer(palette = "Set2") + 
  geom_hline(yintercept = 1/17, linetype = "dotted") +
  theme(legend.position = "top", axis.text.x = element_text(size = rel(.8), angle = 270, hjust = 0, vjust = .5))
```
--> 

The inter-locus balance for this kit should ideally be about $\frac{1}{17} \approx 0.059$. This value is shown by the horizontal dotted line in Figure \@ref(fig:lb). However, the markers in the red dye channel have consistently higher than ideal peaks and those in the yellow channel have consistently lower than ideal peaks. 

### Check stochastic threshold

The stochastic threshold is the value of interest for determining allele drop-out. If a peak is above the stochastic threshold, it is unlikely that an allele in a heterozygous sample "has dropped out" [@butler09]. Allele drop-out occurs when the allele peak height is less than the limit of detection threshold (LDT). As recommended in @butler09, we use an LDT of 50. The stochastic threshold is modeled with a logistic regression.  

```{r sthresh}
set1_do <- calculateDropout(data = set1.slim, ref = ref1.slim, threshold = 50, 
                            method = "1",
                            kit = "ESX17")
table(set1_do$Dropout)
```

In `set1`, there is no dropout, as the samples included are control samples, and thus enough DNA is present during amplification so there are no stochastic effects.  

For a more exciting dropout analysis, we use another data set with more appropriate information. The data `set4` was created specifically for drop-out analysis, and contains 32 samples from three different reference profiles. The `method = "1"` argument computes dropout with respect to the low molecular weight allele in the locus.

```{r sthresh2}
data(set4)
data(ref4)
set4_do <- calculateDropout(data = set4, ref = ref4, threshold = 50, 
                            method = "1",
                            kit = "ESX17")


table(set4_do$Dropout)
```

In the `set4` data, 33 alleles dropped out (`Dropout = 1`), and locus dropout (`Dropout = 2`) occurred in 9 samples (68 alleles). In one sample, all loci dropped out, while only one locus dropped out in three samples. The locus which most commonly dropped out was D22S1045 in seven samples, while loci D19S433 and D8S1179 only dropped out in two samples each. 

The probability of allele drop-out is computed via logistic regression of the dropout score with respect to the method `1`, on the the height of the allele with low molecular weight. The model parameters are also computed using the `calculateT()` function. This function also returns the smallest threshold value at which probability of dropout is less than or equal to a set value, typically 0.01 or 0.05, as well as a conservative threshold, which is the value at which the risk of observing a drop-out probability greater than the specified threshold limit is less than the set value of 0.01 or 0.05. 

```{r sthreshcalc}
set4_do2 <- set4_do %>% filter(Dropout != 2) %>% 
  rename(Dep = Method1, Exp = Height)
do_mod <- glm(Dep ~ Exp, family = binomial("logit"), 
              data =set4_do2)
set4_ths <- calculateT(set4_do2,pred.int = .98)
```

Next, we compute predicted dropout probabilities $P(D)$ and corresponding 95% confidence intervals and plot the results. 

```{r predictdo, fig.cap="Probability of dropout in `set4` for peaks from 100-1500 RFUs. 95% confidence interval for drop-out probability shown in red."}
xmin <- min(set4_do2$Exp, na.rm = T)
xmax <- max(set4_do2$Exp, na.rm = T)
predRange <- data.frame(Exp = seq(xmin, xmax))
ypred <- predict(do_mod, predRange, type = "link", se.fit = TRUE)
# 95% prediction interval
ylower <- plogis(ypred$fit - qnorm(1 - .05 / 2) * ypred$se) # Lower confidence limit.
yupper <- plogis(ypred$fit + qnorm(1 - .05 / 2) * ypred$se) # Upper confidence limit.
# Calculate conservative prediction curve.
yconservative <- plogis(ypred$fit + qnorm(1 - .05) * ypred$se)
# Calculate y values for plot.
yplot <- plogis(ypred$fit)
# combine them into a data frame for plotting 
predictionDf <- data.frame(Exp = predRange$Exp, Prob = yplot, yupper = yupper, ylower = ylower)
# plot
th_dat <- data.frame(x = 500, y = .5,
            label = paste0("At ", round(set4_ths[1], 0), " RFUs,\nthe estimated probability\nof dropout is 0.01."))
ggplot(data = predictionDf, aes(x = Exp, y = Prob)) + 
  geom_line() + 
  geom_ribbon(fill = "red", alpha = .4, aes(ymin = ylower, ymax = yupper))+
  geom_vline(xintercept = set4_ths[1], linetype = "dotted") + 
  geom_text(data = th_dat, inherit.aes = FALSE, aes(x =x, y = y, label = label), hjust = 0) +
  xlim(c(0, 1500)) + 
  labs(x = "Peak Height (RFUs)", y = "Probability of allele drop-out")
```

We can also look at a heat map of dropout for each marker by sample. All the loci in sample BC10.11 dropped-out, while most other samples have no dropout whatsoever.

```{r heatmap, out.width='100%', fig.cap="Dropout for all samples in `set4` by marker."}
set4_do %>% tidyr::separate(Sample.Name, into = c("num", "name", "num2")) %>% 
  mutate(Sample.Name = paste(name, num, ifelse(is.na(num2), "", num2), sep = ".")) %>% 
  ggplot(aes(x = Sample.Name, y = Marker, fill = as.factor(Dropout))) + 
  geom_tile(color = "white") + 
  scale_fill_brewer(name = "Dropout", palette = "Set2", labels = c("none", "allele", "locus")) + 
  theme(axis.text.x = element_text(size = rel(.8), angle = 270, hjust = 0, vjust = .5), 
        legend.position = "top")
```


<a rel="license" href="http://creativecommons.org/licenses/by-nc-nd/3.0/us/"><img alt="Creative Commons License" style="border-width:0; display: block; margin-left: auto; margin-right: auto;" src="https://i.creativecommons.org/l/by-nc-nd/3.0/us/88x31.png" /></a>
