Beta-Binomial Model for the Detection of Rare Mutations in Pooled Next-Generation Sequencing Experiments

Abstract

Against diminishing costs, next-generation sequencing (NGS) still remains expensive for studies with a large number of individuals. As cost saving, sequencing genome of pools containing multiple samples might be used. Currently, there are many software available for the detection of single-nucleotide polymorphisms (SNPs). Sensitivity and specificity depend on the model used and data analyzed, indicating that all software have space for improvement. We use beta-binomial model to detect rare mutations in untagged pooled NGS experiments. We propose a multireference framework for pooled data with ability being specific up to two patients affected by neuromuscular disorders (NMD). We assessed the results comparing with The Genome Analysis Toolkit (GATK), CRISP, SNVer, and FreeBayes. Our results show that the multireference approach applying beta-binomial model is accurate in predicting rare mutations at 0.01 fraction. Finally, we explored the concordance of mutations between the model and software, checking their involvement in any NMD-related gene. We detected seven novel SNPs, for which the functional analysis produced enriched terms related to locomotion and musculature.

1. Introduction

With high-throughput technologies, namely next-generation sequencing (NGS), scientists can analyze whole genome at once, being able to identify diseases caused by multiple genetic disorders. For the detection of rare mutations, targeted exome sequencing is used. A mutation represents any change in DNA sequence, such as substitution of single-nucleotide polymorphism (SNP) or small insertions or deletions of a few DNA bases (INDELs). In NGS analysis, amplified DNA sequences are broken down in subsequences and sequenced in parallel. Therefore, the NGS technology has become a valuable tool and strategy for massively parallel analysis of a large number of genes with successful identification of several Mendelian disease genes (Ng et al., 2010). However a single NGS experiment generates billions of nucleotides of DNA sequences making detection of rare mutations computationally very challenging.

The object in this article is neuromuscular disorders (NMD), which refer to a group of many different syndromes and diseases that either directly or indirectly impair the function of the skeletal muscles hampering locomotion. The disease might affect the muscle itself or alter its nerve where the signal is sent. Alterations often occur in the muscle cell membrane called sarcolemma as a result of one or more genetic mutations in those particular cells.

There are several ways to diagnose NMD and among them genetic testing might be the most promising, not only to ascertain the syndrome but also to predict a person's chance of developing or passing it to his or her offspring. Targeted gene and chromosomal tests are often performed in case of clinical and histological evidence, although about 40% of cases rely on changes in the responsible genes that cannot be picked up. Kaplan releases annually a table listing the different neuromuscular diseases and the mutated genes. The latest version (Kaplan and Hamroun, 2014) comprises 16 different syndromes encompassing 685 disease phenotypes involving 360 different genes, 92 of them still unidentified. The first NGS experiment in NMDs has been done on a patient affected by Charcot–Marie–Tooth neuropathy (Lupski et al., 2010). Since then, several studies have created a validated NGS approach to the diagnosis of NMDs using the sequencing of panel genes (Valencia et al., 2012, 2013; Vasli et al., 2012; Böhm et al., 2013). Historically, Sanger sequencing and deletion duplication analysis were among the most exploited ones, however, these techniques are very low throughput and inefficient with regard to the large heterogeneity of the molecular parts involved.

Notwithstanding its cost–effectiveness, the NGS approach is still very expensive for studies that involve a large number of patients. In this case, an alternative is grouping individuals in pools and sequencing their genome together. Pooling approach was applied for the first time in 2010 (Futschik and Schlötterer, 2010), providing a reliable and cost-effective result for genome-wide allele frequency estimates. Since then, further NGS pooling experiments (Calvo et al., 2010; Wang et al., 2011; Ding et al., 2012; Ferraro et al., 2014) were used for rare mutation detection.

Many software have been developed to detect SNPs and other kinds of polymorphisms (hereafter variant calling software). They all work in the context of NGS data analysis aiming to find nucleotides that differ from their relative reference genome of interest. Various mathematical approaches from frequents to Bayesian were implemented in variant calling software. CRISP (version 071812) (Bansal, 2010) has been exclusively developed for pooling data in NGS experiments. It has been especially created to identify rare variants, first comparing the distribution of allele counts across multiple pools using contingency table, instead of running pool separately like other software do, and then calculating the probability of observing nonreference base calls due to sequencing errors. The authors compared CRISP with SNPseeker (Druley et al., 2009), VarScan (Koboldt et al., 2009), and MAQ (Li et al., 2008) detecting variants in data with minor allele frequency (MAF) at 2% with sensitivity 86.3% and specificity 97% (Cushing et al., 2013).

In SNVer (0.5.3 release) (Wei et al., 2011), binomial models are implemented for the modeling of sequencing error rates and variant distribution simultaneously. Common (>5%) and rare (0.1%–1%) variants are investigated in one integrated model. The Simes method is used to calculate pooled p value for a candidate being a variant from multiple pools. SNVer does not discard loci with low depth of coverage to achieve reliable variant calls. SNVer performance was compared with CRISP indicating binomial models' advantages versus frequentist Fisher's exact test approach (Wei et al., 2011).

The Genome Analysis Toolkit (GATK, version 3.4-46) (DePristo et al., 2011) uses a Bayesian genotype likelihood model to estimate simultaneously the most likely genotype and allele frequency of pooled samples, emitting a genotype of each sample (UnifiedGenotyper walker). This is handled by keeping information of low depth of coverage, which might be discarded in other variant calling software. GATK has already been used for pooled data showing high specificity and sensitivity (Calvo et al., 2012; Pihlstrøm et al., 2014; Savarese et al., 2014), however, its accuracy suffered in finding rare (<2%) variants in various studies (Bonachea, 2013; Gui et al., 2014; Spencer et al., 2014).

FreeBayes (version v0.9.10-3-g47a713e-dirty) (Garrison and Marth, 2012) uses a Bayesian haplotype-based method. It calls variants based on the literal sequences of reads aligned to a particular target, not their precise alignment. This method avoids one of the core problems with alignment-based variant detection, where software has to operate on a single position at a time, the locus, given that identical sequences may have multiple possible alignments. In a very recent rare mutation study (Clevenger et al., 2015), FreeBayes and GATK share very little concordance in SNP detection (5.3% of FreeBayes SNPs and 23.7% of GATK SNPs), and the first one failed in calling a large number of rare mutations, being more prone in finding false-positive SNPs.

Beyond variant calling software, Flaherty et al. (2011) used a beta-binomial model for the detection of rare single-nucleotide mutations at the 0.1% level for synthetic DNA samples. Also, authors analyzed nine clinical samples of H1 N1 influenza A and detected an oseltamivir (antiviral therapy) resistance mutation in the H1 N1 neuraminidase gene at a sample fraction of 0.18%. For the calculation of probability that specific position is mutated, the multireference framework was set up using diluted samples as reference. Proposed method became as background for the development of new variant calling software RVD (Cushing et al., 2013).

Recently, Ferraro et al. (2014) presented a detailed method to detect rare and putative causative mutations in pooled experiments on 98 genes in a population of 128 patients affected by a clinical diagnosis of muscular disease. They proposed a computational procedure consisting of two steps: first, decisional rules are applied to filter data, and then, a supervised classification technique is performed and compared for specificity and sensitivity with four variant calling software.

Apparently, although there are variant calling software available and promising methods to detect rare mutation developed, the results, meaning method/model sensitivity and specificity, are very diverse and typically less than 1. Therefore, it is valuable to test existing models to detect rare mutations for different types of data to better understand their limitations and trade-offs. Research of concordance between software/models is very scarce; however, it could very favorable for selection of techniques in future pooled NGS experiments.

The objective of this article is the detection of rare mutations in pooled NGS experiments, using the multireference beta-binomial model. We apply the model to find mutated positions in 98 genes of 128 patients, with a clinical diagnosis of muscular disease. We bring novelty using the statistical approach for detecting rare mutations (a nonreference rate equal to 0.01) in human exome data. We propose the multireference framework for pooled data with ability being specific up to two patients. We expect that the implemented multireference approach will be suitable for detection of rare mutation comparing results with the variant calling software discussed earlier. The aim to compare the results with a study by Ferraro et al. (2014) determined our selection of software. Finally, we analyze how many mutations are identified by all five methods and check whether they might be involved in any NMD-related gene.

2. Methods

2.1. DNA preparation and pooling

One hundred twenty-eight DNA samples were collected and exome sequenced through Illumina technology (HaloPlex Target Enrichment System For Illumina Sequencing, Protocol version D, August 2012; Agilent Technologies, Santa Clara, CA) from patients affected by NMD.

All patients were grouped in eight pools of 16 patients each (consisting one control sample), called main pools. Other eight pools were created replicating two samples of each pool, as presented in Tables 1 and 2, which we call reference pools. For the comparison reasons, we use the same pool organization as in Ferraro et al. (2014). The Human Genome version 19 (also called GRCh37) was used as reference genome.

Table 1.

Pooling Organization: Main Pools (1–8)

Pool
P1	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16
P2	17	18	19	20	21	22	23	24	25	26	27	28	29	30	31	32
P3	33	34	35	36	37	38	39	40	41	42	43	44	45	46	47	48
P4	49	50	51	52	53	54	55	56	57	58	59	60	61	62	63	64
P5	65	66	67	68	69	70	71	72	73	74	75	76	77	78	79	80
P6	81	82	83	84	85	86	87	88	89	90	91	92	93	94	95	96
P7	97	98	99	100	101	102	103	104	105	106	107	108	109	110	111	112
P8	113	114	115	116	117	118	119	120	121	122	123	124	125	126	127	128

The control samples are indicated in bold.

Table 2.

Pooling Organization: Reference Pools (9–16)

Pool	P1	P1	P2	P2	P3	P3	P4	P4	P5	P5	P6	P6	P7	P7	P8	P8
P9	1	2	17	18	33	34	49	50	65	66	81	82	97	98	113	114
P10	3	4	19	20	35	36	51	52	67	68	83	84	99	100	115	116
P11	5	6	21	22	37	38	53	54	69	70	85	86	101	102	117	118
P12	7	8	23	24	39	40	55	56	71	72	87	88	103	104	119	120
P13	9	10	25	26	41	42	57	58	73	74	89	90	105	106	121	122
P14	11	12	27	28	43	44	59	60	75	76	91	92	107	108	123	124
P15	13	14	29	30	45	46	61	62	77	78	93	94	109	110	125	126
P16	15	16	31	32	47	48	63	64	79	80	95	96	111	112	127	128

The control samples are indicated in bold.

2.2. Beta-binomial model

Binomial mixtures arise naturally if one assumes that the probability of success θ at each trial is not constant but it varies according to mixing distribution g(θ). A common choice is the beta distribution as mixing distribution, giving rise to the beta-binomial distribution.

Read count data at each position are decomposed to reference and nonreference reads. Observed random variable number of nonreference reads r_ij is distributed according binomial distribution with parameter θ_ij, which is true error rate at position j in pool i: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} {r_{ij}} \vert { \theta _{ij}},{n_{ij}} \sim Bin \big ( { \theta _{ij}},{n_{ij}} \big ) , \tag{1} \end{align*} \end{document}

where n_ij is the total read count; i = 1, 2, …, N and j = 1, 2, …, J.

Unobserved random variable θ_ij is distributed as beta distribution as follows: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} { \theta _{ij}} \vert { \mu _j},M \sim Beta \big ({ \mu _j},M \big ) , \tag{2} \end{align*} \end{document}

and \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} {\it Pr} \Big ({r_{ij}} \vert { \theta _{ij}},{n_{ij}} \Big ) = \left( {_{{r_{ij}}}^{{n_{ij}}}} \right) \theta _{ij}^{rij}{ \Big ( 1 - { \theta _{ij}} \Big )^{{n_{ij}} - {r_{ij}}}} \tag{3} \end{align*} \end{document} \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \it Pr \Big ( { \theta _{ij}} \vert { \mu _{ij}},M \Big ) = Bin \Big ({ \mu _j}M , ( 1 - { \mu _j} ) M \Big ) \theta _{ij}^{{ \mu _j}M - 1}{ ( 1 - { \theta _{ij}} ) ^{ ( 1 - { \mu _j} ) M - 1}}. \tag{4} \end{align*} \end{document}

We observe that the variance of nonreference rate \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$ { \frac { { r_ { ij } } } { { n_ { ij } } } } $$ \end{document} is homogenous across positions.

Maximum likelihood function for r|θ is as follows: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} {L_r} ( r \vert \theta , \mu , M ) = \mathop \prod _{i = 1}^N \prod \limits_{j = 1}^{J ( i ) }P ( {r_{ij}} \vert { \theta _{ij}},{ \mu _j},M ). \tag{5} \end{align*} \end{document}

The random variables and parameters are estimated by the maximum likelihood method using the expectation maximization (EM) algorithm. The algorithm is initialized with the plug-in of \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\mu$$ \end{document} and M moment estimates. The θ estimate has closed form solution as follows: \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} \begin{align*} \widehat \theta = { \frac { \widehat { \mu _J } \widehat M + { r_ { ij } } - 1 } { \widehat M + { n_ { ij } } - 2 } } , \tag { 6 } \end{align*} \end{document}

while parameters \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $$\mu$$ \end{document} and M are estimated by a constrained interior-point algorithm. The EM algorithm alternates optimizing complete log-likelihood until the change in the log-likelihood is < \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${10^{ - 4}}$$ \end{document} , as in the study by Flaherty et al. (2011). Binomial distribution was approximated by a Gaussian distribution, and z-test applied to test whether estimated true error rate is significantly different from the observed nonreference rate.

The earlier described model was used in the multireference framework, as presented in Table 3 and described in Section 3.1. We estimate beta-binomial model parameters for reference pools only. Estimated true error rate from reference pools was compared with observed nonreference rate of a main pool. Differences assessed as significant if p value was lower than \documentclass{aastex}\usepackage{amsbsy}\usepackage{amsfonts}\usepackage{amssymb}\usepackage{bm}\usepackage{mathrsfs}\usepackage{pifont}\usepackage{stmaryrd}\usepackage{textcomp}\usepackage{portland, xspace}\usepackage{amsmath, amsxtra}\pagestyle{empty}\DeclareMathSizes{10}{9}{7}{6}\begin{document} $${10^{ - 3}}$$ \end{document} .

Table 3.

Pool Management Scheme for Beta-Binomial Model

Pool	Reference pool	Patient		Pool	Reference pool	Patient
1	10, 11, 12, 13, 14, 15, 16	1	2	5	10, 11, 12, 13, 14, 15, 16	65	66
	9, 11, 12, 13, 14, 15, 16	3	4		9, 11, 12, 13, 14, 15, 16	67	68
	9, 10, 12, 13, 14, 15, 16	5	6		9, 10, 12, 13, 14, 15, 16	69	70
	9, 10, 11, 13, 14, 15, 16	7	8		9, 10, 11, 13, 14, 15, 16	71	72
	9, 10, 11, 12, 14, 15, 16	9	10		9, 10, 11, 12, 14, 15, 16	73	74
	9, 10, 11, 12, 13, 15, 16	11	12		9, 10, 11, 12, 13, 15, 16	75	76
	9, 10, 11, 12, 13, 14, 16	13	14		9, 10, 11, 12, 13, 14, 16	77	78
	9, 10, 11, 12, 13, 14, 15	15	16		9, 10, 11, 12, 13, 14, 16	79	80
2	10, 11, 12, 13, 14, 15, 16	17	18	6	10, 11, 12, 13, 14, 15, 16	81	82
	9, 11, 12, 13, 14, 15, 16	19	20		9, 11, 12, 13, 14, 15, 16	83	84
	9, 10, 12, 13, 14, 15, 16	21	22		9, 10, 12, 13, 14, 15, 16	85	86
	9, 10, 11, 13, 14, 15, 16	23	24		9, 10, 11, 13, 14, 15, 16	87	88
	9, 10, 11, 12, 14, 15, 16	25	26		9, 10, 11, 12, 14, 15, 16	89	90
	9, 10, 11, 12, 13, 15, 16	27	28		9, 10, 11, 12, 13, 15, 16	91	92
	9, 10, 11, 12, 13, 14, 16	29	30		9, 10, 11, 12, 13, 14, 16	93	94
	9, 10, 11, 12, 13, 14, 15	31	32		9, 10, 11, 12, 13, 14, 16	95	96
3	10, 11, 12, 13, 14, 15, 16	33	34	7	10, 11, 12, 13, 14, 15, 16	97	98
	9, 11, 12, 13, 14, 15, 16	35	36		9, 11, 12, 13, 14, 15, 16	99	100
	9, 10, 12, 13, 14, 15, 16	37	38		9, 10, 12, 13, 14, 15, 16	101	102
	9, 10, 11, 13, 14, 15, 16	39	40		9, 10, 11, 13, 14, 15, 16	103	104
	9, 10, 11, 12, 14, 15, 16	41	42		9, 10, 11, 12, 14, 15, 16	105	106
	9, 10, 11, 12, 13, 15, 16	43	44		9, 10, 11, 12, 13, 15, 16	107	108
	9, 10, 11, 12, 13, 14, 16	45	46		9, 10, 11, 12, 13, 14, 16	109	110
	9, 10, 11, 12, 13, 14, 15	47	48		9, 10, 11, 12, 13, 14, 16	111	112
4	10, 11, 12, 13, 14, 15, 16	49	50	8	10, 11, 12, 13, 14, 15, 16	113	114
	9, 11, 12, 13, 14, 15, 16	51	52		9, 11, 12, 13, 14, 15, 16	115	116
	9, 10, 12, 13, 14, 15, 16	53	54		9, 10, 12, 13, 14, 15, 16	117	118
	9, 10, 11, 13, 14, 15, 16	55	56		9, 10, 11, 13, 14, 15, 16	119	120
	9, 10, 11, 12, 14, 15, 16	57	58		9, 10, 11, 12, 14, 15, 16	121	122
	9, 10, 11, 12, 13, 15, 16	59	60		9, 10, 11, 12, 13, 15, 16	123	124
	9, 10, 11, 12, 13, 14, 16	61	62		9, 10, 11, 12, 13, 14, 16	125	126
	9, 10, 11, 12, 13, 14, 15	63	64		9, 10, 11, 12, 13, 14, 16	127	128

2.3. Filtering and software parameters

The read data for the beta-binomial model were filtered out that the coverage of the loci is equal to at least 500. Also, fixing a nonreference or error rate equal to 0.01, a mutation is characterized by one base with a frequency in the interval (1.1%, 8.3%) and two bases below 1%. The latter rule was selected in line with Ferraro et al. (2014) for the comparison reason. INDELs (insertions/deletions) were excluded from the analysis.

All software were set up with same parameters of minimum base quality (10), minimum read mapping (20), and pool size (number of haploid genomes in each pool, 32). Filtering criteria were the same as described for the beta-binomial model. Other specific software parameters were as follows: for CRISP, the maximum number of noncorrespondence between sample and reference was set up at 6, since the length of the reads was at 75 base pairs, and the “discrete” value for the pool was at 1; for GATK, the standard minimum confidence threshold for calling and the standard minimum confidence threshold for emitting were at 20.0; for FreeBayes, beyond INDELs, we excluded complex events and pool was considered “discrete” with evaluation of only the best N SNPs (“-use-best-n-alleles”), with N = 2, this to save computational cost of the analysis given the pool size (Garrison, pers. comm.), and last, the minimum alternate fractions were set up at 0.01.

Being untagged pools, a mix of nontagged individual mutations, there is no possibility to catch up mutations specific to single patients, for our experiments (5), however, given the way our pools were arranged (the first 8 made of a sequential number of 16 individuals each from 1 to 128, the last 8 made of the first two individuals from 1 to 64, respectively, and the last 14 individuals sequentially from 17 till 128), we were able to catch up specificity for couple of individuals. For example, couple of patients 1 and 2 were specified by overlapping pools 1 and 9 software variant calling results and couple 3–4 overlapping pools 1 and 10 (Ferraro et al., 2014).

2.4. Benchmark data set

For the assessment of our results, the benchmark data set was used (Ferraro et al., 2014). As we exclude INDELs from analysis, the data set includes 189 positions biologically confirmed with Sanger sequencing. Of those, 133 were found to be true positives (mutations) and 56 were false positives.

2.5. Structural and functional annotation

We performed a structural and functional annotation of those SNPs found simultaneously by all four software and the beta-binomial model for the couple of individuals 1–2 (overlapping pools 1–9). A Venn diagram performed at http://bioinformatics.psb.ugent.be/webtools/Venn/gave the list of shared SNPs. The list was pruned out of those found in the benchmark data set. Left over SNPs were searched under the SNP repository database dbSNP www.ncbi.nlm.nih.gov/SNP/ and Variant Effect Predictore online software /www.ensembl.org/info/docs/tools/vep/index.html, where their MAF information was recorded. A variant call format file was created through the software Kaviar3 http://db.systemsbiology.net/kaviar/cgi-pub/Kaviar.pl.

The genetic structural annotation with effect prediction analysis was performed using the command-line software SnpEff (Cingolani et al., 2012). This software annotates SNPs (finding their genome localization) and calculates the effect they produce on known genes. Subsequently, the list of SNPs was run in SIFTdbSNP at http://sift.jcvi.org/www/SIFT_dbSNP.html, which predicts whether SNPs that have a putative deleterious effect in amino acid substitution affect their protein function (Kumar et al., 2009). The list was then filtered out of those synonymous variation SNPs (neutral effect) keeping just the missense nonsynonymous ones. A Gene Ontology (GO) attribution term enrichment analysis (functional annotation) was finally performed using FuncAssociate 2.1 (Berriz et al., 2009), which aims to find properties enriched, given a list of genes. It uses Fisher's exact test to identify GO attributes that functionally describe a fraction of the entries in the input list. Enriched GO was then characterized by the online software AmiGO 2 (Carbon et al., 2009).

3. Findings

3.1. Multireference framework

In the analysis, we used untagged pooled NGS data, where DNA samples of different individuals were combined in a single experiment, still providing information with an uncertainty limited up to two patients.

For the detection of possibly mutated positions per two patients, we run 64 beta-binomial models in total. Proposed multireference framework for the pools was organized as follows. For example, to obtain mutated positions for the patients 1 and 2, we use pool P1 as the main and pools P10 to P16 as the reference (Table 3). We discard pool P9 as a reference as it contains information from patients 1 and 2. Hence, we estimate beta-binomial model parameters from seven pools (for example, pools P10–P16) and test the fit to the main pool data. Consequently, each output of a single beta-binomial model provides information about mutated positions up to two patients.

Initially, we had around 400,000 read positions per model. After filtering, the number of possible mutated positions diminished to 5460–14109 depending on the model. The total number of positions J, that is, positions of the main pool that were detected in reference pools, was different in each model.

3.2. Sensitivity and specificity analysis

For the assessment of multireference framework using the beta-binomial model and software, we calculated performance characteristics using a benchmark data set (see Section 2.4). One hundred thirty-three positions were confirmed as mutated, while 56 were verified as false positive using Sanger sequencing. Accuracy, sensitivity, specificity, positive predicted value (PPV), and negative predicted value (NPV) were calculated for the model and software used. For the calculation of the latter characteristics, instead of using seven reference pools, we used four pools with minimum variance.

From Table 4, one can read that the most accurate results were obtained using CRISP software. Beta-binomial model and SNVer gave the most specific results, and CRISP and FreeBayes gave the most sensitive results. The PPV, except GATK, was almost 1 for all software and the model.

Table 4.

Performance Characteristics of Beta-Binomial Model and Software Used

Model/software	Accuracy, %	Sensitivity, %	Specificity, %	PPV, %	NPV, %
GATK	75.1	69.9	87.5	93.0	55.1
SNVer	69.0	57.1	100.0	100.0	47.2
FreeBayes	89.1	85.7	98.0	99.1	72.5
CRISP	94.2	92.5	98.2	99.2	84.6
Beta-binomial	70.8	62.7	100.0	100.0	42.9

GATK, Genome Analysis Toolkit; NPV, negative predicted value; PPV, positive predicted value.

3.3. Concordance analysis

Similarity of the results assessed doing pairwise comparison of identified mutated positions between the software and model. A number of mutually identified mutated positions presented as matrices in Table 5, where on the diagonal the total number of identified mutated positions for specific software and the model. The percentage of pairwise comparison is calculated as row percentage of diagonal.

Table 5.

Pairwise Concordance of Possibly Mutated Positions by the Software Analyzed and Beta-Binomial Model

Software/model	GATK	SNVer FreeBayes		CRISP	Beta-binomial
			Mutual positions, n
GATK	450477	4653	17129	17127	1562
SNVer	4653	8229	7950	7006	113
FreeBayes	17129	7950	41380	23906	273
CRISP	17127	7006	23906	38196	282
Beta-binomial	1562	113	273	282	19343
			Mutual positions, %
GATK		1.0	3.8	3.8	0.3
SNVer	56.5		96.6	85.1	1.4
FreeBayes	41.4	19.2		57.8	0.7
CRISP	44.8	18.3	62.6		0.7
Beta-binomial	8.1	0.6	1.4	1.5

From Table 5, GATK, FreeBayes, and CRISP detected the largest number of positions as mutated. The number of possibly mutated positions identified by beta-binomial model was twice smaller compared to GATK, FreeBayes, and CRISP. However, the smallest number of potentially mutated positions was detected using SNVer and almost all these positions were captured by FreeBayes and CRISP. The largest mutual agreement was observed between GATK and SNVer software. Possibly mutated positions detected by beta-binomial model and GATK were different as we obtained low concordance with other software.

Given that having main and reference pools we can be specific up to two patients, we analyzed similarity of obtained results for one specific couple in more detail. From Table 6 we see that most often there are up to three mutations detected per couple in the benchmark data set. Four couples have six mutations and, specifically, couple of patient 1 and 2 has six mutations. We verified whether all mutations where detected by the model and all software.

Table 6.

Frequency of Mutations Per Couple

Mutations per couple, n	0	1–3	4–6	>6	Total
Confirmed positions	12	42	4	0	133
Nonconfirmed positions	41	11	8	4	56

Therefore, we explored the software and model agreement for the couple of patient 1 and 2 separately. Apparently, all software and the model detected 10 mutated positions for the couple. Of which three positions were listed in the benchmark data set and seven positions were not listed in benchmark data set. Similarity of the results for all patients, all positions detected by SNVer were also detected by other software or beta-binomial model. We provide functional annotation for seven positions not listed in the benchmark data set in the following section.

3.4. Functional annotation

Ten SNPs (specific of couple of individuals 1–2) were concordantly detected from the four variant calling software used in this study and the beta-binomial model (Fig. 1). Three of them were filtered out because they were already present in the benchmark data set. Noticeably, among the seven left over positions (Table 7), four mapped on chromosome 2, while the other three on chromosome 5, 14, and 19. All the SNPs, but one, were found in dbSNP. Their global MAF (GMAF) is relatively low.

FIG. 1.

Overlapping identified mutations by the software and model for patient 1 and 2. BB, beta-binomial model.

Table 7.

Structural Annotation of Novel Polymorphisms Detected by Beta-Binomial Model

CHR	POS	ID	REF/ALT	Effect	Functional	Gene	GMAF	SIFT
chr2	71709120	rs12470028	G/A	INT.	NONE	DYSF	0.0929	NA
chr2	158594953	.	G/A	N.SYN.	MIS.	ACVR1	NA	0.06
chr2	179403946	rs182683829	C/T	N.SYN.	MIS.	TTN	NA	NA
chr2	179421609	rs9808036	C/T	SYN.	SIL.	TTN	0.0843	NA
chr5	155935720	rs11748588	T/C	INT.	NONE	SGCD	0.0028	NA
chr14	64498016	rs200842904	G/A	N.SYN.	MIS.	SYNE2	0.0018	NA
chr19	13482623	rs17846942	G/C	DOW.	NONE	AC005513.1	0.0561	NA

CHR, chromosome; POS, position; ID, dbSNP identification number; REF, reference allele; ALT, alternate allele; GMAF, global minor allele frequency; SIFT, sorting intolerant from tolerant.

Four SNPs were annotated on exons, three of them with a nonsynonymous missense effect and one with a synonymous silent effect. The last one, together with other three SNPs not found in coding protein regions (one in an intron, one downstream of a gene, and the last one of unknown position), was pruned out of the list. SIFTdbSNP was not able to find any nonsynonymous SNPs eventually affecting protein function. The three nonsynonymous SNPs were localized inside genes ACVR1, SYNE2, and TTN1. FuncAssociate 2.1 gave four total GO terms: “locomotion” shared by all the three genes and “actin filament-base movement,” “Z disc,” and “ameboidal cell migration” shared by combinations of two. All GO presented an adjusted p value below than 0.05.

4. Discussion

In this article, we detected rare mutations in pooled NGS experiments using the multireference beta-binomial model. The model was applied to find mutated positions in 98 genes of 128 patients with a clinical diagnosis of muscular disease. We demonstrated that the used statistical approach is relevant in detecting rare mutations in human exome data.

Comparing results with currently available variant calling software (such as GATK, SNVer, FreeBayes, and CRISP), we obtained that the used approach is highly specific outperforming GATK, FreeBayes, and CRISP. This is in line with results by Flaherty et al. (2011), where 99% specificity was demonstrated for synthetic DNA at 0.1% sample fraction. However, the sensitivity of the model is 62.7% for pooled human exome NGS data. The sensitivity of beta-binomial model can be substantially increased by increasing base quality to 30 (Cushing et al., 2013). Also, to raise the power of the method to detect rare mutations, read depth should be very high (Flaherty et al., 2011). Median of total reads in benchmark data set was equal to 3000, and therefore, the increase in filtering criteria would lead to filtering out large proportion of confirmed mutated positions from the analysis. Although sensitivity of beta-binomial model ought to be improved, the PPV is 100. This indicates that the multireference approach applying beta-binomial model is accurate in predicting rare mutations at 0.01 fraction. Comparing our findings with Ferraro et al. (2014), we see that the accuracy of beta-binomial model is comparable with many supervised classification techniques. We obtained better accuracy for almost all software, indicating that the software is improving over time.

So far, very few studies have shown low concordance among different variant calling software when using same input file (B. Chapman, unpublished results at www.bcbio.wordpress.com; O'Rawe et al., 2013; Yu et al., 2013; Clevenger et al., 2015; Li, 2015). We provide the concordance analysis of beta-binomial model and the software confirming low agreement level (Table 5). GATK and beta-binomial model showed the lower concordance percentages when compared with each other and to other software. SNVer, CRISP, and FreeBayes showed much comparable overlapping percentages when compared with each other. Having such good concordance among SNVer, FreeBayes, and CRISP, one might consider using one instead of all three software. Low concordance is also due to the fact that variant calling software behave much worse when dealing with rare mutation detection than common ones (O'Rawe et al., 2013; Gui et al., 2014; Clevenger et al., 2015).

Indeed, by overlapping variant calling results of all the software and beta-binomial model, at least for one couple of individuals, we were able to detect three rare mutation SNPs (among a panel of 10 SNPs), which were not found before by Sanger experiments. The first one we describe here, SNP position 158594953 of chromosome 2 (Table 7), represents a not known DBSNP found by the Exome Aggregation Consortium (ExAC, Cambridge, MA; http://exac.broadinstitute.org, accessed August 2015), with a very low population frequency (5.766 × 10⁻⁵), localized in the exon number 9 of the gene ACVR1, which encodes for the protein activin A receptor, type I, also known as activin receptor-like kinase-2. Activins are growth factor proteins that can stimulate the production of follicle-stimulating hormone. ACVR1 is found in the 2015 NMD Kaplan list (Kaplan and Hamroun, 2014) as one of those genes where a mutation can cause a dystrophy disorder. As a matter of fact, an amino acidic substitution (arginine206histidine, db-SNP: rs121912678) (Shore et al., 2006; Bocciardi et al., 2009; Kaplan et al., 2009) of its protein has shown to cause fibrodysplasia ossificans progressiva (FOP, OMIM 35100), which provokes a progressive ossification of skeletal muscle tissues.

The second nonsynonymous missense rare mutation is DBSNP rs182683829 (Table 7). GMAF information for this SNP is not known but Variant Effect Predictor was able to find an EAMAF (NHLBI-ESB European American population allele frequency), which showed to be very low (0.0004) confirming its rareness. It is localized in the exon 181 of gene TTN (present in Kaplan list) (Kaplan and Hamroun, 2014) that encodes protein Titin. Titin is responsible for the passive elasticity of muscles. There are different pathologies generated by mutation correlated to this gene, two of them NMD related: tibial muscular dystrophy (TMD, OMIM 600334) and limb-girdle muscular dystrophy type 2J (LGMD2J, OMIM 608807). The first one is known to be associated with the presence of a particular INDEL (11-BP DEL/INS [dbSNP: rs281864927]) (Hackman et al., 2002) or the presence of an SNP (such as leucine 34315 proline [db-SNP: rs267607156]) (Hackman et al., 2002) and is a condition that affects the muscles at the front of the lower leg.

Last significant nucleotide substitution found here is the DBSNP rs200842904. The GMAF is 0.0018. Its location is in exon 44 of gene SYNE2 (Table 7, included in Kaplan list) (Kaplan and Hamroun, 2014), which encodes for the spectrin repeat-containing nuclear envelope protein-2 or nesprin-2, and is responsible for the structural integrity of the organelle hosting DNA in the human cells, the nucleus. A mutation (threonine 89 methionine [dbSNP: rs36215895]) (Zhang et al., 2007) found in nesprin-2 is already known to be responsible for the Emery-Dreifuss muscular dystrophy 5 (EDMD5, OMIM 612999), which is a condition that chiefly affects muscles used for movement (skeletal muscles) and heart (cardiac) muscle.

Functional annotation of these three genes enriched the “locomotion” (GO:0040011) pathway. AmiGO 2 defines this GO as “Self-propelled movement of a cell or organism from one location to another.” Other enriched GO pathways came out from three different combinations of two genes: (1) “actin filament-based movement” (GO:0030048), defined as “Movement of organelles or other particles along actin filaments, or sliding of actin filaments past each other, mediated by motor proteins”; (2) “Z disc” (GO:0030018), defined as “Platelike region of a muscle sarcomere to which the plus ends of actin filaments are attached”; (3) “ameboidal cell migration” (GO:0001667), defined as “Cell migration that is accomplished by extension and retraction of a pseudopodium.” All of these GOs seem to be related to muscle functionality. Thus, according to our structural and functional analysis, we suggest to consider these three SNPs for further experimental checks to investigate whether they might be involved in any NMD-related disorders.

We conclude that the results of multireference beta-binomial model are very promising for real NGS pooled data. However, we should admit that the multireference framework is challenging. One should properly think how to set up pooled NGS experiments that reference samples/pools available for multireference analysis. Another improvement is on the beta-binomial side, where more heterogeneous and less quality data might be considered.

Footnotes

Acknowledgments

This work has been partially funded by the Italian Flagship Project Interomics and by the Project PON02 00619. For Audrone Jakaitiene, the contribution to this article has been performed at the High Performance Computing and Networking Institute, National Research Council, Via Pietro Castellino 111, 80131 Naples, Italy. The authors thank Mr. Trerotola Giuseppe, working at the High Performance Computing and Networking Institute, National Research Council, Via Pietro Castellino 111, 80131 Naples, Italy, for his technical support in realizing the experiments.

Author Disclosure Statement

No competing financial interests exist.

References

Bansal

2010. A statistical method for the detection of variants from next-generation resequencing of DNA pools. Bioinformatics, 26, i318–i324.

Berriz

, Beaver

, Cenik

, et al. 2009. Next generation software for functional trend analysis. Bioinformatics, 25, 3043–3044.

Bocciardi

, Bordo

, Di Duca

, et al. 2009. Mutational analysis of the ACVR1 gene in Italian patients affected with fibrodysplasia ossificans progressiva: confirmations and advancements. Eur. J. Hum. Genet., 17, 311–318.

Böhm

, Vasli

, Malfatti

, et al. 2013. An integrated diagnosis strategy for congenital myopathies. PLoS One, 8, e67527.

Bonachea

E.M.

2013. Use of a targeted next generation sequencing approach for the study of a cardiac valve malformation with complex polygenic heritability. Presented at the 2013 AAP National Conference and Exhibition. American Academy of Pediatrics, Orlando, FL, USA.

Calvo

, Compton

, Hershman

, et al. 2012. Molecular diagnosis of infantile mitochondrial disease with targeted next-generation sequencing. Sci. Transl. Med. 25, 118ra10.

Calvo

S.E.

, Tucker

E.J.

, Compton

A.G.

, et al. 2010. High-throughput, pooled sequencing identifies mutations in NUBPL and FOXRED1 in human complex I deficiency. Nat. Genet., 42, 851–858.

Carbon

, Ireland

, Mungall

, et al. 2009. AmiGO: online access to ontology and annotation data. Bioinformatics, 25, 288–289.

Cingolani

, Platts

, Coon

, et al. 2012. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of drosophila melanogaster strain w1118; iso-2; iso-3. Fly, 6, 80–92.

10.

Clevenger

, Chavarro

, Pearl

, et al. 2015. Single nucleotide polymorphism identification in polyploids: a review, example, and recommendations. Mol. Plant, 8, 831–846.

11.

Cushing

, Flaherty

, Hopmans

, et al. 2013. RVD: A command-line program for ultrasensitive rare single nucleotide variant detection using targeted next-generation DNA resequencing. BMC Res. Notes, 6, 206.

12.

DePristo

, Banks

, Poplin

, et al. 2011. A framework for variation discovery and genotyping using next-generation DNA sequencing data. Nat. Genet., 43, 491–498.

13.

Ding

, Bashashati

, Roth

, et al. 2012. Feature-based classifiers for somatic mutation detection in tumour-normal paired sequencing data. Bioinformatics, 28, 167–175.

14.

Druley

T.E.

, Vallania

F.L.

, Wegner

D.J.

, et al. 2009. Quantification of rare allelic variants from pooled genomic DNA. Nat. Methods, 6, 263.

15.

Ferraro

M.B.

, Savarese

, Di Fruscio

, et al. 2014. Prediction of rare single-nucleotide causative mutations for muscular diseases in pooled next-generation sequencing experiments. J. Comput. Biol., 21, 665–675.

16.

Flaherty

, Natsoulis

, Muralidharan

, et al. 2011. Ultrasensitive detection of rare mutations using next-generation targeted resequencing. Nucleic Acids Res. 40, e2.

17.

Futschik

, and Schlötterer

2010. The next generation of molecular markers from massively parallel sequencing of pooled DNA samples. Genetics, 186, 207–218.

18.

Garrison

, and Marth

2012. Haplotype-based variant detection from short-read sequencing. arXiv Preprint arXiv, 1207.3907.

19.

Gui

, Bao

, Tang

, et al. 2014. Targeted next-generation sequencing on hirschsprung disease: a pilot study exploits DNA pooling. Ann. Hum. Genet., 78, 381–387.

20.

Hackman

, Vihola

, Haravuori

, et al. 2002. Tibial muscular dystrophy is a titinopathy caused by mutations in TTN, the gene encoding the giant skeletal-muscle protein titin. Am. J. Hum. Genet., 71, 492–500.

21.

Kaplan

, Xu

, Seemann

, et al. 2009. Classic and atypical fibrodysplasia ossificans progressiva (FOP) phenotypes are caused by mutations in the bone morphogenetic protein (BMP) type I receptor ACVR1. Hum. Mutat., 30, 379–390.

22.

Kaplan

J.C.

, and Hamroun

2014. The 2015 version of the gene table of monogenic neuromuscular disorders (nuclear genome). Neuromuscul. Disord., 24, 1123–1153.

23.

Koboldt

D.C.

, Chen

, Wylie

, et al. 2009. VarScan: variant detection in massively parallel sequencing of individual and pooled samples. Bioinformatics, 25, 2283–2285.

24.

Kumar

, Henikoff

, and Ng

2009. Predicting the effects of coding non-synonymous variants on protein function using the SIFT algorithm. Nat. Protoc., 4, 1073–1081.

25.

2015. Toward better understanding of artifacts in variant calling from high-coverage samples. Bioinformatics, 30, 2843–2851.

26.

, Ruan

, and Durbin

2008. Mapping short DNA sequencing reads and calling variants using mapping quality scores. Genome Res. 18, 1851–1858.

27.

Lupski

J.R.

, Reid

J.G.

, Gonzaga-Jauregui

, et al. 2010. Whole-genome sequencing in a patient with Charcot-Marie-Tooth neuropathy. N. Engl. J. Med., 362, 1181–1191.

28.

S.B.

, Nickerson

D.A.

, Bamshad

M.J.

, et al. 2010. Massively parallel sequencing and rare disease. Hum. Mol. Genet., 19, R119–R124.

29.

O'Rawe

, Jiang

, Sun

, et al. 2013. Low concordance of multiple variant-calling pipelines: practical implications for exome and genome sequencing. Genome Med. 5, 1.

30.

Pihlstrøm

, Rengmark

, Bjørnarå

, et al. 2014. Effective variant detection by targeted deep sequencing of DNA pools: an example from Parkinson's disease. Ann. Hum. Genet., 78, 243–252.

31.

Savarese

, Di Fruscio

, Mutarelli

, et al. 2014. Motorplex provides accurate variant detection across large muscle genes both in single myopathic patients and in pools of DNA samples. Acta Neuropathol. Commun. 2, 100.

32.

Shore

, Xu

, Feldman

, et al. 2006. A recurrent mutation in the BMP type I receptor ACVR1 causes inherited and sporadic fibrodysplasia ossificans progressiva. Nat. Genet., 38, 525–527.

33.

Spencer

, Tyagi

, Vallania

, et al. 2014. Performance of common analysis methods for detecting low-frequency single nucleotide variants in targeted next-generation sequence data. J. Mol. Diagn., 16, 1525–1578.

34.

Valencia

C.A.

, Ankala

, Rhodenizer

, et al. 2013. Comprehensive mutation analysis for congenital muscular dystrophy: a clinical PCR-based enrichment and next-generation sequencing panel. PLoS One, 8, e53083.

35.

Valencia

C.A.

, Rhodenizer

, Bhide

, et al. 2012. Assessment of target enrichment platforms using massively parallel sequencing for the mutation detection for congenital muscular dystrophy. J. Mol. Diagn., 14, 233–246.

36.

Vasli

, Böhm

J.S.

, Le Gras

, et al. 2012. Next generation sequencing for molecular diagnosis of neuromuscular diseases. Acta Neuropathol. 124, 273–283.

37.

Wang

, Pradhan

, Ye

, et al. 2011. Estimating allele frequency from next-generation sequencing of pooled mitochondrial DNA samples. Front. Genet. 2, 51.

38.

Wei

, Wang

, Hu

, et al. 2011. SNVer: a statistical tool for variant calling in analysis of pooled or individual next-generation sequencing data. Nucleic Acids Res. 39, e132.

39.

, and Sun

2013. Comparing a few SNP calling algorithms using low-coverage sequencing data. BMC Bioinform. 14, 274.

40.

Zhang

, Bethmann

, Worth

, et al. 2007. Nesprin-1 and −2 are involved in the pathogenesis of Emery Dreifuss muscular dystrophy and are critical for nuclear envelope integrity. Hum. Mol. Genet. 16, 2816.

Pool
P1	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16
P2	17	18	19	20	21	22	23	24	25	26	27	28	29	30	31	32
P3	33	34	35	36	37	38	39	40	41	42	43	44	45	46	47	48
P4	49	50	51	52	53	54	55	56	57	58	59	60	61	62	63	64
P5	65	66	67	68	69	70	71	72	73	74	75	76	77	78	79	80
P6	81	82	83	84	85	86	87	88	89	90	91	92	93	94	95	96
P7	97	98	99	100	101	102	103	104	105	106	107	108	109	110	111	112
P8	113	114	115	116	117	118	119	120	121	122	123	124	125	126	127	128

Pool	P1	P1	P2	P2	P3	P3	P4	P4	P5	P5	P6	P6	P7	P7	P8	P8
P9	1	2	17	18	33	34	49	50	65	66	81	82	97	98	113	114
P10	3	4	19	20	35	36	51	52	67	68	83	84	99	100	115	116
P11	5	6	21	22	37	38	53	54	69	70	85	86	101	102	117	118
P12	7	8	23	24	39	40	55	56	71	72	87	88	103	104	119	120
P13	9	10	25	26	41	42	57	58	73	74	89	90	105	106	121	122
P14	11	12	27	28	43	44	59	60	75	76	91	92	107	108	123	124
P15	13	14	29	30	45	46	61	62	77	78	93	94	109	110	125	126
P16	15	16	31	32	47	48	63	64	79	80	95	96	111	112	127	128

Pool
P1	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16
P2	17	18	19	20	21	22	23	24	25	26	27	28	29	30	31	32
P3	33	34	35	36	37	38	39	40	41	42	43	44	45	46	47	48
P4	49	50	51	52	53	54	55	56	57	58	59	60	61	62	63	64
P5	65	66	67	68	69	70	71	72	73	74	75	76	77	78	79	80
P6	81	82	83	84	85	86	87	88	89	90	91	92	93	94	95	96
P7	97	98	99	100	101	102	103	104	105	106	107	108	109	110	111	112
P8	113	114	115	116	117	118	119	120	121	122	123	124	125	126	127	128

Pool	P1	P1	P2	P2	P3	P3	P4	P4	P5	P5	P6	P6	P7	P7	P8	P8
P9	1	2	17	18	33	34	49	50	65	66	81	82	97	98	113	114
P10	3	4	19	20	35	36	51	52	67	68	83	84	99	100	115	116
P11	5	6	21	22	37	38	53	54	69	70	85	86	101	102	117	118
P12	7	8	23	24	39	40	55	56	71	72	87	88	103	104	119	120
P13	9	10	25	26	41	42	57	58	73	74	89	90	105	106	121	122
P14	11	12	27	28	43	44	59	60	75	76	91	92	107	108	123	124
P15	13	14	29	30	45	46	61	62	77	78	93	94	109	110	125	126
P16	15	16	31	32	47	48	63	64	79	80	95	96	111	112	127	128

Pool
P1	1	2	3	4	5	6	7	8	9	10	11	12	13	14	15	16
P2	17	18	19	20	21	22	23	24	25	26	27	28	29	30	31	32
P3	33	34	35	36	37	38	39	40	41	42	43	44	45	46	47	48
P4	49	50	51	52	53	54	55	56	57	58	59	60	61	62	63	64
P5	65	66	67	68	69	70	71	72	73	74	75	76	77	78	79	80
P6	81	82	83	84	85	86	87	88	89	90	91	92	93	94	95	96
P7	97	98	99	100	101	102	103	104	105	106	107	108	109	110	111	112
P8	113	114	115	116	117	118	119	120	121	122	123	124	125	126	127	128

Pool	P1	P1	P2	P2	P3	P3	P4	P4	P5	P5	P6	P6	P7	P7	P8	P8
P9	1	2	17	18	33	34	49	50	65	66	81	82	97	98	113	114
P10	3	4	19	20	35	36	51	52	67	68	83	84	99	100	115	116
P11	5	6	21	22	37	38	53	54	69	70	85	86	101	102	117	118
P12	7	8	23	24	39	40	55	56	71	72	87	88	103	104	119	120
P13	9	10	25	26	41	42	57	58	73	74	89	90	105	106	121	122
P14	11	12	27	28	43	44	59	60	75	76	91	92	107	108	123	124
P15	13	14	29	30	45	46	61	62	77	78	93	94	109	110	125	126
P16	15	16	31	32	47	48	63	64	79	80	95	96	111	112	127	128