Sage Journals: Discover world-class research

Abstract

The need to identify large numbers of missing people from mass disasters, armed conflicts, and human rights abuses has become increasingly prevalent. DNA methods for disaster victim identification (DVI) typically rely on the identification of postmortem (PM) samples through comparison to antemortem (AM) samples from close relatives through mitochondrial or short tandem repeat (STR) profiling. There are several limitations to these tests: STRs cannot identify relatives further out than second degree and do not always provide full profiles in degraded samples; and mitochondrial DNA requires a matrilineal relative. The ForenSeq^® Kintelligence Kit was designed to interrogate 10,230 forensically relevant single-nucleotide polymorphisms (SNPs) and allows up to three samples to be sequenced simultaneously. This sequencing plexity is required to type enough SNPs to detect relationships up to fifth degree when searching DNA databases such as GEDmatch PRO™, however DVI scenarios rarely require this level of kinship matching and need a higher throughput solution. In this article, we present a novel approach for the ForenSeq Kintelligence Kit, with an increased sequencing plexity of 12 PM samples and 32 AM samples. Mock PM samples tested included bones, dental remains, and artificially degraded and low input samples. An offline model of the ForenSeq Kintelligence kinship algorithm was used to compare profiles and determine kinship. Despite lower call rates and decreased heterozygosity, relationships could be determined out to the third degree with 100% specificity and 100% sensitivity. Using the ForenSeq Kintelligence workflow in this way can facilitate higher-throughput DNA analysis and still yield a high success rate to aid cases of DVI.

Introduction

Identification of victims of mass fatality incidents (MFI) is required for both humanitarian and legal reasons: giving families searching for missing family members resolution and justice for families when death is not accidental as in civil and criminal cases. This process can be challenging given the large numbers of victims and the impact of the disaster on the integrity of the bodies of the victims.

Global events such as the terrorist attack on the World Trade Center on September 11, 2001 or the boxing day tsunami on December 26, 2004 highlighted the need for efficient procedures and methods for recovering and cataloging the remains of victims, storing information related to the remains, and identification.^1–6

The short tandem repeats (STRs) have been used successfully for many years for human identification, including disaster victim identification (DVI),^1,7–11 primarily due to their highly polymorphic nature.^12–17 In highly degraded samples however, larger STR markers are less likely to be amplified with commonly used capillary electrophoresis (CE) based kits, thereby limiting the number of markers that can be typed to make an identification.

Although there have been considerable developments to reduce the size of the amplicons targeted with CE,¹⁸ this technology is limited by its reliance on size-based separation.^18,19 Next generation, or massively parallel sequencing (NGS or MPS) STR assays already accommodate smaller amplicon sizes and allow for the analysis of more markers within one assay.^20–25 Mitochondrial DNA can be sequenced when nuclear DNA is extremely degraded, but only the maternal lineage can be matched.

STR data is appropriate for cases where antemortem (AM) DNA samples are available from the missing person or very close family members such as first-degree relatives (parent, child, or sibling), but it is less successful if only more distant family members are available for comparisons, due to the number of false positive identifications that can occur with unrelated individuals.^26–29

Although single-nucleotide polymorphism (SNP) markers have less discriminatory power than STRs for individual identification, they have become an important tool in the context of DVI since the advent of MPS, given that many more markers can now be multiplexed.

The SNP amplicons are shorter than STR amplicons (usually <150 bp), making them well suited for the analysis of degraded DNA. They also have a lower mutation rate than STRs, which is useful for kinship analysis where decreased mutation rate over the course of generations will lead to a reduced chance of false exclusion.^30–36

The ForenSeq^® Kintelligence Kit was developed to interrogate 10,230 forensically relevant SNPs for the purpose of assisting in violent crime and cold case investigations, including missing persons identification.^37–41 Currently, up to three samples are analyzed per NGS run, which allows a sufficient number of SNPs typed for detecting relationships up to fifth degree when searching the GEDmatch PRO database to assist law enforcement in solving cold cases.^38,41,42

DVI requires a higher throughput solution due to the often high numbers of postmortem (PM) samples from an MFI, and the number of AM samples included for comparison. In addition to increasing throughput, a software solution is necessary that can: 1.

Identify relationships in the absence of a pedigree or information of a relationship, which can be important when trying to assign multiple remains to a single individual or when a relative of a victim is not available.

Reduce the false positive rate in identifying relationships, which is often elevated when calculating likelihood ratios (LRs) using STRs.²⁸

Maintain privacy for victims of MFIs, which currently requires knowledge of a programming language such as R, knowledge of an LR software such as Familias,⁴³ DNA View^©, or The Mass Fatality Identification System (M-FISys^©)^43–48; or a private account on Bonaparte⁴⁷ (can only identify up to second-degree relationships on a small set of markers). Alternatively, sample genotypes may be uploaded to a public database such as GEDmatch PRO or FamilyTreeDNA to query for relatives, which would require approval from relatives to upload their genotypes to a public database.

The Kintelligence kinship algorithm is currently only available on GEDmatch PRO, where users can upload Kintelligence GEDmatch PRO reports and query the database.⁴⁹ In this study, we utilized the Kintelligence kinship algorithm localized on a private computer to identify relationships among samples prepared with the Kintelligence kit, sequenced on the MiSeq FGx^®, and analyzed with the Universal Analysis Software (UAS).

We demonstrated that this algorithm can efficiently identify relationships, with high sensitivity and specificity, up to and including third degree for degraded/low input mock PM samples sequenced at a plexity of 12 and for mock reference or AM samples sequenced at a plexity of 32 for a high-throughput solution for DVI.

Materials and Methods

Human DNA samples and mock PM samples

NA24385 DNA from the Coriell Institute for Medical Research was used as positive control. Ninety-three mock AM samples were selected from the 1000 Genomes Project, CEPH collection, the Personal Genome Project, and Innogenomics, Inc. (New Orleans, LA) (Supplementary Table S1).

Thirteen mock PM samples consisted of contemporary tooth samples, bone samples, and ancient bone (Supplementary Table S2). Two series of degraded DNA were also purchased from Innogenomics, Inc. (Supplementary Table S2). Mock PM sample DNA extracts consisted of five contemporary tooth samples extracted with the dental forensic kit DFK^® (Supplementary Table S2; Innogenomics, Inc.⁵⁰), seven contemporary bone samples kindly provided by Dr. Rachel Houston (Department of Forensic Science, Sam Houston State University [SHSU], Huntsville, TX) in conjunction with the willed body program at the Southeast Texas Applied Forensic Science Facility (Supplementary Table S2⁵¹), and one DNA extract from an ancient bone of Eastern European origin provided by the laboratory of Dr. Mitchel Holland (Department of Biochemistry & Molecular Biology, Pennsylvania State University [PS], State College, PA) extracted with a total demineralization method (Supplementary Table S2^51–53).

SHSU bone DNAs were extracted using either the PrepFiler™ forensic DNA extraction kit (Thermo Fisher, Waltham, MA) for samples SHSU 1, SHSU 3, SHSU 4, SHSU 6, and SHSU 7, or demineralization protocol for bone samples SHSU 2 and SHSU 5.^52,53 The degradation index and DNA concentration of the SHSU bone DNA samples was determined using Quantifiler™ Trio DNA Quantification Kit (Thermo Fisher).

In addition, intentionally degraded DNA was used as either mock AM or mock PM samples. Two series of degraded DNA were purchased from Innogenomics, Inc.: DNA was extracted from whole blood using organic methods from two different male donors and sheared by sonication at 50°C for times ranging from 0 to 16 h (Samples 1231 and 3551, Supplementary Tables S1 and S2).⁵⁰

The InnoQuant HY real-time quantitative PCR (qPCR) kit (Innogenomics, Inc.) was used according to the manufacturer's instructions to assess the amount of human DNA and the degradation state.^50,54 Buccal samples were collected from volunteers from a family with a known pedigree (V004, V016-021), herein referred to as the Verogen Family (Supplementary Table S2), each of whom signed an informed consent form authorizing the use of de-identified samples for research use and publication.

DNA was extracted and purified from buccal swabs with the QIAmp DNA Investigator kit (Qiagen, Germantown, MD). Two DNA samples (V004 and V016) were artificially degraded using high temperature as follows: 5 replicates of purified buccal DNA from each individual were subjected to 21 cycles of heating and chilling at 98°C for 1 h followed by 4°C for 10 min.

DNA grade water (Fisher Bioreagents; Thermo Fisher) was added to the subsequently dried DNA to bring the DNA into solution. Degradation indices and DNA concentration were determined for all Verogen family DNA samples using the PowerQuant^® System (Promega, Madison, WI) following the manufacturer's instructions. Degradation indexes varied for replicates with values of 1, 2.1, 2.6, 5.1, and 20 for sample V004 and 1, 1.5, 2.0, 2.2, and 2.9 for sample V016 (Supplementary Table S2).

Before library preparation, intact DNA samples were quantified using the QuantiFluor ONE dsDNA system (Promega) for input into library preparation. Mock PM DNA samples were quantified utilizing the qPCR method. Unless otherwise noted, the DNA was diluted to 40 pg/μL for 1 ng total DNA added to the library preparation reaction (Supplementary Tables S1 and S2). The Positive Control DNA NA24385 was serially diluted to 20, 10, 4, and 2 pg/μL for total DNA inputs of 500, 250, 100, and 50 pg to mimic low input PM samples.

The purchased, artificially degraded samples had DNA concentrations sufficient to add 1 ng of DNA to the library preparation reactions (Supplementary Table S2). Not all the degraded Verogen Family samples had sufficient DNA concentration for input of 1 ng DNA into the library preparation reactions.

For sample V016, degraded replicates with degradation index (DI) of 2.0. 2.2, and 2.9 had sufficient DNA concentrations to add 1 ng to the library preparation reactions (Supplementary Table S2). The PS bone DNA concentration was estimated to be 390 pg based on the mtDNA quantification of ∼1400 mtDNA copies/μL (Personal communication; 4th November, 2020 [M. Holland]).

ForenSeq Kintelligence library preparation

ForenSeq Kintelligence libraries were prepared using the ForenSeq Kintelligence Kit (Verogen, San Diego, CA) following the manufacturer's instructions,⁵⁵ and libraries were quantified using the QuantiFluor ONE dsDNA system (Promega). The ForenSeq Kintelligence Kit supplies six unique dual indexed (UDIs) adapters.⁵⁶ To sequence libraries with plexities higher than six, additional UDIs were utilized (Verogen).

These additional UDIs are not supported in the UAS ForenSeq Kintelligence analysis pipeline, and so required the use of a custom informatics pipeline (refer to the Primary, Secondary, and Tertiary Data Analysis section). Degraded replicates of sample V004 with DIs of 2.1, 20, 5.1, and 2.6, had 600, 600, 700, and 250 pg total DNA input added to the library preparation reactions, respectively (Supplementary Table S2). To mimic low input samples, 50 pg of V004 with DI of 1 and 50 and 250 pg of V016 with DIs of 1 and 1.5, respectively, were added to the library preparation reactions (Supplementary Table S2).

Sample pooling and sequencing

After amplification of targets and purification, libraries were normalized to 0.75 ng/μL following the method described in the ForenSeq Kintelligence Kit Reference Guide.⁵⁵ If library yields were lower than 0.75 ng/μL, the library was pooled neat without dilution. Mock AM libraries generated from commercially obtained intact DNAs showed library yields greater than 0.75 ng/μL with three exceptions (Supplementary Table S3).

Some libraries generated from mock PM, low DNA input, and commercially degraded DNA samples also had yields greater than 0.75 ng/μL (Supplementary Table S3). Libraries were pooled by pipetting 8 μL of each library into a 1.7 mL microcentrifuge tube. Libraries generated from mock AM samples were pooled at sample plexities of 3, 12, 16, 24, 30, or 32 for denaturation and sequencing (Supplementary Table S3).

Libraries generated from mock PM DNA samples were pooled at sample plexities of 3 or 12 for denaturation and sequencing (Supplementary Table S3). Pooled libraries were denatured with freshly diluted NaOH (HP3) by incubation at room temperature for 5 min and then followed by a dilution with HT1 according to the ForenSeq Kintelligence Kit Reference Guide.⁵⁵

Human sequencing control, a library consisting of 33 STRs, 16 of which are analyzed, that serves as a positive sequencing control for the MiSeq FGx instrument,^49,55 was also similarly denatured and diluted with HT1 according to the ForenSeq Kintelligence Kit Reference Guide.⁵⁵ The denatured library pools were sequenced on the MiSeq FGx instrument with the MiSeq FGx Reagent Kit following the manufacturer's recommendations.^55,57

The higher plexity library pools containing UDI adapters not supported by the UAS were sequenced in Research Use Only mode on the MiSeq FGx.⁵⁷ Sequencing utilized 151 cycle paired-end reads for all libraries. The sequencing runs include 2 eight cycle indexing reads required to demultiplex the libraries utilizing the indices present in the UDI adapters.^56,58

To generate sequencing results for mock AM libraries with very high reads per sample for simulation studies, 30 libraries from mock AM samples (Supplementary Table S3) were sequenced on the NextSeq 500 instrument using NextSeq 500/550 High Output Kit v2.5 (300 Cycles) kit (Illumina, San Diego, CA) following the manufacturer's recommendations.^59,60

Primary, secondary, and tertiary data analysis

All library samples prepared with the six UDI adapters⁵⁶ supplied in the ForenSeq Kintelligence Kit sequenced on the MiSeq FGx were analyzed with the ForenSeq UAS v2.2 (Verogen) for allele and genotype calling as previously described.⁶¹

Data from library samples prepared with UDI adapters not supported in the UAS, the sequencing runs performed with the Research Use Only mode⁵⁷ were analyzed using the same bioinformatic pipeline utilized by the UAS but on a separate server through command line tools. First, samples are demultiplexed based on the supplied index sequences found on the UDI adapters by demultiplexing binary base call (BCL) files and generating FASTQ files.^58,62

Reads 1 and 2 are aligned to the primer sequences using the Smith-Waterman–Gotoh⁶³ algorithm. Reads aligned to specific primer pairs are assigned to the loci corresponding to those pairs. Alignments were then written in the binary alignment map format. At the position of each SNP, matches to the reference base call and matches to the alternate base call were counted, filtering on a minimum base quality of 30.

Number of reads were then summed for each type of call (reference or alternate) at each locus, requiring a minimum coverage of greater than 10 reads. The SNP genotypes were then determined by filtering on the analytical threshold (AT) and interpretation threshold (IT), which were both set to 3%.³⁸ The ATs and ITs were determined by multiplying 3% by the sum of read counts at that locus. When low coverage occurred, a minimum of 650 reads was used to calculate the threshold values. The resulting AT and IT values were then compared with the total read counts for the reference and alternate alleles for each locus. Genotypes were determined for each locus if the call passed both ATs and ITs. Genotyped results were then written in the variant call format.

Downsampling analysis

Thirty mock AM samples sequenced on the NextSeq sequencer (Illumina) (described in the Sample Pooling and Sequencing section) were demultiplexed using a custom build of the ForenSeq UAS secondary analysis pipeline. We calculated the number of reads to simulate a 16 sample Kintelligence run and a 30 sample Kintelligence run by dividing 25,000,000 total estimated number of reads sequenced on a MiSeq FGx by 16 or 30, respectively, then multiplying by the average percent aligned (96%).

This resulted in 1,500,000 reads for the 16-sample simulation and 800,000 reads for the 30-sample simulation. We randomly chose the desired number of reads (1.5 M for 16-sample and 800 K for 30-sample) from the FASTQ files, which is defined as downsampling, using seqtk (https://github.com/lh3/seqtk). Seqtk outputs new FASTQ files containing the desired number of reads, which were processed as described in the Primary, Secondary, and Tertiary Data Analysis section.

GEDMatch data simulations for whole genome kinship algorithm testing

GEDMatch database profiles were downloaded and analyzed as previously described.⁴¹ Briefly, a set of 1000 anonymized samples, termed “query samples,” were randomly selected from GEDMatch. These samples were then queried for relatives in the GEDMatch database using the standard microarray segment matching algorithm and any hits, termed “target samples,” were selected based on shared centiMorgans (cM) values calculated by the GEDMatch one-to-many tool.

This search resulted in 2954 target samples together with the matched query samples. These results included query-target pairs that had zero shared cMs, representing true unrelated pairs. Therefore, the number of target samples is larger than the number of query samples due to the inclusion of both related and unrelated target samples. Relationship degree was determined by comparing the resulting shared cM values generated from the one-to-many tool with the expected range of shared cM per degree provided by DNA Painter.⁶⁴

The loci typed for the samples included in the query and target set were first filtered for the 10,230 SNPs in the Kintelligence panel, then randomly filtered to 80%, 60%, 40%, and 20% call rates, resulting in 8000, 6000, 4000, and 2000 loci, respectively, that were called in each query-target pair. The whole genome kinship coefficient was calculated for each query-target sample pair for each level of reduced locus call rate using the Kintelligence kinship algorithm.

Pairs with a whole genome kinship coefficient greater than 0.031 were considered related; pairs with a whole genome kinship coefficient of ≤0.031 were considered unrelated. Sensitivity and specificity were calculated by comparing these results with the one-to-many tool query results. In other words, the one-to-many query results were considered the truth set and the results generated by the Kintelligence kinship algorithm were considered the test set.

Calculation of LRs and ForenSeq Kintelligence kinship values

The LR was calculated as follows:

where D represents the genotypes, H_r represents the hypothesis that the individuals are related, and H_u represents the hypothesis that the individuals are unrelated. The related hypothesis was signified by a pedigree where the unidentified individual was tested as the relative. The unrelated hypothesis was signified by a Hardy–Weinberg equilibrium calculation. An LR value was calculated per locus and then multiplied across loci, which resulted in a final LR for the relationship. To improve computational efficiency, each locus LR was converted to logarithm and loci LRs were summed.

However, due to the high plexity of this platform and the often high degradation level and/or low input of PM samples in MFI cases, stochastic effects during PCR cause allele drop out and can introduce a situation where the allele combination between two individuals at a locus appears impossible.

This typically occurs in parent-offspring relationships where each individual in the relationship is genotyped as homozygous for different alleles but one or both are actually heterozygous (e.g., parent has genotype AA and child as genotype BB when the parent or both parent and child are actually AB). Genotyping error and de novo mutation are also possible causes of such a case.

This situation results in an LR of 0 and can result in an overall LR of 0 given that the final LR is the product of the loci LRs. To ensure locus LRs are not 0, a modification was made to pedprobr. The LR for the locus in these cases was calculated as follows: $L R = \frac{0.001}{p^{2} q^{2}}$

where 0.001 represents a genotyping error rate, p is the allele frequency of the allele 1, and q is the allele frequency of allele 2.⁶⁵

This method is useful when the relationship between two individuals is unknown, thereby not requiring a pedigree. The PC-AiR method first takes a set of genotyped individuals and separates them into two nonoverlapping subsets: one set containing unrelated individuals that represent ancestries of all individuals (unrelated subset), the other set containing individuals that have at least one relative within the first subset (related subset).

To build the unrelated subset, we made a modification to the original PC-AiR method to improve computational efficiency in building the model. Samples with none or the fewest relatives are added to the unrelated subset, whereas those with more relatives are excluded from the unrelated subset.

This is performed by calculating a kinship value⁶⁶ per pair and categorizing each individual as having a relative or not based on stringent thresholds: a relative is considered if the kinship value is greater than 0.01 and not related if the kinship value is less than −0.025.

Samples with <5% missing SNP data were excluded. Next, principal component analysis was performed on the unrelated subset, then values were predicted along components of variations for all individuals in the related subset based on genetic similarities with individuals in the unrelated subset. The resulting components represented a model that can be used in place of static population frequencies to identify matches in a set of unknown individuals. The model used in this study was built using the GEDMatch database.

PC-Relate uses the principal components from PC-AiR and separates genetic correlations into two components: one for the sharing of alleles that are identical by descent from recent common ancestors and another for the sharing of alleles due to more distant common ancestors.

The components from PC-AiR were used to estimate allele frequencies based on the individual's ancestral background using linear regression instead of static population frequencies, such as those from gnoMAD. For two individuals, i and j, a kinship coefficient, ${\hat{ϕ}}_{i j}$ , was then calculated using the estimated allele frequencies, $û$ , from the PC-AiR model as follows:

where s is an SNP in S SNPs that were typed in both individuals, $g_{i s}$ and $g_{j s}$ are the number of reference alleles in i and j at SNP s, respectively, and $û_{i s}$ and $û_{j s}$ are the expected allele frequencies calculated by PC-AiR for i and j at SNP s, respectively. This algorithm is termed “whole genome kinship,” as it considers the entire genome as one segment of relatedness, signified by a whole genome kinship coefficient.

This whole genome kinship coefficient was used to identify relationships when referring to the whole genome kinship algorithm. When the number of SNPs shared between two individuals was <6000, a whole genome kinship coefficient of >0.031 between two individuals was required to be considered relatives in this study.

The threshold of 0.031 was determined by querying both the mock PM samples and the degraded Verogen family samples (Supplementary Table S2) against the GEDmatch database using the whole genome kinship coefficient algorithm and identified a threshold that met our expected specificity. To do this, we first filtered out kits that belonged to known family members.

We assumed that the remaining kits were unrelated and expected that a query of these mock PM samples against the database would return hits we knew were considered false matches. This assumption allowed us to determine the number of false positive matches based on a given whole genome kinship coefficient. We then filtered for hits with an SNP overlap of ≤6000, then iterated kinship coefficient thresholds starting with 0.00 and increased by 0.01 for each iteration, and calculated specificity for each threshold.

Because this application to DVI would result in a much smaller database as compared with GEDmatch, we can expect fewer false positive matches when querying a small, local database. This is because GEDmatch has greater than 1 million kits, and a private database will likely never reach that size, resulting in a lowered rate of false positives when querying such a small database.

Based on that assumption, we accepted a threshold of 0.031, which had a specificity of 99%. Setting a threshold that resulted in 100% specificity would result in zero false positives but might be too restrictive and classify higher order relationships as false negatives.

Snedecor et al. introduced an additional step, “windowed kinship,” to identify distant relationships more accurately. Windowed kinship consisted of calculating windows of kinship across the genome to find shared kinship segments. This was performed by enumerating all possible windows within each chromosome and calculating a kinship coefficient for all windows.

These windows were then filtered by a minimum kinship coefficient threshold and included in the shared cMs calculation. The filtered segments were then iterated and stretches of SNPs sharing at least one allele and two alleles were categorized separately. Shared cM of each window was then calculated across all segments by subtracting the cM of the last SNP by the cM of the first SNP in each window. Total shared cM was then calculated by summing the shared cM across all windows.

Thresholds for windowed kinship were determined previously.⁴¹ Briefly, total shared cM and the longest segment of cM were used to identify relationships when referring to the windowed kinship algorithm. When the number of SNPs shared between two individuals is between 6000 and 8000, the shared cM value must be above 180 and the longest segment of cM must be above 30 to be considered a relationship.

When the number of SNPs shared between two individuals is between 8000 and 9000, the shared cM value must be above 150 and the longest segment of cM must be above 30 to be considered a relationship. When the number of SNPs shared between two individuals is 9000 or more, the shared cM value must be above 140 and the longest segment of cM must be above 30 to be considered a relationship.

The whole genome kinship coefficient can be used to filter at any number of SNPs shared. However, Snedecor et al. observed a higher specificity when filtering on shared cM and longest segment cM (e.g., using windowed kinship) when the SNP overlap was greater than 6000, particularly for higher degrees of relationships.⁴¹

However, in cases where there were fewer than 6000 SNPs, overall sensitivity and specificity was lower using windowed kinship on close relationships. Windowed kinship is necessary for distant relationships; however, if there are not enough SNPs typed, it is difficult to have enough information to reliably call windowed segments at third-, fourth-, and fifth-degree relationships (Supplementary Fig. S1).⁴¹

More simply, the number of SNPs typed between two individuals (SNP overlap) decides when to use the whole genome kinship algorithm (<6000 SNPs overlap) and when to use the windowed kinship algorithm (greater than 6000 SNPs overlap). And, when one algorithm is decided on based on that SNP overlap, a value or a set of values are used to filter the data to identify relationships, depending on which algorithm was chosen.

The cutoffs for both whole genome kinship and windowed kinship were chosen to ensure a high sensitivity but more importantly, a high specificity as demonstrated in Snedecor et al. Supplementary Table S4.⁴¹ Lowering these thresholds may capture more relationships (i.e., increase sensitivity) but will introduce more false positive hits, particularly for more distant relationships (e.g., fourth- and fifth-degree) and when the database of samples is large (tens to hundreds of thousands of samples).

Results

Downsampled high coverage data—high plexity demonstration

To show feasibility, we first simulated high plexity with a set of 30 libraries generated from high-quality DNA samples (Supplementary Table S1) that were sequenced on the NextSeq 500 sequencer. The range of total reads per sample produced by the NextSeq was 8,086,090–32,707,490 with an average of 23,186,251 reads.

To simulate high plexity, we randomly selected reads from FASTQ files for each sample generated by the NextSeq sequencing run until a desired number of reads was met, which is termed downsampling. We simulated sequencing plexities of 16 and 30 on the MiSeq FGx with the MiSeq FGx Reagent Kit by downsampling the data to 1.5 M reads and 800,000 reads for each sample, respectively.

To determine the SNP call rate across samples, the number of typed SNPs was determined for each sample at each simulated sequencing plexity. The distribution of typed SNPs for the two simulated plexities for the 30 libraries is presented in Figure 1. The average recovery rate was 8586 (7781–8848) and 7234 (6375–7516) SNPs for the simulated 16 and 30 plexity sequencing runs, respectively.

Fig. 1.

The distribution of the number of SNPs typed in a downsampled dataset of Kintelligence libraries sequenced on a NextSeq 500 instrument. Downsampling is defined as the random selection of reads from a FASTQ file until the desired number of reads is achieved. Sixteen-sample run (16plex) data was generated by downsampling the raw reads in the FASTQ files to 1.5 million reads per sample, whereas the 30-sample run (30plex) data was generated by downsampling the raw reads to 800,000 reads per sample. Subsequent genotyping was performed, and the number of SNPs typed per sample was summed. For the 16-sample run, the minimum was 7781, the first quartile was 8472, the median was 8630, the third quartile was 8708, and the maximum was 8848. For the 30plex, the minimum was 6375, the first quartile was 7179, the median was 7299, the third quartile was 7382, and the maximum was 7516. SNPs, single-nucleotide polymorphisms.

Relationship determination with the Kintelligence kinship algorithm and LRs relies significantly on the number of typed SNPs in each sample in a 1-to-1 comparison, also called the SNP overlap. The higher the number of SNPs typed in both samples, the more true relationships will be identified (higher sensitivity). The average SNPs in common across all combinations of samples in both simulated runs were 6998 SNPs (6058–7322).

The simulations demonstrated that sequencing Kintelligence libraries resulted in fewer total reads for each sample, resulting in a smaller, stable set of SNPs. The SNP overlap was <8000, which may not be sufficient to identify higher order relationships (fourth- and fifth-degree)⁴¹ but should be enough for identifying third-degree relationships.

Impact of sequencing Kintelligence libraries at high plexities on SNP locus detection and heterozygosity

Simulated data described in the previous section demonstrated that increasing the throughput of Kintelligence libraries would decrease the number of SNPs typed per sample. Bearing this in mind, the number of SNPs typed, amount of heterozygosity lost, allele concordance, and number of SNPs in common between AM and PM samples sequenced at high throughput were evaluated.

Total coverage and SNP locus detection on mock AM and PM samples

Kintelligence libraries generated from mock AM samples (Supplementary Tables S1 and S3) were sequenced at the recommended plexity of 3 samples per sequencing run and again at four higher plexities of 12, 16, 24, and 32 samples per run. Total coverage was highest for the samples sequenced at three samples per run, as expected (Fig. 2A).

Fig. 2.

Total coverage and the number of SNPs typed in Kintelligence library samples sequenced on the MiSeq FGx in 3 (3plex), 12 (12plex), 16 (16plex), 24 (24plex), and 32 (32plex) sample runs for (A, C) mock AM samples (libraries generated with 1 ng intact DNA), and 3 (3plex) and 12 (12plex) sample runs for (B, D) mock PM samples (libraries generated with cremated, interred, and burnt bones, dental remains, degraded DNA samples from whole blood, and low input DNA samples [50, 100, 250, and 500 pg]). Total coverage is displayed in the range of Millions (M) of sequence reads. AM, antemortem; PM, postmortem.

All higher plexities demonstrated a large drop in total coverage; however, among the higher plexity runs, the coverage did not vary as much. The number of loci with reads below the AT increased as plexity increased for mock AM samples (Fig. 2C). In addition, the distribution of SNPs dropping below the AT widened as the sequencing plexity increased with an average of 10,111 (7653–9753).

The SNPs typed per sample were sequenced on a 32-sample run. Importantly, although coverage dropped by a large amount between 3plex and 32plex runs, the number of SNPs typed remained high. The minimum number of SNPs typed remained above 7000 at the highest plexity of 32, indicating that sequencing Kintelligence libraries at high plexity showed higher numbers of typed SNPs compared with the simulated results (Fig. 1).

The MFI samples are often degraded and can contain low levels of genomic DNA. Stochastic effects in PCR amplification lead to elevated numbers of loci falling below the AT, resulting in lower coverage and numbers of typed SNPs. Increasing the number of samples sequenced per run compounds these issues. To evaluate how increasing sequencing plexity compared with the standard Kintelligence protocol with 3 samples per run for degraded and low input DNA, libraries from 31 mock PM samples were generated (Supplementary Table S2).

These libraries were sequenced at the standard (3plex) and higher (12plex) plexities. Total coverage decreased when 12 samples per run were sequenced, as compared with the total coverage when 3 samples per run were sequenced (Fig. 2B). The number of SNPs typed above AT per sample was not vastly different between the two runs, although some samples sequenced at higher plexity had much fewer SNPs typed (Fig. 2D; 3plex average of 9313 SNPs, 7215–9991; 12plex average of 8796 SNPs, 4603–9903).

Similar to the AM samples, although coverage decreased when 12 samples per run were sequenced, the number of SNPs typed remained high and fairly comparable to the number of SNPs typed when 3 samples per run were sequenced. These results indicate that increasing the number of samples per run for degraded and/or low DNA input samples will increase the number of loci below AT, but the biggest impact on numbers of typed loci is the degradation state and/or quantity of DNA.

Allele concordance and heterozygosity of mock AM and PM samples

Increased plexity can also lead to elevated allele dropout and subsequent decreased heterozygosity, particularly for PM samples where the DNA is degraded and/or of low quantity. We evaluated the impact of higher plexity sequencing on loss of sister alleles for heterozygous loci by comparing SNP genotypes for samples sequenced at the recommended plexity of 3 libraries to the genotypes for the same libraries sequenced at a plexity of 32.

For each sample, genotypes were considered concordant if both runs returned the same alleles, and discordant otherwise. Allele concordance for each sample was calculated by dividing the number of concordant loci by the total number of loci typed in both sequencing runs. For mock AM samples, allele concordance decreased by an average of 1.9% between libraries sequenced at a plexity of 3 compared with 32 (0.50–2.8%; Fig. 3A, left y-axis).

Fig. 3.

Allele concordance and heterozygosity of (A) mock AM and (B) mock PM samples sequenced at 3 samples per run (3plex) and 32 samples per run (32plex), and 3plex and 12 samples per run (12plex) for mock AM and mock PM, respectively. Percent concordance was calculated by comparing the 32plex or 12plex with the respective 3plex run and classifying the same allele calls at each locus as concordant, and different allele calls as discordant (below AT) (A. Mock AM samples shown in blue bars. B. Mock PM samples shown in green bars). Both alleles must match to be considered concordant for autosomal loci. A percentage was calculated by dividing the number of matching loci by the total number of loci that both runs called for each sample. The range of loci called was (A) 7808 to 9667 with an average of 8610; and (B) 4580 to 9795 with an average of 8627. Percent heterozygosity was determined by summing the number of heterozygous loci and dividing that value by the total number of called loci (A) Mock AM samples: 3plex shown as orange line, 32plex shown as gray line. (B) Mock PM samples: 3plex shown as blue line, 1plex shown as yellow line. Mock PM samples in (B) were DNA extracted from dental remains (Tooth), blood of varying degradation levels (sample names starting with 1231 and 3551), a buried bone (PS2), a cremated bone (SHSU 1), embalmed bones (SHSU 2 and SHSU 3), burnt bones (SHSU 5, SHSU 6, SHSU 7), and a series of low DNA input of Coriell sample NA24385 (50, 100, and 250 pg). AT, analytical threshold; SHSU, Sam Houston State University.

Heterozygosity was determined by summing the number of heterozygous loci per sample and dividing that value by the total number of loci called. Samples sequenced at a plexity of 32 demonstrated the highest difference in heterozygosity compared with the standard plexity of 3 (6.8% average difference; Fig. 3A, right y-axis).

Allele concordance for mock PM samples decreased on average by 1.2% (Fig. 3B, left axis). Samples in Figure 3B were sorted based on decreasing degradation index, where applicable (Supplementary Table S2). Dental remains and cremated and burnt bones demonstrated the lowest allele discordance, with 0.6% average discordance.

The PS2 sample demonstrated the smallest allele concordance (n = 70.0%). Difference in heterozygosity between 3plex and 12plex of the sample aligned with the allele discordance (Fig. 3B, right axis). The 3plex run of PS2 contained 2515 loci (26% of SNPs typed) with the Imbalance flag whereas the 12plex run had 537 loci (7% of SNPs typed) with the Imbalance flag, which explains the large difference in heterozygosity between the 3plex and 12plex.

SNP loci detected in common (SNP overlap) across mock AM and PM samples

To evaluate kinship, an AM sample is compared with a PM sample to test for relatedness, becoming a 1-to-1 comparison. The number of SNPs typed in both individuals (SNP overlap) significantly impacts the determination of relationships, particularly distant relationships (e.g., fourth- and fifth-degree).⁴¹

We, therefore, calculated the SNP overlap among all combinations of samples between the mock PM samples sequenced at a plexity of 12 and the mock AM samples sequenced at a plexity of 32. The average overlap between combinations of mock PM and mock AM samples was 8020 (4071–9727).

These data indicate that, with an average of around 8000, the Kintelligence libraries sequenced at higher plexities (12 for PM and 32 for AM samples per run) will be able to identify first-, second-, third-, most fourth-degree relationships, and some fifth-degree relationships. Pairs that exhibited <8000 SNPs overlapping will likely not have enough data to match fourth- or fifth-degree relationships (Supplementary Table S4 published in Snedecor et al.⁴¹).

Evaluation of the Kintelligence kinship algorithm on anonymized GEDMatch samples with reduced SNP call rates

Given the range of SNPs typed, we investigated the ability of the lowered number of loci to determine kinship among a set of known relatives. Data presented by Snedecor et al. demonstrated that to identify higher order relationships (e.g., fourth- and fifth-degree), windowed kinship was more accurate when the SNP overlap was at least 6000 SNPs with the best performance at 9000 overlapping SNPs.⁴¹ As shown in Figure 2C and D, mock AM and mock PM samples sequenced at high plexity were less likely to reach this call rate.

We investigated the ability of the whole genome kinship algorithm (considers the entire genome as a single window) to identify relationships when the SNP overlap was <9000. One thousand anonymized query profiles were chosen from the GEDMatch database at random and were searched across the database to select 2954 “related” target profiles for relationships spanning first- to fifth degree as well as unrelated samples with zero shared cMs. The profiles were randomly filtered to 8000, 6000, 4000, and 2000 loci typed for each query-target pair.

The whole genome kinship algorithm was applied to test one-to-one relationships for every query-target pair at all SNP numbers for a total of 446,182 comparisons. Sensitivity and one minus the specificity are plotted in Figure 4. Ideally, all data points should fall in the top left corner of each plot. This indicates a low false positive rate and a high true positive rate.

Fig. 4.

Sensitivity and 1-specificity of whole genome kinship coefficients for anonymized GEDMatch samples. The loci of 1000 query samples and 2954 target samples were randomly filtered to 80%, 60%, 40%, or 20% call rates, resulting in 8000, 6000, 4000, or 2000 loci typed per sample. Query and target samples were paired, and the whole genome kinship coefficient was calculated for each pair using the Kintelligence kinship algorithm. A whole genome kinship coefficient threshold of greater than 0.031 was implemented to classify a relationship between two individuals, and any pairs with a kinship coefficient of ≤0.031 were considered unrelated. Sensitivity and specificity were calculated by comparing the results of the whole genome kinship coefficient with the results generated by the GEDMatch one-to-many tool, with the one-to-many tool representing the truth data. Axes are the same for all graphs.

If data points remain on the left-most part of the plot, this indicates that there was a low false positive rate but a higher false negative rate, where true relationships were classified as unrelated. First- and second-degree relationships exhibited a high specificity and sensitivity (Fig. 4), even for profiles with 2000 typed SNPs (average sensitivity of 70.6% and average specificity of 100%).

Sensitivity began to decrease for third-degree relationships, but the whole genome kinship algorithm performed better for relationship detection⁴¹ than the windowed kinship algorithm for profiles with 4000–8000 SNPs typed (Fig. 4). Fourth- and fifth-degree matches demonstrated low sensitivity and specificity for all levels of SNPs typed (Fig. 4).

Importantly, the specificity for all degrees remained above 99.98%, with 6000 and 8000 SNP ranges remaining above 99.9975% (Fig. 4). These results indicated that the whole genome kinship algorithm has a very low false positive rate. Overall, these simulated results show that this modified Kintelligence method for higher throughput sequencing coupled with the Kintelligence whole genome algorithm should work well for DVI and can be implemented on a local server.

Evaluation of high plexity Kintelligence on a known family

We tested the ability of the Kintelligence kinship algorithm to identify true relationships for Kintelligence libraries generated from a known family set of DNA samples sequenced at a high plexity. To determine the highest plexity that can maintain accuracy in identifying true relationships while excluding known unrelated pairs, these libraries were sequenced at plexities of 12, 16, 24, or 32 samples on an MiSeq FGx.

We chose the Utah/CEPH 1463 family from the Coriell Institute for Medical Research (Supplementary Fig. S2A). Total coverage decreased as plexity increased (Fig. 5A). The number of SNPs typed in each sample per run was determined to ensure enough SNPs were typed to perform kinship and LR calculations (Fig. 5B). The SNP call rate decreased as the number of samples per run increased (Fig. 5B) and the results were comparable to our initial high-throughput testing (Fig. 2C), indicating the reproducibility of high plexity sequencing of Kintelligence libraries.

Fig. 5.

(A) Distribution of total coverage and (B) distribution of the number of SNPs typed across samples in the Utah/CEPH 1463 family when sequenced with the MiSeq FGx using the ForenSeq^® Kintelligence kit at four different samples per run: 12 samples per run (12plex), 16 samples per run (16plex), 24 samples per run (24plex), and 32 samples per run (32plex). Total coverage is displayed in the range of Millions (M) of sequence reads.

We next calculated kinship metrics, including whole genome kinship coefficients, shared cM, longest segment cM, and LRs for all combinations of members in the pedigree. However, most members in this family were related, with only five pairs of individuals being unrelated.

To increase the number of unrelated controls, we included 100 randomly selected 1000 Genomes Project samples (Supplementary Table S4). To simulate a PM versus AM comparison as it would be performed for DVI, the 12-sample Kintelligence run samples were treated as mock PM samples. Each PM sample profile was then paired to each of the other sample profiles from all the sequencing runs (12, 16, 24, and 32 samples per run) as well as the 100 unrelated profiles from the 1000 Genomes Project to calculate all kinship metrics for each pair to identify related pairs.

Although the minimum SNP overlap was greater than 8900 among samples that were sequenced at all plexities (the libraries were generated with 1 ng of high-quality DNA), which would require the use of windowed kinship, both windowed kinship and whole genome kinship algorithms performed equally well, resulting in 100% sensitivity and 100% specificity (Supplementary Table S5).

Whole genome kinship coefficients clearly discriminate each of the expected relationships, with no false positive relationships detected (Fig. 6A). Taken together, all calculated metrics of kinship were able to accurately identify all relationships from the Utah/CEPH 1463 family, even when comparing profiles generated on the run with 12 samples to profiles generated on the run with 32 samples.

Fig. 6.

The distribution of the (A) kinship coefficient and (B) logLR for all combinations of pairs taken from the Utah/CEPH 1463 family and 100 randomly selected samples from the 1000 Genomes Project, which represented unrelated controls. The Utah/CEPH 1463 family includes grandparents, parents, and siblings. Each unique graph represents a comparison of the 12-sample run versus itself (12_vs_12), the 16-sample (12_vs_16), the 24-sample (12_vs_24), and 32-sample (12_vs_32) runs to simulate a PM versus AM comparison and to determine the highest plexity that is able to maintain accurate identification of a close relationship for AM samples. Samples were paired, and the kinship coefficient and logLR were calculated for each pair. A pair was considered related if the kinship coefficient was greater than 0.031, unrelated otherwise, which is represented by the black vertical line in A. A pair was considered related if the logLR was greater than 0, unrelated otherwise, which is represented by the black vertical line in B. logLR, log base 10 likelihood ratio.

The LRs demonstrated a similar trend (Fig. 6B). These data demonstrate that at a maximum plexity of 32 samples per run for mock AM Kintelligence libraries, profiles contain enough SNP data to accurately identify relationships up to second degree.

Assessment of increased plexity and the Kintelligence kinship algorithm of degraded/low DNA input samples with the Verogen family

Increased throughput on high-quality samples from the Utah/CEPH 1463 family demonstrated 100% sensitivity and 100% specificity in identifying true relationships (Fig. 6). This family contains only up to second-degree relationships (grandparents), and the DNA samples used to generate the Kintelligence libraries were all high quality and quantity. In addition, the Impact of Sequencing Kintelligence Libraries at High Plexities on SNP Locus Detection and Heterozygosity section shows that a plexity of 12 returned a high number of SNPs typed in degraded and/or low input samples; however, the samples included were not related and therefore could not be used to test the performance of the kinship algorithm at this higher plexity (as compared with the current standard 3plex).

The Verogen family was included to provide higher order relationships (up to fifth degree), some of which were degraded, to model DVI and missing persons casework samples. Sample replicates taken from the Verogen family (Supplementary Fig. S2B), which included all degrees of relationship up to fifth degree, were prepared with the Kintelligence kit and sequenced at a plexity of 30 samples.

Replicates from two members, V004 and V016, were artificially degraded, resulting in degradation indices of 2.1, 2.6, 5.1, and 20 for four replicates of V004 and 1.5, 2.0, 2.2, and 2.9 for four replicates of V016 (Supplementary Tables S2 and S3). Low DNA input libraries of V004 and V016 (50 pg) were also included in the 12-sample mock PM sequencing run.

Total coverage for the mock AM/intact samples was much lower (Fig. 7A) than the previous mock AM samples (Figs. 2A and 5A). This trend was likely due to the increased plexity and several samples within the degraded/low input series having low degradation indices and being sequenced at a 12plex instead of a 30plex (Supplementary Table S2).

Fig. 7.

Distribution of (A) total coverage and (B) the number of SNPs typed for the Verogen family that included first, second-, third-, fourth-, and fifth-degree relationships. Degraded and low input samples were sequenced at 12 samples per run (12plex), and intact samples were sequenced at 30 samples per run (30plex). Total coverage is displayed in the range of Millions (M) or thousands (K) of sequence reads.

As demonstrated in Figure 5A, mock AM samples sequenced at 12plex demonstrated high coverage when the samples were not degraded. If PM samples have low degradation indices and high DNA input amounts, we expect the total coverage to be higher as compared with highly degraded and/or low input PM samples.

These data confirm that degradation appears to affect total coverage more than increased plexity, as shown by lower average and smaller minimum of total coverage in the degraded/low input series as compared with the total coverage of the intact series in Figure 7A.

The mock AM samples show a higher number of typed SNPs with a tight distribution compared with the mock PM samples, where the number of typed SNPs was significantly lower and distributed across a wide range (Fig. 7B). Further, the mock AM sample libraries generated from intact DNA had fewer typed SNPs compared with the mock AM sample generated from commercially available DNA samples sequenced at a similar plexity (Fig. 7B compared with Figs. 2C and 5B).

The average percent difference in heterozygosity between the intact mock AM V004 samples sequenced at a plexity of 30 and corresponding degraded/low DNA input mock PM V004 samples sequenced at a plexity of 12 was 14.6%. The average loss of heterozygosity was higher than observed in Figure 3B, likely due to the combination of both degradation and low DNA input that ranged from 250 to 700 pg for the degraded V004 samples (Supplementary Table S2).

These results demonstrate that increased degradation levels, lowered DNA input, and, to a lesser extent, increased plexity caused a decrease in the number of SNPs typed in mock PM samples (Figs. 2B and 7A). We next sought to determine whether this decrease in SNPs typed was deterministic (specific loci dropout consistently across samples) or stochastic (random loci dropout across samples).

If locus dropout were deterministic, we would expect to see a core set of SNPs typed in common across several samples, as certain loci are more susceptible to dropout than others, resulting in a high SNP overlap between samples. Conversely, if locus dropout were stochastic, we would expect to see a small number of SNPs typed in common across several samples, as random locus dropout would mean that typed loci would not necessarily be the same between samples.

This would result in a significantly lower SNP overlap between samples. To determine whether locus dropout was stochastic or deterministic, we calculated the SNP overlap across eight artificially degraded samples from the Verogen family. The number of SNPs typed per sample, the SNP overlap between eight pairs of samples, and the SNP overlap across the eight artificially degraded samples are shown in Figure 8.

Fig. 8.

Venn diagram of the SNPs typed in eight low input or degraded samples from the Verogen family as indicated at the edges of each oval. Also included are the number of SNPs typed in common (SNP overlap) between the samples (values within the overlap of two ovals) and an overall SNP overlap across all samples, as indicated in the overlap of all ovals.

The SNP overlap between each pair of samples remained high and close to the lowest value between the pairs. The overall SNP overlap, 3007, represented 72% of the SNPs typed in the sample with the lowest number of SNPs typed, 4173. These results indicate that locus dropout was deterministic and not stochastic, where a stable set of SNPs typed can be expected across samples.

One-to-one comparisons were performed by pairing each degraded/low input mock PM sample with each intact mock AM sample to calculate whole genome kinship coefficients, shared cM, longest segment cM, and LRs (Supplementary Table S6). All relationships up to and including third degree were identified with 100% sensitivity and 100% specificity when using both the whole genome and windowed Kintelligence kinship algorithms (Supplementary Table S6). We were able to identify 8 of the 30 possible fourth-degree relationships with the Kintelligence windowed kinship algorithm but only two with the Kintelligence whole genome algorithm (Supplementary Table S6).

No fifth-degree relationships passed the kinship thresholds with either of the algorithms. These results confirm that first- to third-degree relationships can be determined for mock PM and mock AM sample Kintelligence libraries sequenced at high plexity with high specificity.

Discussion

The ForenSeq Kintelligence Library Prep Kit can be used to solve violent crimes and missing persons cases using forensic genetic genealogy though uploads to databases such as GEDMatch PRO³⁸ or FamilyTreeDNA.⁶⁷ An unidentified remains case in Oregon state was resolved in this way,⁴² where DNA was extracted from a molar tooth, amplified using Kintelligence, and the results uploaded to GEDmatch PRO.

The GEDMatch PRO algorithm was specifically designed to work with the 10,230 SNPs in the Kintelligence kit, but it requires upload to a public database to search for relatives, and a minimum of 6000 SNPs typed in the sample.³⁸ The current configuration of the Kintelligence kit allows for a maximum of three samples to be sequenced at a time on the MiSeq FGx, ensuring enough SNPs are typed for GEDMatch PRO upload.

If a relative is known but does not exist in one of the two databases, relatives of the missing person would be required to upload their profiles to the databases, which creates privacy concerns. In order for Kintelligence to be adopted for DVI, a higher throughput approach is required, and where data do not need to be uploaded to an online database.

This work focused on sequencing more Kintelligence libraries simultaneously, with the goal of maximizing plexity for both AM and PM samples, while maintaining accuracy in identifying relationships up to third degree. Based on simulations, we predicted that the optimal plexity of sequencing for this degree relationship determinations was 12 mock PM or 32 mock AM libraries per run (Fig. 1), which was confirmed by performing sequencing runs of mock AM and PM samples on the MiSeq FGx (Figs. 2, 5, and 8).

Allele dropout was slightly more elevated at 32plex for AM samples and 12plex for PM samples, as compared with the corresponding 3plex runs (Fig. 3). However, the average allele concordance remained above 95% for AM samples sequenced at 32plex and 94% for PM samples sequenced at 12plex. Differences in heterozygosity between 3plex and 32plex or 12plex for AM or PM samples, respectively, correlated with allele concordance, with PM samples demonstrating a larger difference in heterozygosity between 3plex and 12plex (Fig. 3).

Although more alleles fell below threshold for libraries sequenced at higher plexities, enough SNPs were typed in common between the mock AM and PM samples to allow up to third-degree relationships even for highly degraded and low DNA input samples (Fig. 4 and Supplementary Table S6).

We have shown that the overlap of SNPs called between PM and AM samples will decrease as the quality of the DNA decreases. This overlap significantly affects the performance of the kinship algorithms, reducing sensitivity at the higher degrees of relationship (Fig. 4). In the Verogen family analysis, which included both degraded and low input samples, we were able to identify all relationships up to and including third-degree relationships (Supplementary Table S6).

The average number of SNPs typed was too low for windowed kinship to be appropriate (<6000) for fourth- and fifth-degree relationships, meaning this kinship algorithm is unsuitable for this application. The whole genome kinship algorithm is not as sensitive as the windowed kinship algorithm at identifying fourth- and fifth-degree relationships but can identify third-degree relatives when the number of overlapping SNPs typed is <6000 (Supplementary Table S6).⁴¹

If PM samples are expected to be fourth- or fifth-degree relatives in the context of an MFI for example, we recommend maximizing the SNP overlap by decreasing the number of samples sequenced concurrently, especially if the input is low and/or the sample is highly degraded.

Alternatively, if PM samples are expected to be first-, second-, or third-degree relatives, a higher number of samples sequenced concurrently (up to 12) even for degraded and/or low input samples should not affect the ability of the kinship algorithm to identify these relationships. For AM samples, the number of SNPs typed must be maximized to ensure the highest level of SNP overlap. Up to 32 samples sequenced concurrently enabled all third-degree relationships to be identified in this work.

Conclusions

Modifying the configuration of the ForenSeq Kintelligence kit allowed the sequencing of libraries at plexities of up to 12 PM samples and 32 AM samples. The Kintelligence whole genome kinship algorithm identified, with excellent sensitivity and specificity, all relationships up to third degree, with a false positive rate of zero thanks to highly tuned thresholds.

This algorithm can provide privacy for victims and their families by being installed on a private server, with an additional method for LR calculations, enabling offline searching of relatives in the context of MFIs and correct determination of relationships. By coupling higher plexity sequencing and an offline kinship matching algorithm, we propose this novel ForenSeq Kintelligence method as a solution for DVI.

Footnotes

Acknowledgments

The authors would like to thank Yonmee Han and Richelle Barta for expert technical assistance in generating some of the ForenSeq Kintelligence libraries described in this study, many thanks to Rachel Houston (Sam Houston State University) and Mitchell M. Holland and Sydney Gaston Sanchez (The Pennsylvania State University) for providing bone samples used in this study.

Authors' Contributions

Conceptualization, K.M.S., G.P.; methodology, S.M.R., J.A., J.S., and K.M.S.; laboratory processing, J.A., C.Z.; software, S.M.R.; formal analysis, S.M.R., J.S., J.A., C.Z., and K.M.S.; data curation, J.A.; writing—original draft preparation, S.M.R., J.A., L.D., and K.M.S.; writing—review and editing, S.M.R., J.A., J.S., C.Z., L.D., G.P., and K.M.S.; visualization, S.M.R.; supervision, K.M.S.; project administration, K.M.S., G.P. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Buccal swabs were collected from volunteers who each signed an informed consent form authorizing the use of de-identified samples for research use and publication.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available within this article, tables, figures, and Supplementary Tables and figures.

Author Disclosure Statement

All authors are employed by Verogen. Verogen is actively developing this application into a commercial product.

Funding Information

No funding was received for this article.

Supplementary Material

References

Lau

, Tan

. After the Indian Ocean tsunami: Singapore's contribution to the international disaster victim identification effort in Thailand. Ann Acad Med Singap, 2005; 34(5):341–351.

Beauthier

J-P

, Valck

, Lefevre

, et al. Mass disaster victim identification: The tsunami experience. Open Forensic Sci J, 2009; 2(1):54–62.

Lessig

, Rothschild

. International standards in cases of mass Disaster Victim Identification (DVI). Forensic Sci Med Pathol, 2012; 8(2):197–199.

Brough

, Morgan

, Rutty

. The basics of disaster victim identification. J Forensic Radiol Imaging, 2015; 3(1):29–37.

de Boer

, Roberts

, Delabarde

, et al. Disaster victim identification operations with fragmented, burnt, or commingled remains: Experience-based recommendations. Forensic Sci Res, 2020; 5(3):191–201.

INTERPOL. INTERPOL Disaster Victim Identification Guide: Part “A”. 2018. Available from: https://www.interpol.int/en/How-we-work/Forensics/Disaster-Victim-Identification-DVI [Last accessed: December 1, 2022].

Holland

, Cave

, Holland

, et al. Development of a quality, high throughput DNA analysis procedure for skeletal samples to assist with the identification of victims from the World Trade Center attacks. Croat Med J, 2003; 44(3):264–272.

Budowle

, Bieber

, Eisenberg

. Forensic aspects of mass disasters: Strategic considerations for DNA-based human identification. Leg Med (Tokyo), 2005; 7(4):230–243.

Wright

, Mundorff

, Chaseling

, et al. A new disaster victim identification management strategy targeting “near identification-threshold” cases: Experiences from the Boxing Day tsunami. Forensic Sci Int, 2015; 250:91–97.

10.

Akhteruzzaman

, Hasan

, Hossain

, et al. Disaster victim identification by DNA analysis: The tazreen fashions garment fire incident experience in Bangladesh. Biores Commun. 2015; 1(2):116–120.

11.

Bertoglio

, Grignani

, Di Simone

, et al. Disaster victim identification by kinship analysis: The Lampedusa October 3rd, 2013 shipwreck. Forensic Sci Int Genet, 2020; 44:102156.

12.

Gettings

, Kiesler

, Faith

, et al. Sequence variation of 22 autosomal STR loci detected by next generation sequencing. Forensic Sci Int Genet, 2016; 21:15–21.

13.

Wendt

, King

, Novroski

NMM

, et al. Flanking region variation of ForenSeq™ DNA Signature Prep Kit STR and SNP loci in Yavapai Native Americans. Forensic Sci Int Genet, 2017; 28:146–154.

14.

Novroski

, King

, Churchill

, et al. Characterization of genetic sequence variation of 58 STR loci in four major major population groups. Forensic Sci Int Genet, 2016; 25:214–226.

15.

Devesse

, Ballard

, Davenport

, et al. Concordance of the ForenSeq™ system and characterisation of sequence-specific autosomal STR alleles across two major population groups. Forensic Sci Int Genet, 2018; 34:57–61.

16.

Phillips

, Devesse

, Ballard

, et al. Global patterns of STR sequence variation: Sequencing the CEPH human genome diversity panel for 58 forensic STRs using the Illumina ForenSeq DNA Signature Prep Kit. Electrophoresis, 2018; 39(21):2708–2724.

17.

Delest

, Godfrin

, Chantrel

, et al. Sequenced-based French population data from 169 unrelated individuals with Verogen's ForenSeq DNA signature prep kit. Forensic Sci Int Genet, 2020; 47:102304.

18.

Butler

, Shen

, McCord

. The development of reduced size STR amplicons as tools for analysis of degraded DNA. J Forensic Sci, 2003; 48(5):1054–1064.

19.

Hill

, Kline

, Coble

, et al. Characterization of 26 miniSTR loci for improved analysis of degraded DNA samples. J Forensic Sci, 2008; 53(1):73–80.

20.

Ambers

, Churchill

, King

, et al. More comprehensive forensic genetic marker analyses for accurate human remains identification using massively parallel DNA sequencing. BMC Genomics, 2016; 17(Suppl 9):750.

21.

Fattorini

, Previderé

, Carboni

, et al. Performance of the ForenSeq(™) DNA Signature Prep kit on highly degraded samples. Electrophoresis, 2017; 38(8):1163–1174.

22.

Watherston

, McNevin

, Gahan

, et al. Current and emerging tools for the recovery of genetic information from post mortem samples: New directions for disaster victim identification. Forensic Sci Int Genet, 2018; 37:270–282.

23.

Watherston

, Watson

, Bruce

, et al. Efficient DNA profiling protocols for disaster victim identification. Forensic Sci. 2021; 1(3):148–170.

24.

Senst

, Caliebe

, Scheurer

, et al. Validation and beyond: Next generation sequencing of forensic casework samples including challenging tissue samples from altered human corpses using the MiSeq FGx system. J Forensic Sci, 2022; 67(4):1382–1398.

25.

Senst

, Caliebe

, Drum

, et al. Recommendations for the successful identification of altered human remains using standard and emerging technologies: Results of a systematic approach. Forensic Sci Int Genet, 2023; 62:1–14.

26.

Brenner

, Weir

. Issues and strategies in the DNA identification of World Trade Center victims. Theor Popul Biol, 2003; 63(3):173–178.

27.

Birus

, Marcikić

, Lauc

, et al. How high should paternity index be for reliable identification of war victims by DNA typing?. Croat Med J, 2003; 44(3):322–326.

28.

Alonso

, Martin

, Albarran

, et al. Challenges of DNA profiling in mass disaster investigations. Croat Med J, 2005; 46(4):540–548.

29.

Graham

. Disaster victim identification. Forensic Sci Med Pathol, 2006; 2(3):203–207.

30.

Tillmar

, Sturk-Andreaggi

, Daniels-Higginbotham

, et al. The FORCE panel: An all-in-one SNP marker set for confirming investigative genetic genealogy leads and for general forensic applications. Genes, 2021; 12(12):1968.

31.

Turner

, Nagraj

, Scholz

, et al. Evaluating the impact of dropout and genotyping error on SNP-based Kinship analysis with forensic samples. Front Genet, 2022; 13:882268.

32.

Pontes

, Sousa

, Medeiros

. SNPs and STRs in forensic medicine. A strategy for kinship evaluation. Arch Med Sadowej Kryminol, 2017; 67(3):226–240.

33.

Phillips

, García-Magariños

, Salas

, et al. SNPs as supplements in simple Kinship analysis or as core markers in distant pairwise relationship tests: When do SNPs add value or replace well-established and powerful STR tests?. Transfus Med Hemother, 2012; 39(3):202–210.

34.

Sanchez

, Phillips

, Børsting

, et al. A multiplex assay with 52 single nucleotide polymorphisms for human identification. Electrophoresis, 2006; 27(9):1713–1724.

35.

Westen

, Matai

, Laros

, et al. Tri-allelic SNP markers enable analysis of mixed and degraded DNA samples. Forensic Sci Int Genet, 2009; 3(4):233–241.

36.

Zhang

, Wang

, Cheng

, et al. Complex kinship analysis with a combination of STRs, SNPs, and indels. Forensic Sci Int Genet, 2022; 61:102749.

37.

Verogen. ForenSeq Kintelligence Kit Datasheet: Document # VD2020054. Rev. A. 2021. 2021. Available from: https://verogen.com/wp-content/uploads/2021/03/forenseq-kintelligence-datasheet-vd2020054-a.pdf [Last accessed: December 1, 2022].

38.

Verogen. SNP Typing in Universal Analysis Software and Kinship Estimation with GEDmatch PRO: Document # VD2020058 Rev. A. 2021. 2021 Available from: https://verogen.com/wp-content/uploads/2021/03/snp-typing-uas-kinship-estimation-gedmatch-pro-tech-note-vd2020058-a.pdf [Last accessed: December 1, 2022].

39.

Verogen. High-Quality Outcomes from Low-Quality Samples with ForenSeq Kintellingence Application Note: Document # VD2021002 Rev. B. 2021. 2021. Available from: https://verogen.com/wp-content/uploads/2021/07/high-quality-outcomes-low-quality-samples-app-note-vd2021002-b.pdf [Last accessed: December 3, 2022].

40.

Peck

, Koeppel

, Gorden

, et al. Internal validation of the ForenSeq Kintelligence Kit for application to Forensic Genetic Genealogy. bioRxiv. 2022:2022.10.28.514056.

41.

Snedecor

, Fennell

, Stadick

, et al. Fast and accurate kinship estimation using sparse SNPs in relatively large database searches. Forensic Sci Int Genet, 2022; 61:102769.

42.

DNA Labs International is the first accredited lab to use the ForenSeq Kintelligence System to aid in the identification of unidentified remains [press release]. PRWeb: Cision, February 23, 2022.

43.

Kling

, Tillmar

, Egeland

. Familias 3—Extensions and new functionality. Forensic Sci Int Genet, 2014; 13:121–127.

44.

Budimlija

, Prinz

, Zelson-Mundorff

, et al. World Trade Center human identification project: Experiences with individual body identification cases. Croat Med J, 2003; 44(3):259–263.

45.

Gjertson

, Brenner

, Baur

, et al. ISFG: Recommendations on biostatistics in paternity testing. Forensic Sci Int Genet, 2007; 1(3-4):223–231.

46.

Bradford

, Heal

, Anderson

, et al. Disaster victim investigation recommendations from two simulated mass disaster scenarios utilized for user acceptance testing CODIS 6.0. Forensic Sci Int Genet, 2011; 5(4):291–296.

47.

Slooten

. Validation of DNA-based identification software by computation of pedigree likelihood ratios. Forensic Sci Int Genet, 2011; 5(4):308–315.

48.

Issan

, Avlas

, Daniel

. Investigation of different ways in which the CODIS 7.0 may be used in mass disaster identification. J Forensic Sci, 2019; 64(3):852–856.

49.

Verogen. Universal Analysis Software Kintelligence Module V2 Reference Guide: Document # VD2022002 Rev. A. 2022. Available from: https://verogen.com/wp-content/uploads/2021/05/universal-analysis-software-v2-0-reference-guide-vd2019002-e.pdf [Last accessed: December 3, 2022].

50.

Carrasco

, Inostroza

, Didier

, et al. Optimizing DNA recovery and forensic typing of degraded blood and dental remains using a specialized extraction method, comprehensive qPCR sample characterization, and massively parallel sequencing. Int J Legal Med, 2020; 134(1):79–91.

51.

Holt

, Stephens

, Walichiewicz

, et al. Human mitochondrial control region and mtGenome: Design and forensic validation of NGS multiplexes, sequencing and analytical software. Genes, 2021; 12:599.

52.

Loreille

, Diegoli

, Irwin

, et al. High efficiency DNA extraction from bone by total demineralization. Forensic Sci Int Genet, 2007; 1(2):191–195.

53.

Jakubowska

, Maciejewska

, Pawłowski

. Comparison of three methods of DNA extraction from human bones with different degrees of degradation. Int J Legal Med, 2012; 126(1):173–178.

54.

InnoGenomics. InnoQuant HY Human and male DNA quantification and degration assessment kit using 7500 real-time PCR system user guide v1.5 2016. Available from: https://innogenomics.com/wp-content/uploads/files/InnoQuant_HY_Using_7500_Real_Time_PCR_System_User_Guide_v1_5.pdf [Last accessed: December 1, 2022].

55.

Verogen. ForenSeq Kintelligence Kit Reference Guide: Document # VD2030053 Rev. B. 2021. Available from: https://verogen.com/wp-content/uploads/2021/03/forenseq-kintelligence-reference-guide-vd2020053-b.pdf [Last accessed: December 3, 2022].

56.

Illumina

. Understanding unique dual indexes (UDI) and associated library prep kits 2021. Available from: https://support.illumina.com/bulletins/2018/08/understanding-unique-dual-indexes—udi—and-associated-library-p.html [Last accessed: December 3, 2022].

57.

Verogen. MiSeq FGx Sequencing System Reference Guide. Revision F. 2021. 2019. Available from: https://verogen.com/wp-content/uploads/2021/02/miseq-fgx-system-reference-guide-vd2018006-f.pdf [Last accessed: December 3, 2022].

58.

Illumina

. Illumina. An introduction to Next-Generation Sequencing Technology. 2017. 2017. Available from: https://www.illumina.com/content/dam/illumina-marketing/documents/products/illumina_sequencing_introduction.pdf [Last accessed: December 3, 2022].

59.

Illumina

. NextSeq System Denature and Dilute Libraries Guide: Document # 15048776 v09 2018. Available from: https://support.illumina.com/content/dam/illumina-support/documents/documentation/system_documentation/nextseq/nextseq-denature-dilute-libraries-guide-15048776-09.pdf [Last accessed: December 3, 2022].

60.

Illumina

. NextSeq 500 System Guide: Document # 15046563 v06 2019. Available from: https://support.illumina.com/content/dam/illumina-support/documents/documentation/system_documentation/nextseq/nextseq-500-system-guide-15046563-07.pdf [Last accessed: December 3, 2022].

61.

Jager

, Alvarez

, Davis

, et al. Developmental validation of the MiSeq FGx Forensic Genomics System for targeted next generation sequencing in Forensic DNA casework and database laboratories. Forensic Sci Int Genet, 2017; 28:52–70.

62.

Illumina

. Illumina. Illumina Sequencing Technology. 2010. 2010. Available from: https://www.illumina.com/documents/products/techspotlights/techspotlight_sequencing.pdf [Last accessed: December 3, 2022].

63.

Gotoh

. An improved algorithm for matching biological sequences. J Mol Biol, 1982; 162(3):705–708.

64.

DNA Painter. Available from: https://dnapainter.com/tools/sharedcmv4 [Last accessed: December 6, 2022].

65.

Galván-Femenía

, Barceló-Vidal

, Sumoy

, et al. A likelihood ratio approach for identifying three-quarter siblings in genetic databases. Heredity, 2021; 126(3):537–547.

66.

Conomos

, Miller

, Thornton

. Robust inference of population structure for ancestry prediction and correction of stratification in the presence of relatedness. Genet Epidemiol, 2015; 276–293.

67.

Verogen. Verogen and Gene by Gene Form Groundbreaking Partnership to Accelerate Adoption for Forensic Investigative Genetic Genealogy: Businesswire; 2022. Available from: https://www.businesswire.com/news/home/20220815005116/en/Verogen-and-Gene-by-Gene-Form-Groundbreaking-Partnership-to-Accelerate-Adoption-of-Forensic-Investigative-Genetic-Genealogy [Last accessed: December 10, 2022].

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

0.07 MB

0.05 MB

0.03 MB

0.02 MB

0.15 MB

0.05 MB

Evaluation of a High-Throughput Dense Single-Nucleotide Polymorphism PCR Multiplex Next Generation Sequencing Application for Human Remains Identification

Abstract

Introduction

Materials and Methods

Human DNA samples and mock PM samples

ForenSeq Kintelligence library preparation

Sample pooling and sequencing

Primary, secondary, and tertiary data analysis

Downsampling analysis

GEDMatch data simulations for whole genome kinship algorithm testing

Calculation of LRs and ForenSeq Kintelligence kinship values

Results

Downsampled high coverage data—high plexity demonstration

Impact of sequencing Kintelligence libraries at high plexities on SNP locus detection and heterozygosity

Total coverage and SNP locus detection on mock AM and PM samples

Allele concordance and heterozygosity of mock AM and PM samples

SNP loci detected in common (SNP overlap) across mock AM and PM samples

Evaluation of the Kintelligence kinship algorithm on anonymized GEDMatch samples with reduced SNP call rates

Evaluation of high plexity Kintelligence on a known family

Assessment of increased plexity and the Kintelligence kinship algorithm of degraded/low DNA input samples with the Verogen family

Discussion

Conclusions

Footnotes

Acknowledgments

Authors' Contributions

Institutional Review Board Statement

Informed Consent Statement

Data Availability Statement

Author Disclosure Statement

Funding Information

Supplementary Material

References

Supplementary Material