If you don't remember your password, you can reset it by entering your email address and clicking the Reset Password button. You will then receive an email that contains a secure link for resetting your password
If the address matches a valid account an email will be sent to __email__ with instructions for resetting your password
Forensic Genetics Research Group, Department of Forensic Sciences, Oslo University Hospital, NorwayDepartment of Forensic Medicine, Institute of Clinical Medicine, University of Oslo, Norway
Forensic Genetics Research Group, Department of Forensic Sciences, Oslo University Hospital, NorwayDepartment of Forensic Medicine, Institute of Clinical Medicine, University of Oslo, Norway
MPSproto requires calibration of stutters and noise.
•
MPSproto increases the sensitivity of MPS-STR mixture analyses.
•
MPSproto is a good alternative to EuroForMix.
•
The gamma model in MPSproto performed best overall.
Abstract
We have developed MPSproto as an extension of EuroForMix to improve handling of stutter artefacts and other typing errors that commonly occur in MPS-STR data. MPSproto implements two models for read depth: gamma and negative binomial. It differs from EuroForMix in that calibration is required before mixtures are interpreted. In this study a mixture dataset (2–4 persons) was revisited, where EuroForMix interpretation of MPS-STR mixtures using the LUS+ format was first described; the performance of this model was compared to the MPSproto models. Results indicated that, overall, the MPSproto models performed better than the conventional EuroForMix model, and the gamma model implemented in MPSproto performed best. Differences were highlighted and further investigated to establish causality. Goodness of fit tests showed that the MPSproto models were adequate for the sequence reads when a low analytical threshold was applied.
] for an overview of different software) are important for the interpretation of mixtures of DNA from multiple contributors. However, all models of existing software packages were designed to use data generated with capillary electrophoresis (CE) based technology. During the last decade a new technology known as Massively Parallel Sequencing (MPS) has gained attention and is already viewed as a viable alternative to CE by the forensic community [
]. Whereas CE provides information based on the number of short tandem repeats (STR), accompanied by a peak height (RFU) measurement, MPS also provides the actual underlying sequence of the repeats, which is accompanied by the number of reads (read depth). For instance, a length-based homozygote genotype of ‘10’/‘10’ could be composed of two different sequences that are isometric but heterozygous, e.g., ‘[ATCG]10’/‘[ATCG]6 ATCT [ATCG]3’. The detection of increased variability via sequencing, along with the potential to examine many more autosomal markers, leads to increased discriminatory power compared to the traditional CE approach.
However, if the DNA amount is low the discriminatory power may not necessarily be increased; in this scenario it would be more likely to observe an allele of a homozygous genotype than both alleles of a heterozygous genotype, especially if a high analytical threshold is used [
]. Hence to leverage the potential for increased discriminatory power with MPS, the analytical threshold must remain low. Yet, use of a low analytical threshold may increase the challenges of STR mixture interpretation – first, by the inclusion of more sequence errors (typically observed with MPS technology), and secondly, the alleles from a major contributor may produce stutter artifacts in the range of a minor donor, making it hard to distinguish between the donor alleles and stutter [
]. This latter issue is circumvented by developing suitable probabilistic genotyping models, where there is no need to definitively distinguish stutters and alleles [
]. For instance, the ‘n0’ stutter is composed of both backward and forward stutters in the same sequence strand, hence there is no net change in repeat length.
] established negative binomial models for read depth because reads are discrete counts. In contrast, EuroForMix is based on the gamma model which assumes a continuous outcome. Furthermore, Vilsen et al. [
] suggested that parameters for inter-locus balance and stutter models could be estimated based on a calibration dataset before attempting to interpret complex mixtures. As part of the stutter model, Vilsen et al. [
] used the concept of block length of missing motif (BLMM) to describe stutter size. For instance, if an allele is ‘[ATCG]6 ATCT [ATCG]3’, and a stutter product ‘[ATCG]5 ATCT [ATCG]3’ was produced from it, then the BLMM value of the stutter product is 6, which is also the longest uninterrupted stretch (LUS) value. If instead the stutter product ‘[ATCG]6 ATCT [ATCG]2’ was produced, the BLMM value would be 3, whereas the LUS value would still be 6. It is well established that the expected stutter proportion increase with the value of BLMM [
]. Here sequences are converted to a format of type “x_y” where x is the regular CE designation based on the number of repeats, and y is the longest uninterrupted sequence (LUS). EuroForMix uses this format to identify backward or forward stutters where the LUS information is utilized. This improves the discriminatory power when DNA quantities are high. However, for minor contributors evaluated in a prior study, the performance was similar to that achieved with ordinary CE nomenclature [
]. In that study, which was based on ForenSeq DNA Signature Prep Kit (Verogen, Inc.) data developed using the manufacturer’s protocol, the largest influence on performance was related to how the data were pre-filtered: use of a static analytical threshold with no stutter prefiltering performed best, since potential alleles from minor contributors were not removed. Also, that study considered an analytical threshold of 30 reads since EuroForMix had not been adapted to work with lower values. Hence, this provides a motivation to develop a model which simultaneously lowers the analytical threshold and applies a probabilistic model to evaluate stutters and potential sequence errors.
Multiple types of stutters have been observed for MPS-STRs: n-1, n-2, n+1, n0, and n-1 from a region other than the LUS [
]. The bracket format facilitates analyst visualization since the compact form is more easily read and interpreted. Recently, the lusSTR tool was developed to automatically convert sequences into this format [
] we used the bracket format to identify many of the afore-mentioned stutter types, and we quantified stutter proportion sizes using a beta regression with BLMM as an explanatory variable. In the present work we continue to use the bracket format in an adapted model for MPS-STRs and show how a calibrated stutter model, similar to [
], can be directly integrated into probabilistic genotyping.
This paper is structured as follows: The Materials and methods section introduces the dataset used for model calibration and mixture evaluation. Subsequently, the mathematical details of MPSproto are presented, and we describe how calibration of the model was carried out. Then a mixture evaluation study is described, in which several kinds of models were compared, implemented as part of either MPSproto or EuroForMix. The Results section then describes the outcomes of the comparison study, including model performance and goodness of fit.
2. Materials and methods
2.1 Data
MPS-STR data developed using the ForenSeq™ DNA Signature Prep Kit (Verogen, Inc.) with DNA Primer Mix B according to the manufacturer’s protocols were considered. The ForenSeq UAS version 1.3 (Verogen, Inc.) was used to obtain the sequences and read depths. This version has a hard-coded detection threshold, resulting in sequences with a minimum of 11 reads in the exported data.
A total of 30 undegraded samples were used for model calibration of MPSproto: each was run with template quantities of 2 ng, 1 ng and 0.5 ng, for a total of 90 single-source profiles. Libraries were pooled in sets of 32 for sequencing (30 single-source libraries plus a positive and negative amplification control per run).
The dataset used for the mixture evaluation study was the same as that used by Bleka et al. [
]. ForenSeq libraries were pooled in sets of 32 for MiSeq FGx sequencing (30 mixture libraries plus a positive and negative control per run). A list of the mixtures is given in Table S2 for ease of reference.
Allele frequencies were also based on the African American dataset [
], where the alleles were sequence-based. New allele observations were imputed with the minimum allele frequency 0.00146 (this was the minimum observed frequency), and frequencies were normalized afterwards (ensuring sum equal to 1). Allele frequencies are provided, together with the mixture dataset from http://www.euroformix.com/MPSproto. The coancestry coefficient (also known as the “theta-correction”) was set to 0.01 for all calculations.
All sequences were evaluated using a re-implementation of lusSTR, LUSstrR (v0.2.2) [
]. The “forward bracket format” was used for MPSproto and the LUS+ format was used for EuroForMix. MPSproto requires that sequences are converted into bracket format of type “” where each block consists of a motif which is repeated times, and the blocks are separated with a whitespace. If , then no brackets are used. For instance, for an allele sequence with =ATCG, =ATG, , and , the bracketed format of the allele would be “”.
2.2 The MPSproto model
2.2.1 Models for sequence reads
The MPSproto model was based on an extension of EuroForMix with marker efficiency parameters included, inspired by Vilsen et al. [
]. The implemented version we used, v0.8.1, supports two different distributions for the read depths: gamma (GA) and negative binomial (NB). GA assumes that read depths are continuous, as the model was originally designed for CE-based data (i.e., for peak heights). NB is a discrete distribution, based on read depth, and is proposed as a suitable choice for MPS data. We parameterized the NB model such that and had the same interpretation as for the GA model: the expectation and the coefficient of variation of a heterozygous read depth in a single-source profile without any degradation or stutter effects.
The models for signal intensity (read depth) for sequence at marker , all observed above analytical threshold (marker specific), are defined as follows:
-GA: gamma(shape=, scale=)
-NB: negative-binomial(mu=, size=)
where, without assuming any stutter model, we have:
(1)
for the GA model, and for the NB model we have:
(2)
Here the models assume that the evidence is comprised of K individuals (K-person mixture), where each individual contributes with mixture proportion parameter . is the allele contribution (0, 1 or 2) decided by the assumed genotypes vector : , where and is the indicator function. The marker efficiency parameters are restricted such that = 1. Formulae for the probability density and cumulative functions and of the read depths, and the derivation of the size argument are included in the Appendix B, Appendix B.
2.2.2 Stutter model
MPSproto extends the stutter model of EuroForMix by allowing a general set of stutter types and has already been described in [
]. We define as the expected stutter proportion of a particular stutter type coming from donor allele . Likewise, we define to be the stutter product coming from donor allele . Then, considering only stutter type , the shape or mu argument (GA or NB model) is modified to:
where indicates the identity without assuming any stutter model (Eq. 1 for the GA model or Eq. 2 for the NB model). The first term is the “donates” part whereas the second term is the “receives” part. When several stutter types () are defined, the formula extends to:
As an example, consider contributions of allele ‘4’, ‘5’, and ‘6’ (i.e., ) and assume that only backward (BW) and forward (FW) stutter types are used. The expected stutter proportions are estimated at , and , the BW stutter at allele , , and the FW stutter at allele , . The formula above hence gives:
In MPSproto we model the expected stutter proportion of allele for stutter type with a logit-link function:
(3)
where are model coefficients of the regression and is the BLMM for stutter and is the BLMM for a different motif (also for stutter ). The latter is used for n0 stutters where two parts of the sequence stutter simultaneously, but there is no net change of the size of the stutter product compared to its parent allele.
2.2.3 Noise model
Sequences in the calibration dataset which do not originate from any of the known contributors or defined stutters types are classified as noise. The noise model consists of two parts and is carried out separately for each marker: 1) the number of unique noise sequences and 2) the read depth for the corresponding noise sequences.
(A)
The number of (unique) noise sequences, , follows a geometric distribution with parameter :
(B)
The read depth of a noise sequence, , is proportional to a Type I Pareto distribution with parameter [
where is the analytical threshold for a particular marker. We included an upper limit cap of read depth for noise, such that the normalizing constant is calculated as
Importantly, the noise parameters, and , are marker specific and must be decided based on a “noise dataset”, i.e., neither contributor alleles nor stutter.
2.2.4 Dropout and degradation model
In MPSproto, dropout is modelled the same way as described for EuroForMix: the dropout probability of an unobserved allele is the integral from zero to for GA (continuous) or the sum from zero to for NB (discrete). As with EFM, the compound Q-allele (indicated as ‘99’ in the MPSproto program) was used to designate a dropout allele, representing all the non-observed alleles from the allele outcome (defined by the allele frequency information).
MPSproto also supports degradation modeling, where is scaled with , where is the degradation slope parameter and is the fragment length of allele at marker . To enable this, sequences in the evidence profile must be in the following format: “RU:Sequence”, where “RU” is the repeat unit (CE format), “Sequence” is the bracket format, and the two are separated by a colon. For instance, sequence [ATCG]10 must be described in the format “10:[ATCG]10”. The corresponding kit (e.g., “ForenSeq”) must also be specified, and the fragment length of each common RU must be defined. For the “ForenSeq” kit, the fragment lengths for each CE allele can be obtained from https://verogen.com/wp‐content/uploads/2018/08/ForenSeq‐DNA‐Prep‐Guide‐VD2018005‐A.pdf
2.2.5 Implementation
MPSproto is implemented as an R-package available from github (https://github.com/oyvble/MPSproto) or the webpage http://www.euroformix.com/MPSproto. Numerical tests of the likelihood functions (for both GA and NB models) are included in the package, where calculations from the C++ code are compared to calculations from the R code.
2.3 Calibration of MPSproto
The model calibration of MPSproto consists of three main steps:
1)
Estimation of marker amplification efficiencies
2)
Identification of stutter types and estimation of regression coefficients
3)
Estimation of the noise model
The stutter and noise calibrations are conducted on a per-marker basis where all sequences are assumed to be above the defined analytical threshold, which can be specified differently for each marker. In this paper was used for all markers analysed with MPSproto.
2.3.1 Estimating marker amplification efficiency
To estimate the marker amplification efficiency parameters , the read depth sum per marker per sample is used. For each calibration sample , the sum of reads per marker is calculated as . If (no marker dropout), we assume that the sum is gamma distributed (for the GA model), with as the shape argument and as the scale argument. On the other hand, if (marker dropout is present), then the dropout probability is calculated as
where is the cumulative density function of the gamma distribution (applied with same shape and scale argument, see Appendix). The estimated marker amplification efficiency parameters are obtained by maximizing the likelihood
2.3.2 Stutter calibration
The following candidate stutter types are considered:
•
BW1: Backward stutter (n-1) where the longest repeat is stuttering.
•
FW1: Forward stutter (n+1) where the longest repeat is stuttering.
•
DBW1: Double backward stutter (n-2) where the longest repeat is stuttering.
•
BW2: Backward stutter (n-1) where a repeat other than the longest repeat is stuttering.
•
FWBW: Longest repeat produces forward stutter, and a repeat other than the longest repeat produces backward stutter (n0)
•
BWFW: Longest repeat produces backward stutter, and a repeat other than the longest repeat produces forward stutter (n0)
This stutter naming format is required by MPSproto and the analysis is carried out separately per marker. All stutter relations are identified based on the bracketed format: the internal function getStutteredSequence is applied on all donor alleles and compared with artefacts. Stutter artefacts are not considered in the calibration if multiple parental allele candidates are observed.
This is only true for the calibration model. All stutters are considered by the interpretation model, with no removal.
For instance, if artefact [ATCG]11 together with donor alleles [ATCG]10 and [ATCG]12 is observed, then the artefact could either be a forward stutter of [ATCG]10, a backward stutter of [ATCG]12, or both, and hence would be removed. Similarly, as carried out by [
], beta regression models were applied per stutter type to estimate the expected stutter proportion in Eq.3. Here, the slope parameter(s) (and ) are only included if estimated as non-negative and significantly different from zero (p-value < 0.05 applied).
For the calibration performed in the present study, a minimum of three stutter artefacts were required per stutter type. Artefacts which occurred fewer than three times were instead used to estimate the noise model (see Section 2.3.3).
2.3.3 Estimating noise model
Sequences that are not known donor alleles or are not recognized as stutter are assigned as noise. The procedure to estimate the noise model parameters ( follows, carried out for each marker separately:
(A)
Let the number of noise sequences per sample be given as . An additional observation is always appended as to make the model more robust for new observations. The parameter is then estimated as , where is the mean of the observations.
(B)
Let the read depth of the noise sequences be given vectorized as . An additional observation is always appended as to make the model more robust for new observations. If (no observations) two additional observations, and, are appended to enable the model to be fitted. The parameter is then estimated using the maximum likelihood estimate, obtained as .
2.3.4 Implementation and functionalities
The entire calibration procedure can be performed by running the function calibrateModel in MPSproto. This eases the calibration steps for users. Additionally, a tutorial is available at https://github.com/oyvble/MPSproto/blob/master/doc/MPSproto_tutorial.html.
2.4 Inference of MPS evidence and LR calculations
2.4.1 Definition of the likelihood function and the calculation of LR
For a new evidence profile, the unknown model parameters are =, whereas the stutter model parameters, marker amplification efficiency and noise model are based on the calibration. As with EuroForMix, the LR calculation is based on optimizing the likelihood function over the model parameters for a specific proposition :
where for a given marker , is the evidence information (sequences and read depths), and is the joint genotype combination of contributors, traversed over the joint genotype outcome defined from proposition . More specifically:
where is the assumed number of noise sequences for a given genotype combination (not derived from a contributor or a stutter). If a sequence is assumed to be a noise sequence, then . The calculation of the genotype probabilities, , includes the coancestry coefficient to take sub-population structure into account, with formulae described in [
]. Finally, the LR for evaluating the evidence given two hypotheses, and , is calculated using the maximum likelihood approach:
2.4.2 Mixture evaluation study
LRs for the 2 and 3-person mixtures were calculated by comparing the following propositions, where the person of interest (POI) is the hypothesized contributor:
-
: “POI + () unknown unrelated are contributors”
-
: “ unknown unrelated are contributors”
where is specified as the number of known contributors to the constructed mixture. For 4-person mixtures, the major contributor (always designated as ‘Ref1’) was conditioned under both propositions such that the hypotheses were modified to:
-
: “POI + Ref1 + 2 unknown unrelated are contributors”
-
: “Ref1 + 3 unknown unrelated are contributors”
LR calculations were performed using each of the following models:
1.
MPSproto applied with gamma model (GA)
2.
MPSproto applied with negative binomial model (NB)
3.
EuroForMix (EFM)
For MPSproto interpretations, a calibrated model was first estimated from the single-source data described in Section 2.1 using the steps described in Section 2.3. For LR calculations an analytical threshold of T = 11 was applied to match the value used for the calibration. The degradation model was not used. A tutorial describing how to use MPSproto to calculate the LR can be found at https://github.com/oyvble/MPSproto/blob/master/doc/MPSproto_tutorial_advanced.html.
The EuroForMix evaluation originally performed in Bleka et al. [
] was used, where the drop-in probability equals 0.05 and the hyperparameter lambda equals 0.01. The degradation model was always applied. No stutter filters were applied, and both forward and backward stutter models in EuroForMix utilized the LUS+ format. If EuroForMix failed to converge for any of the propositions, the model configuration was reduced to calculate with the backward stutter model only; if this also failed, the model was further reduced to apply no stutter model. A static analytical threshold of T = 30 was mainly applied since this was also used in [
]. However, we also investigated the effect when lowering the threshold to T = 20 or T = 11.
LR calculations were carried out for all true contributors as the POI ( true), except for the 4-person mixtures where Ref1 was a known contributor under both propositions as describe above. This comprised a total 160 true comparisons. Examples where models obtained very different LR values were inspected in detail: the LR per marker and model fits were further investigated.
Non-contributor tests ( true) were conducted for additional reference profiles (‘Ref5’,…, ‘Ref10’). Ref3 and Ref4 were also non-contributors for the 2-person mixtures, and Ref4 was also a non-contributor for the 3-person mixtures. This comprised a total of 420 non-contributor comparisons.
The following maximum likelihood function optimizer parameters were used for both MPSproto and EuroForMix: Seed= 1234; number of required optimizations (nDone)= 3; optimization accuracy (steptol)= 0.001.
2.4.3 Model comparisons
For the model comparison we plotted the pairwise LR values for each model and highlighted differences, using Receiver Operating Characteristics (ROC) curves [
]. For a certain threshold , a true positive event is when when is true. If instead was true, the event would be a false positive. The true positive rate is defined as the proportion of true positive events out of all true events, whereas the false positive rate is defined as the proportion of false positive events out of all true events. We also calculated the area under the ROC curve (AUC) for each model using the pROC R-package, including the 95 % confidence interval (function delong).
Goodness of fit for the fitted models was quantified using the Kolmogorov-Smirnov test by comparing an empirical cumulative density function against the cumulative uniform distribution. This was done for each optimized likelihood ( and separately) of each comparison and model. The empirical cumulative values were obtained by sorting the cumulative probabilities obtained from the R-package function validMLE (MPSproto) or validMLEmodels (EuroForMix). A small p-value (e.g. <0.01) from the test indicates a non-adequate model for the read depth.
3. Results
3.1 Calibration of MPSproto for the mixture evaluation study
All 90 single-source profiles were used to estimate the marker amplification efficiency parameters (Table 1). Only five marker dropouts were registered.
Table 1Gives a summary of the parameter estimates from the MPSproto calibration. The estimated parameters varied across the markers, illustrating that assignment of universal parameters would not be ideal.
Marker
A
p
lambda
BW1
FW1
DBW1
BW2
FWBW
BWFW
CSF1PO
0.324
0.989
15.6
-4.44/0.133
NA
NA
NA
NA
NA
D10S1248
0.629
0.91
3.85
-3.91/0.114
NA
-3.83
NA
NA
NA
D12S391
0.574
0.623
3.13
-3.23/0.123
NA
-5.47/0.139
-4.99/0.227
-3.78
NA
D13S317
0.822
0.85
4.08
-5.73/0.201
-4.02
NA
NA
NA
NA
D16S539
1.14
0.765
3.13
-3.79/0.147
-4.4
-7.06/0.246
NA
NA
NA
D17S1301
0.543
0.989
15.6
-4.25/0.2
NA
-3.87
NA
NA
NA
D18S51
0.968
0.778
2.41
-3.69/0.101
-6.74/0.16
-6.02/0.11
NA
NA
NA
D19S433
1.12
0.204
3.55
-3.74/0.102
NA
NA
NA
NA
NA
D1S1656
0.305
0.867
4.26
-4.3/0.17
NA
NA
NA
NA
NA
D20S482
2.44
0.867
3.12
-2.13
-3.97
-4.87
NA
NA
NA
D21S11
0.533
0.746
3.87
-3.54/0.0843
-3.61
NA
-3.36
NA
NA
D22S1045
1.44
0.948
6.04
-3.6/0.118
-5.57/0.146
-6.2/0.136
NA
NA
NA
D2S1338
1.21
0.446
1.85
-2.63/0.0505
NA
-5.37/0.0982
-6.92/0.427
-4.38
-4.37
D2S441
1.14
0.901
3.39
-3.56
-4.9
NA
NA
NA
-4.31
D3S1358
2.08
0.286
2.43
-4.15/0.127
-6.85/0.177
-7.05/0.178
NA
-4.69
-4.83
D4S2408
0.993
0.892
4.01
-5.27/0.208
-6.86/0.23
NA
NA
NA
NA
D5S818
0.29
0.968
12.4
-2.85
NA
NA
NA
NA
NA
D6S1043
1.09
0.85
3.03
-3.67/0.1
-7.42/0.259
-4.59
NA
NA
NA
D7S820
0.791
0.467
2.24
-5.17/0.222
-4.35
NA
NA
NA
NA
D8S1179
1.48
0.875
3.21
-3.39/0.122
-6.83/0.217
-6.3/0.178
NA
NA
NA
D9S1122
1.55
0.812
3.46
-4.77/0.207
-4.19
-7.53/0.231
NA
NA
NA
FGA
1.05
0.655
3.56
-3.03/0.0873
-5.78/0.112
-5.79/0.135
NA
NA
NA
PENTA E
0.256
0.968
6.17
-2.5
NA
NA
NA
NA
NA
TH01
2.77
0.728
2.97
-4.19/0.193
-4.81
-5
-5.25
NA
NA
TPOX
0.894
0.938
3.81
-5/0.167
NA
NA
NA
NA
NA
VWA
0.382
0.958
3.09
-3.56/0.12
NA
NA
NA
NA
NA
PENTA D
0.181
0.431
2.83
-5.7/0.189
NA
NA
NA
NA
NA
Table 1: Estimated parameters used for the calibrated MPSproto model. A is the estimated marker amplification efficiency, p is geometrical model parameter for the number of noise alleles and lambda is pareto model parameter for the noise size. The remaining columns indicate the stutter coefficients (intercept and slope) of the different stutter types separated by a backslash. Slope was only included if significant. “NA” means that the stutter type was not modelled.
In total, 3240 sequences unambiguously identified as stutter (as defined in Section 2.3.2) were used to estimate stutter regression coefficients (Table 1). In addition, 1181 sequences were identified as “noise” and were used to estimate the noise model parameters (Table 1). Eleven of these were originally denoted as stutter, but, as noted in Section 2.3.2, were re-categorized as noise since there were too few observations (< 3) for a given stutter type and marker. The vast majority of the noise sequences (1163, or 98.5 %) differed by a single base compared to one of the parent alleles. The number of observations per marker and per stutter/noise type is given in Table S1 (supplementary).
3.2 LR for true comparisons
LR values obtained using the different models were compared in terms of “ban” unit (for instance, “3 bans” is same as log10LR=3). Results of all true calculations can be found in Supplementary Table S3. Overall, the GA model returned much greater LR values than the EFM (T = 30) model (Fig. 1A). There were a few exceptions (datapoints to the right of the line of equality in Fig. 1A), but only one with more than 2 bans in difference: Comparison Ref3 to mixture 4Q_0.067ng_10‐2‐2‐1 gave 3.6 bans for EFM but only 0.6/− 0.11 for GA/NB. By contrast, there were 56 instances where the GA model resulted in LRs more than 6 bans greater than for EFM (see Supplemental Table S4). This magnitude of difference was observed across the entirety of the LR range of EFM. For instance, for Ref3 compared to mixture 3Q_1.35ng_20‐5‐1, EFM obtained 0.2 bans, whereas an LR of more than 15 bans were observed for GA. The cause of this difference was mainly due to a high number of dropouts: 14 allele dropouts for EFM, but only 3 for GA. Several of the other comparisons which returned a difference greater than 15 bans also had a large difference in the number of dropouts. Details are highlighted in Appendix B, Appendix B.
Fig. 1LR comparisons between different models: (A) comparison of GA versus EFM, and (B) comparison of GA versus NB. The colored points are Hp true comparisons whereas the black points are Hd true comparisons. The different color levels indicate the amount of DNA (ng) for the POI. The shape of the points indicates the number of contributors in the mixture: 2,3 and 4-person mixtures are indicated as squares, circles and triangles, respectively. Data points with values outside the figure range are indicated with a cross. The diagonal whole line is the identity line (y = x), whereas the dotted diagonal lines are ± 2 and ± 6 shifts of y-intercept. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
Comparing the results obtained with the GA and NB models, it was observed that LRs for GA were almost always larger than for NB, however the values were quite similar for most comparisons (Fig. 1B): 76 % returned a difference of less than 2 bans, whereas for 20 % the GA model produced LRs 2–6 bans higher than NB. There were only two instances where the LR using the NB model was more than 2 bans higher than for GA: Ref2 compared to 2Q_0.075ng_2‐1 returned 13.7 and 10.9 bans respectively, and Ref2 compared to 4Q_0.67ng_2‐1‐1‐1 returned 13.2 and 9.4 bans. Conversely, there were four instances in which the GA model resulted in an LR more than 6 bans higher than that obtained with the NB model (see Table S5 and Supplemental Section 1 for details). A comparison of the LRs obtained using the NB model and EFM is illustrated in Fig S1 in the supplementary material, showing lesser LR differences compared to the GA vs EFM comparison.
LRs were large (nearly always > 10 bans) for all three models when there was no allelic dropout for the POI (Fig. 2). However, they were drastically reduced when there was at least one allele dropout, reflecting a reduced quantity of DNA and loss of information. At any given level of allelic dropout the LRs produced by the different models were similar, although there was a tendency for EFM to return larger LR values when at least seven dropouts were present. Overall, the GA model typically returned larger LR values compared to the NB model across the full dropout range.
Fig. 2LR as a function of allele mismatches (dropout) for true contributors. The occurrence counts of each model are provided below each bar. The numbers of mismatches are grouped together in intervals, denoted by the numbers in square brackets below the corresponding group. The horizontal lines indicate 0, − 1 and − 2 bans.
As Fig. 2 indicates, the LR trended towards 1 as the number of dropouts increased and was sometimes less than 1. The fewest number of dropouts that resulted in an LR below 1 was inspected for each model: this was nine for the NB model (Ref3 compared to 3Q_0.34ng_20-5-1 gave log10LR=−0.27), 20 for the EFM model (Ref4 compared to 4P_0.67ng_20-5-5-1 gave log10LR=−0.04) and 22 for the GA model (Ref2 compared to 4P_0.067ng_10-2-2-1 gave log10LR=−1.148). A list of all instances which the LR was less than 1 for any of the models can be found in Table S6 in the Supplementary. We also created the same figure for non-contributors (Fig S2 in Supplementary). Here the lowest number of mismatches that were obtained for non-contributors with LR above 1 was found to be 18 for EFM (T = 30) and 10 for both MPSproto models. The GA model gave consistently larger LR values for the non-contributors compared to the NB model across the full dropout range.
3.3 ROC analysis and LR for true comparisons
The ROC analysis indicated a slightly better AUC value for the NB model compared to the other two models (Fig. 3). The MPSproto models (GA and NB) resulted in higher true positive rate compared to EFM (T = 30) for false positive rate values up to approximately 0.1. For this range, the GA model returned higher true positive rates than the NB model. However, for higher false positive rates, the GA model returned fewer true positive rates than both the NB and EFM (T = 30) model, which behaved very similarly to each other.
Fig. 3The figure shows the true positive and false positive rates for the 160 and 420 comparisons. Each model follows separate colored lines, whereas the symbols indicate specific LR thresholds. The AUC, with corresponding 95 % CI in brackets, are included in the legend for the corresponding method. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
With an LR threshold equal to 1, the number of false negatives/false positives for the three models, GA, NB and EFM (T = 30) were 7/66, 10/36 and 10/39, respectively (Table 2). When the LR threshold was increased to 10, the numbers for the three models were 9/0, 15/3 and 21/0, respectively. Consequently, LRs were low when false positives occurred (Fig. 1): For the NB model three circumstances were obtained where the LR was greater than 10 (references Ref4, Ref7 and Ref10 compared to sample 3P_0.067ng_20-5-1), with largest being LR= 63.4 (Ref7). All remaining non-contributor LR values were below 10 for all models. A list of all instances of non-contributor LRs above log10LR= 0.5, for all models, can be found in Table S7 in the Supplementary material.
Table 2Gives an overview over the true positive rate (TPR), false positive rate (FPR), the number of false negatives (FN) and the number of false positives (FP) for different LR thresholds (Thresh) and models. The NB and GA models both used analytical threshold T = 11.
The Kolmogorov-Smirnov p-values between for the two MPSproto models were uniformly distributed (see marginal distributions in the margin of the figures), indicating adequate models (Fig S3). Also, the p-values seemed independent (no pattern observed), indicating different model properties. On the other hand, the p-values for EFM (T = 30) clustered closer to 0, indicating a poor model fit for most of the samples (Fig S4).
3.5 Reducing the analytical threshold for EuroForMix
Reducing the analytical threshold for EFM further down to T = 20 or T = 11 (the latter value used for both MPSproto models) also reduced the LR differences compared to the MPSproto models, especially for comparisons in the range of log10LR < 10 (Fig S5 in supplementary). However, the reduction of the EFM analytical threshold in general caused a further reduction in the KS p-values, and hence even worse model fits (Fig S6 in supplementary). Nevertheless, the reduced thresholds of T = 20 increased the sensitivity of the EFM model since fewer false negatives were obtained (Table 2). The AUC was also highest for this model (Fig S7 in supplementary). However, reducing the threshold further to T = 11 reduced the sensitivity. We found that when the LR threshold was increased to 10, the false positives vanished for all almost all models (a few were obtained for the NB model). For this LR threshold, the number of false negatives was only nine for the GA model, followed by EFM with a lowered threshold: 14 for T = 20 and 15 for T = 11.
Although most of the LR values increased using EFM when applying T = 20 over T = 30, there was one example were this was not the case: Ref2 compared with 2Q_1.5ng_40-1 recorded a reduction of approximately 6 bans (log10LR=−6.7 instead of log10LR=0.1 obtained with T = 30). Here, around 9 bans were obtained for the MPSproto models. The reason for the large reduction was due to artefacts that could not be modelled by EFM; instead, they were therefore modelled as originating from an unknown contributor. For this example, model validation failed for EFM, but not for MPSproto (GA): the KS p-value under was 0.004/5e-06 for EFM T = 30/T = 20.
4. Discussion
4.1 Summary and overall findings
In this study the development and testing of a new probabilistic genotyping software, MPSproto, is described to evaluate MPS-STR mixtures. MPSproto is an extension of EuroForMix, and was inspired by the work of Vilsen et al. [
], who included marker amplification efficiency parameters to handle the large inter-locus balance that has been observed for MPS-STRs. MPSproto was already described and applied to an illustrative example in the discussion of Agudo et al. [
]; who focused on developing a framework to evaluate complex sequence stutters, describing their sizes using the bracket format nomenclature. The software includes two distribution options for the sequence read depths: the gamma (GA) model which is continuous, and the negative binomial (NB) model which is discrete.
In a departure from the methods used by EuroForMix, MPSproto requires calibration prior to implementation and evaluation of questioned DNA-profiles. The calibration dataset must be comprised of single-source profiles from donors whose alleles are known, and ideally should not include any degraded profiles, since the purpose of calibration is to characterize inter-locus variation. To simplify MPSproto calibration for users, the calibrateModel function in the software is used to perform all necessary calibration steps.
To evaluate the performance of the two MPSproto models compared to EuroForMix, the analysis described by Bleka et al. [
] was revisited for EuroForMix interpretation of 60 autosomal STR mixtures sequenced using the ForenSeq DNA Signature Prep Kit. Different analytical thresholds were applied: T = 11 for MPSproto versus T = 30 for EuroForMix, since for the latter, this was the threshold used in [
] and EuroForMix has not been adapted to work with lower values. However, we investigated the effect of reducing thresholds applied to EuroForMix: T = 20 and T = 11.
With this dataset, the MPSproto GA model performed considerably better than EuroForMix (T = 30), mainly because fewer dropouts were observed due to the lower analytical threshold applied to the former model. However, lowering the EuroForMix threshold reduced the differences, especially when T = 20 was applied instead. However, a worse model fit was returned for EuroForMix, since this also increased the number of artefacts that could not be properly modelled.
Overall, the GA model performed better than the NB model, since higher LRs for true contributors were obtained for the former. Additionally, the former model obtained a higher true positive rate for false positive rate values up to approximately 0.1. There were some situations where the GA model far outperformed the NB model. For example, when Ref2 was interrogated as a contributor to mixture 2P_0.75ng_20–1, the LR resulting from use of the GA model was log10LR= 10.7 versus log10LR= 4.0 for the NB model. With this mixture, Ref2 was a minor contributor with 11 dropouts, and the LR difference between the two models was largest for the D2S441 marker where Ref2 had one allele dropout. In the evaluation of this mixture, a lower penalty for dropout was applied when the GA model was used compared to that applied with the NB model. With further analysis we estimated the dropout probability of the Ref2 allele to be approximately 21 % for the GA model and only 2 % for the NB model (see Section 4.2 for more details). This indicates that the GA model is less strict and would be preferable to evaluate minor contributors compared to the NB model when an analytical threshold of 11 reads is applied (see Section 4.2 for further discussion of thresholds). However, the consequence is that non-contributors would be “excluded” to a lesser degree (albeit with low probative LRs).
Findings regarding model performance in this study do not necessarily support the suggestion by Vilsen et al. [
] that the NB distribution provides a better model than GA when the DNA quantities are low: The GA model seems to provide a higher dropout probability than for the NB model. However, it is possible that our conclusion in this respect could be different if lower analytical thresholds were applied to the mixture dataset; indeed, the effect of applying alternative analytical thresholds was not explored here. Regardless, there may be situations where it could be better to use the NB model – for instance, if the model fit using GA was poor, and the model fit diagnostics indicated that NB was better. For this reason, it may be worthwhile to compute the LR for both models when minor contributors are evaluated as the POI, followed by an examination of the model diagnostics to select the best fit model.
In the present study, LR values for non-contributors mostly did not exceed approximately 10, regardless of the program or model employed for interpretation (there were three observations above LR=10 for the NB model). All LR values greater than 1 were from individuals unrelated to the mixture donors (no check of relatedness amongst non-contributors was performed). Estimating the number of contributors to an unknown sample can be challenging and can impact interpretation results – typically in the form of false negatives for true donors when the contributor number is underestimated, and an increased number of false positives for non-contributors when the contributor number is overestimated [
]. In this study the ground truth number of contributors to the constructed mixtures was used for all evaluations (MPSproto and EuroForMix); with donor DNA inputs as low as 2 pg (see Supplemental Table S2), the contributor number used may therefore have been higher than the apparent number of contributors due to high levels of allelic dropout.
4.2 Model considerations
4.2.1 Stutters
As with EuroForMix, MPSproto uses stutter proportions instead of stutter ratios. The difference between the programs is in the way that stutter proportions are assigned: EuroForMix estimates stutter proportion parameter(s) when the questioned profile is evaluated, whereas MPSproto requires a calibration of the parameter in advance using a single-source dataset. The calibration of the stutter model for MPSproto is somewhat challenging, since it requires 1) the donor alleles to be known, 2) heterozygous allele pairs with sufficient separation (to not mask stutter effects), and 3) a sufficient number of stutter products (per sequence and across different alleles).
Rather than modelling stutter proportions directly, it may also be possible to take a fitted linear model of a stutter ratio, and convert it to a stutter proportion, . Vilsen et al. [
] suggested model for the expected stutter ratio. Conversion back to the expected stutter proportion could be accomplished using the transformation .
4.2.2 Noise
As with stutter, MPSproto requires calibration of the noise parameters in advance using a single-source dataset. MPSproto models the number of noise sequences per marker using a geometrical distribution, and the noise read depth using a discretized pareto distribution. These models provided a good fit for the calibration dataset used in this study. The pareto distribution is a useful choice since it can consider higher noise levels than for instance the exponential distribution, due to heavier tails.
Importantly, as the noise parameter developed during calibration is dependent on the analytical threshold selected, the analytical threshold used for evaluation should be the same as that used for the calibration.
The noise model may be affected by the number of libraries multiplexed for sequencing. For example, the expectation is that a sequencing pool composed of 90 libraries will result in reduced quantity of noise signals/sequences at reduced read depths as compared to a sequencing a pool of 30 libraries. Therefore, a calibration may be needed for each distinct protocol used for data generation. In the calibration dataset we observed that most of the noise sequences (98.5%) appeared to be errors of a single base compared to a parental allele sequence. Hence, there is potential to model the single base error sequencies separately from other type of noise, which may help remedy the above-mentioned issue.
] proposed a one-inflated (zero-truncated) negative binomial distribution for the noise read depth, including all sequences down to singletons (i.e., analytical threshold T = 1). Our approach to handling noise sequences differs: we recommend avoiding an analytical threshold that is too low because of the way that the MPSproto noise model is defined (it may no longer fit). An alternative approach is to apply corrections: either of the type employed by FDStools [
‘FDSTools: a software package for analysis of massively parallel sequencing data with the ability to recognise and correct STR stutter and other PCR or sequencing noise’.
Probabilistic genotyping (PG) systems are well established for CE based applications to deal with stutters. However, interpretation for MPS-STR data lags because of the current lack of availability of such systems, and accordingly has relied on defining thresholds to facilitate interpretation. Experience with PG systems applied to CE has shown that evaluation of evidence that is based upon thresholds, where binary decisions are required to designate alleles, are always suboptimal [
]. This is because information is lost or thrown away. This in turn results in reports that either understate the value of the evidence when true alleles that support are removed (i.e., dropouts for the person of interest); or overstate the value of the evidence when true alleles that support are removed (i.e., dropouts for the unknown contributors). Given the increased complexity, it can be argued that the inherent advantages of MPS will never be properly realized without integration of PG systems for evidence evaluation. Use of a threshold intended to eliminate noise cannot be avoided with the current implementation of MPSproto, but steps have been taken to model noise, reducing the threshold to only 11 reads – which is substantially lower than that used in conventional systems.
4.2.4 Marker amplification efficiencies
On the basis of prior findings regarding variability in marker amplification efficiencies with the ForenSeq DNA Signature Prep Kit [
], MPSproto includes a function to adjust the marker amplification efficiency estimates towards those obtained from a questioned DNA profile by utilizing a prior. To enable this, standard deviations of the amplification efficiencies must be specified in addition to the mean (per marker), and the user can choose between normal distribution or log-normal distribution for the prior. We did not use this functionality in this study, although it would be worthwhile for future investigations, where the prior is, for example, determined from the calibration dataset.
4.2.5 Degradation
Similar to EuroForMix, the MPSproto model also includes the possibility to consider a degradation model for evaluation of questioned DNA profiles. In this study, both the calibration and mixture datasets were constructed without degradation exhibited. When a degradation model was applied in MPSproto, the degradation slope estimates were all close to one (not shown), confirming the pristine state of the test samples. It is highly recommended that the calibration dataset exhibits little or no degradation, so that the estimated marker efficiencies are not affected before the full model is applied to questioned DNA profiles (which may be highly degraded).
EuroForMix was used with the degradation model turned on since more of the data is explained with rather than without using it. The underlying reason is that the marker amplification efficiencies of samples in the mixture dataset decay with the fragment length.
4.2.6 Data format
The format required for MPSproto interpretation of MPS-STR sequences does not necessitate the use of any specific alignment or data analysis software. Any program that can be used to perform marker identification and sequence allele calling from FASTQ files (such as STRaitRazor [
‘FDSTools: a software package for analysis of massively parallel sequencing data with the ability to recognise and correct STR stutter and other PCR or sequencing noise’.
]), along with bracketing of the sequences, could be used as part of a pipeline to produce the data for MPSproto input. In this paper, sequences were obtained from ForenSeq typing results using the ForenSeq UAS v1.3, and these were further converted into a bracket format (forward strand; using the lusSTR program [
Massively parallel sequencing of forensic STRs: considerations of the DNA commission of the International Society for Forensic Genetics (ISFG) on minimal nomenclature requirements.
We revisited the interpretation where reference Ref2 was compared to sample 2P_0.75ng_20‐1 since this was one of the examples where the two MPSproto models differed most (log10LR=6.7 in difference). The largest LR difference was observed for D2S441 which was estimated as 1.14 in marker amplification efficiency (Table 1). One of the heterozygous alleles of Ref2 dropped out. Fig. 4 illustrates how the dropout estimate became larger for the GA model (dropout probability 21 %) than for the NB model (dropout probability equal 2 %), and hence was penalized less in the LR calculation. The reason for this is that the probability density function of GA is skewed towards zero whereas for NB it is not. The read depth distribution for the GA model is also wider than for the NB model, which means that the former also has a wider heterozygous balance distribution compared to the latter.
Fig. 4Comparison of distribution of read depths between the two MPSproto models. The estimated parameters were taken from the fitted model of sample 2P_0.75ng_20–1 where POI=Ref2 (minor). GA estimated Mx= 3.5 %, , whereas NB estimated Mx= 3.9 %, . We used marker efficiency parameter equal 1.14 (as estimated for D2S441; see Table 1). The vertical dashed black line is analytical threshold T = 11. The shaded gray area is part of the density function (continuous), whereas the vertical red lines are part of the probability mass function (discrete), used to calculate the dropout probabilities. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
] demonstrated that concurrent interpretation of PCR replicates in EuroForMix can produce improved results. We will examine the effects of using PCR replicates in a future study.
MPSproto can also be used to evaluate conventional CE-STR data (and accordingly, a CE module is included for the calibration step). Many stutter types are supported for this module: backward, forward, double-backward, double forward, triple backward and 2-bp stutters (half stutters). Here CE allele lengths are used to determine the BLMM. We will compare the performance between MPSproto and EuroForMix for the interpretation of CE data in a future study.
The current implementation MPSproto does not support combination of CE-STR with MPS-STR data for simultaneous evaluation. For this to be possible the model would need to be extended so that typing results for each marker can be differentiated according to CE or MPS origin, as a given marker could have observations from either CE typing, MPS typing, or both. For instance, the ForenSeq DNA Signature Prep Kit types 22 autosomal STRs that are also typed using the PowerPlex® Fusion 6C System Kit (Promega Corporation).
EuroForMix has been used to interpret SNP typing results [
] showed better performance with machine learning approaches). Accordingly, MPSproto can also be used to interpret SNP data. For this application, the MPSproto stutter model would be turned off and not defined, and the degradation model would not be used. However, calibration is still needed to define the marker amplification efficiency and noise parameters to be used. Provided STR and SNP markers are neither linked nor in linkage disequilibrium, the LRs for the two systems could be multiplied together. The performance of MPSproto with SNPs will be described in a future paper.
The current version of MPSproto does not have a graphical user interface (GUI). As most forensic practitioners are not used to running programs from the command line, a near-term priority is to develop a user-friendly GUI for both model calibration and profile interpretation. Additionally, as MPSproto requires specific formatting for the input sequence data, we aim to integrate the LUSstrR conversion tool [
] into MPSproto to ease translation of sequence strings into the proper format. Other MPS-STR assays such as the PowerSeq GY System (Promega Corporation), the Precision ID GlobalFiler NGS STR Panel (Thermo Fisher Scientific), or custom assays, could be used with MPSproto, provided that the sequences are converted to the proper format. However, to unlock the degradation model for these assays, the fragment lengths of common RU alleles must be defined in the MPSproto kit file; at present, only information about the ForenSeq DNA Signature Prep Kit is included.
Finally, there is no relationship testing module included in the current version of MPSproto, but we aim to include this in a future version.
5. Conclusion
In this paper we found that the developed MPSproto models are good alternatives to EuroForMix for the interpretation of MPS-STR mixtures typed using ForenSeq DNA Signature Prep Kit. MPSproto enables increased sensitivity by lowering the analytical threshold whilst demonstrating a characteristically good fit to the data. Although, it requires a model calibration of stutters, marker efficiency and noise before the interpretation.
Conflict of interest
The authors declare that they have no conflict of interest.
Acknowledgements
This work was funded under Agreement No. HSHQDC-15-C-00064 awarded to Battelle National Biodefense Institute (BNBI) by the Department of Homeland Security (DHS) Science and Technology Directorate (S&T) for the management and operation of the National Biodefense Analysis and Countermeasures Center (NBACC), a Federally Funded Research and Development Center. The views and conclusions contained in this document are those of the authors and should not be interpreted as necessarily representing the official policies, either expressed or implied, of DHS or the U.S. Government. DHS does not endorse any products or commercial services mentioned in this presentation. In no event shall DHS, BNBI or NBACC have any responsibility or liability for any use, misuse, inability to use, or reliance upon the information contained herein. In addition, no warranty of fitness for a particular purpose, merchantability, accuracy or adequacy is provided regarding the contents of this document. All research involving living individuals, their data, or their biospecimens was conducted in compliance with the Federal Policy for the Protection of Human Subjects (The Common Rule, codified for DHS as 6 CFR 46), DHS Management Directive 026-04, and any other applicable statutory requirements. Research involving human subjects has only been initiated after the following has occurred: the need for IRB review has been determined, IRB approval has been obtained as applicable, and DHS Compliance Assurance Program Office certification or concurrence has been issued. Notice: This manuscript has been authored by Battelle National Biodefense Institute, LLC under Contract No. HSHQDC-15-C-00064 with DHS. The US Government (USG) retains and the publisher, by accepting the article for publication, acknowledges that the USG retains a non-exclusive, paid up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for USG purposes.
Appendix A
Mathematical description of the distributions
Gamma distribution:
Continuous outcome: .
Probability density function:
Cumulative distribution function:
Negative binomial distribution (re-parameterized version):
Discrete outcome: .
Probability density function:
Cumulative distribution function:
Deriving the relation between size parameter and the coefficient of variation
We assume following parameterization: and .
Then the size parameter can be written as .
Since the coefficient of variation , then .
Computational sidenote: Very large values of can cause computational issues which is solvable by restricting it below a value . It can be derived that the restriction is equivalent to . Hence smaller values of would require a large value of . We used a value of in the C++ function which avoids the maximum likelihood optimization to crash using parallelization with OpenMP.
‘FDSTools: a software package for analysis of massively parallel sequencing data with the ability to recognise and correct STR stutter and other PCR or sequencing noise’.
Massively parallel sequencing of forensic STRs: considerations of the DNA commission of the International Society for Forensic Genetics (ISFG) on minimal nomenclature requirements.
00122-3&id=gr1.jpg){kind=link}
00122-3&id=gr2.jpg){kind=link}
00122-3&id=gr3.jpg){kind=link}
00122-3&id=gr4.jpg){kind=link}