Whole genome characterization of sequence diversity of 15,220 Icelanders

Understanding of sequence diversity is the cornerstone of analysis of genetic disorders, population genetics, and evolutionary biology. Here, we present an update of our sequencing set to 15,220 Icelanders who we sequenced to an average genome-wide coverage of 34X. We identified 39,020,168 autosomal variants passing GATK filters: 31,079,378 SNPs and 7,940,790 indels. Calling de novo mutations (DNMs) is a formidable challenge given the high false positive rate in sequencing datasets relative to the mutation rate. Here we addressed this issue by using segregation of alleles in three-generation families. Using this transmission assay, we controlled the false positive rate and identified 108,778 high quality DNMs. Furthermore, we used our extended family structure and read pair tracing of DNMs to a panel of phased SNPs, to determine the parent of origin of 42,961 DNMs.


Background & Summary
Characterization of genetic diversity is of paramount importance to population genetics. Advances in sequencing technologies have been instrumental in the creation of whole-genome sequencing data sets on the population level [1][2][3][4][5] . Sequencing data augmented with chip-typed samples have proven a valuable approach in genome-wide association studies, by combining the low cost of chip typing with the comprehensive characterization of sequence diversity via whole-genome sequencing [1][2][3]5 . In contrast, the analysis of DNMs and other recent mutations requires the sequencing of families in which these recent mutations have occurred. DNM detection is hampered by the fact that they are typically only observed in very few individuals (most often a single individual) and sequencing unrelated individuals is of limited value. Furthermore, analysis of singletons and DNMs is complicated since current sequencing technologies and calling algorithms identify a non-trivial number of false-positives 6 . The analysis of rare variants requires more extensive quality control than the analysis of common variants.
Families with well understood inheritance patterns provide useful tools for this purpose, such as monozygotic twin concordance and the transmission of DNM alleles to offspring of probands in threegeneration families. These methods provide genome-wide validation assays with no ascertainment beyond the family relationships of the probands. A useful feature of the Icelandic population are comprehensive genealogical records going back to 1,600 (and in some instances back to 740). The deCODE genetics genealogical databases currently includes records of 819,410 Icelanders 7 .
Here, we have sequenced 12,803 Icelanders in addition to the previously described 2,417 (ref. 7) to an average 34X coverage. This sequencing effort resulted in 39,020,168 autosomal variants passing GATK filters and for the entire variant set we were able to impute 37,127,995 variants into a set of 151,677 chip genotyped Icelanders.
We particularly focused on obtaining quality DNM calls. We calibrated our DNM calling with segregation of alleles in three-generation families, and subsequently validated our DNM calls with concordance among monozygotic twins. This resulted in data set of 108,778 high quality DNM calls, of which we were able to determine the parent of origin for 42,961 by observing their segregation in threegeneration families and read tracing DNMS with phased and imputed variants.

Methods
The following methods are an extension of those reported by Jónsson et al. 8 Chip genotyping and long range phasing We performed chip genotyping on 47,555 additional individuals relative to the set presented by Gudbjartsson et al. 7 . We used the OmniExpress24 chip for the majority of the samples (33,710; Table 1). This update resulted in a set of 151,677 individuals with chip genotypes. Finally, we long range phased (LRP) the data set using the algorithm described by Kong et al. 9 , resulting in a panel of 151,677 individuals with LRP chip genotypes.
Preparation of samples for whole-genome sequencing DNA was derived from two sample types: buccal swab and whole blood. Three different sample preparation kits from Illumina were employed: TruSeq DNA (Method A), TruSeq Nano (Method B) and TruSeq PCR-Free (Method C). Samples were prepared for sequencing according to the manufacturers' instructions (Illumina). In short, either 50 ng (Method B) or 1 μg (Methods A and C) of genomic DNA, isolated from either frozen blood samples or buccal swabs, was fragmented to a mean target size of 300-400 bp using a Covaris E220 instrument. End repair, generating blunt ended fragments was performed followed by size selection using different ratios of AMPure XP magnetic purification beads. 3′-Adenylation and ligation of indexed sequencing adaptors containing a T nucleotide overhang was performed, followed either by AMPure purification (Method C) alone or purification followed by PCR enrichment (10 cycles) using appropriate primers (Methods A and B). The quality and concentration of all sequencing libraries was assessed using either the Agilent 2,100 Bioanalyzer (12-samples) or the LabChip GX (96-samples) instrument from Perkin Elmer. Sequencing libraries were diluted and stored at −20°C. Further quality control of sequencing libraries was done by multiplexing and pooling either 24 or 96 samples and sequencing each pool on an Illumina MiSeq instrument to assess optimal cluster densities, library insert size, duplication rates and library diversities. All steps in the workflow were monitored using an in-house laboratory information management system (LIMS) with barcode tracking of all samples and reagents.

DNA whole genome sequencing
Sequencing libraries were hybridized to the surface of paired-end (PE) flowcells using the Illumina cBot. Paired-end sequencing-by-synthesis (SBS) was performed on Illumina sequencers, including GAII x , HiSeq 2,000/2,500 or HiSeq X instruments, respectively. Read lengths depended on the instrument and/or sequencing kit being employed and varied from 2 × 76 cycles to 2 × 150 cycles of incorporation and imaging. Real-time analysis involved conversion of image data to base-calling in real-time. The largest number of samples (approximately 12,000) were sequenced on HiSeq X instruments with read lengths of 2 × 150 cycles, using either v1 or v2 flowcells and sequencing chemistries, respectively ( Table 2).

Reference
The reference sequences used to map reads was the human genome assembly GRCh38, not including alternate assemblies (GCA_000001405.15_GRCh38_no_alt_analysis_set.fna) in addition to sequences determined to represent common contaminants in our sequencing pipeline. These sequences are the sequences of the bacteriophage PhiX (Enterobacteria_phage_phiX174_sensu_lato_uid14015/NC_001422. fna), the bacteria Ralstonia Pickettii (Ralstonia_pickettii_12D_uid58859) and two sequences from the human microbiome (Coprobacillus_D7_uid32495/NZ_EQ999972.fna, Coprobacillus_D7_uid32495/ NZ_EQ999922.fna).

Alignment of sequences
We aligned the raw sequences against the reference described in the previous section with BWA version 0.7.10 mem 10 . The sequences in the BAM files were realigned around indels with GenomeAnalysisTK-Lite/2.3.9 (ref. 11) using a public set of known indels and a set of indels previously discovered in the Icelandic data. We marked PCR duplicates with Picard tools 1.117.

Filtering of BAM files
We organized the sequencing data of each individual into one BAM file per lane on the flowcell. Much of the sequencing data generated for each individual prior to the HiSeqX machines derived from multiple lanes. In contrast, the sequencing data generated for each individual from HiSeqX machines largely derives from single lanes. In the following, we will refer to BAM files created before and after the merging per individual as BM-and AM-BAM files, respectively.
We performed a pileup using samtools BM-BAM files at all sites where the sample is homozygous according to their chip genotypes. We combined these pilups for each sample and calculated a mismatch rate as the number of bases not matching the chip genotype divided by the number of chip typed bases. We exlucded BM-BAM files with mismatch rates above 2% from the analysis.  We merged BM-BAM files into single AM-BAM files based on unique individual and sample type source combinations. The total number of AM-BAM files, prior to the subsequent quality filtering was 17,326. For contamination estimation, we used number and fraction of read pairs between pairs of heterzygote SNP markers supporting three or more haplotypes as a surrogate for the contamination level. We omitted AM-BAM files from subsequent analysis where more than 0.1% of SNP pairs had two or more reads supporting the third haplotype.
We removed AM-BAM files from the analysis if any of the underlying BM-BAM files failed any of the following criteria:  Table 3. The counts of machine and flowcell combinations for the AM-BAM files are listed in Table 4.

Variant calling
We applied a filtering step after merging the alignments and before the variant calling. We required the alignment to contain at least 45 matching bases (not necessarily consecutive) and removed parts of reads extending from template and into the adapter. For the adapter removal, we hard clip bases from the forward read extending further than the rightmost alignment of the reverse read and vice versa for the bases in the reverse read extending further than the leftmost alignment in the forward read.
Subsequently, we merge BAM files across all individuals in segments of 50,000 bases using SAMtools (version 1.3) 12 . These merged BAM files were used as input to GATK unified genotype caller 11 , resulting in a pooled genotype calling of 15,220 individuals. This resulted in 39,020,168 autosomal variants that passed the GATK recommended filters 1 , of which 31,079,378 and 7,940,790 were SNPs and indels, respectively (Table 5). Of the 15,220 individuals used for the sequence genotyping, we restricted the DNM identification to 14,688 individuals with over 20x coverage.

Variant phasing
To improve genotype quality and to phase genotypes, we used an iterative algorithm based on the IMPUTE HMM model 13 and LRP haplotypes 1,7 . This method is based on the principle that individuals that share haplotypes spanning the marker of interest are more likely to carry the same allele at the marker on the background of the shared haplotypes. For example, a pair of carriers of a rare sequence variant that share a haplotype are more likely to carry the rare variant on the shared haplotype than the two haplotypes that they do not share.

Variant imputation
Variants were imputed based on the IMPUTE HMM model 13 as previously described 1,7 , where chip genotyped individuals that share haplotypes with individuals in the set of sequenced individuals (the training set) are imputed to have the alleles on the background of the shared haplotypes.

Genotype imputation information
The informativeness of genotype imputation was estimated by the ratio of the variance of imputed expected allele counts and the variance of the actual allele counts:  where θ ∈ {0, 1} is the allele count. Var(E(θ/chip data)) was estimated by the observed variance of the imputed expected counts and Var(θ) was estimated by p(1 − p), where p is the allele frequency. The number of sequence variants with imputation information greater than 0.8 per variant type is listed in Table 5.

Extraction of de novo candidates
A schematic overview of the DNM identification procedure is shown in Fig. 1. Here we define allelic balance as the fraction of reads supporting the alternative allele out of the reads supporting the reference and alternative alleles combined. We used the genotype from GATK to define possible carriers. We define likely carriers of a DNM as those that meet the following requirements: read depth over 12, allele balance deviation from 0.5 less than 0.25 and genotype likelihood difference greater than 20 between the highest and second-highest scoring genotype (GQ). We restricted our DNM analysis to the primary assembly of the autosomes and chromosome X (hg38). We extracted DNM candidates from the variants satisfying the following criteria: • The proband must be an alternative allele carrier.
• For homozygous alternative allele carriers, we only consider candidates with 1 or no reads supporting the reference allele • Minimum depth in the parent of 12 reads for the autosomes and 6 reads for the X chromosome in males • No more than one read supporting the alternative allele in the parent • Maximum allelic balance for parent of 0.05 • Minimum depth of 12 reads for proband • Minimum allelic balance for proband of 0.15 • Maximum of 10 possible and 3 likely carriers beyond the descendants of the parent-pair • Maximum of 10% average soft clipping per read covering the DNM candidate   Furthermore, we removed probands from the DNM analysis with more than: 10% average soft clipping per read, 1.5% average fraction of Ns in the read alignment and 300 DNM candidates. This removed 15 probands, of which two are three-generation probands and one was a monozygotic twin proband.

QC-filtering using transmission in three-generation families
The DNM candidates identified in three-generation family probands were used to tune the DNM calling by tracing the transmission of alleles to offspring of the proband. We interpret discrepancy between haplotype sharing and segregation of a DNM as an indication that it is not present in the germline. More specifically, we dichotomized the DNMs according to whether they are consistent with haplotype transmission from the proband to his or her offspring. We restricted to DNMs where two or more offspring share distinct haplotypes from the proband in three-generation families.
We imposed a biological constraint on the DNM candidates on the X chromosome, i.e., male carriers must be homozygous for the alternative allele, otherwise we consider it inconsistent. To avoid inconsistency calls due to low quality genotypes in the offspring, we treated genotype calls of the offspring as missing if they do not meet the following requirements: at least 2 reads supporting an alternative allele, allelic balance over or equal to 0.1 and depth over or equal to 10 reads.
This evaluation of the DNM candidates gave us a binary outcome y, that we incorporate into the following Generalized Additive Model (GAM) Y∼sðPsroband À allelic À balanceÞ þ sðoxoG À FoxoG À metricÞ þ sðTrio À NPOSSÞ þ gatk À filter: s() denotes a smooth term and the covariates are described in the following listing: • Proband_allelic_balance, allelic balance of the proband.  Using the fitted model, we predicted response for all DNM candidates. We used response strictly greater than 0.8 to define high-quality DNMs (Fig. 2).

Phasing of de novo mutations
We used two complementary approaches to phasing DNMs, i.e., using haplotype sharing in threegeneration families and read pair tracing DNMs with phased variants. In the former approach, we determined the parent of origin as in our previous analysis 15 . For example, if an offspring of the proband is a carrier of the DNM and carries the paternal chromosome of the proband, we assigned the mutation to the father. Conversely, if the offspring is not a DNM allele carrier then we would assign the DNM to the mother. If there were conflicting phase assignments among multiple siblings, we treated the DNM as unphased. We restricted the haplotype sharing analysis to segments of genetic length of at least 0.8 cM and with 200 consecutive markers present in the LRP panel.
In the latter approach, we used the phased and imputed panel of markers as described above to serve as a reference panel. We counted the number of read pairs traversing each biallelic combination between the DNM variant and the neighboring imputed and phased variants. We also recorded the number of read pairs supporting three or more haplotypes and reads not consistent with any of the biallelic combinations. We aggregated these numbers of read pairs per parent pair and DNM site. Using this summary, we phased DNMs by assigning the DNM allele to the parental chromosome with the read support. We only considered DNMs with read support for only one parent-of-origin. Finally, we aggregated the phase information from both methods into a consensus phase. Of the 4,566 DNMs that were phased by both methods, we excluded 53 that were phased inconsistently (Table 6).

Code availability
We used publicly available software (URLs listed below) in conjunction with the above described algorithms in the sequencing processing pipeline (sections Variant phasing and Variant imputation).

Data Records
We have deposited two data records to the European Variation Archive (EVA), the 39,020,168 sequence variants passing the GATK filters in 15,220 sequenced Icelanders (Data Citation 1) and the 108,778 high quality DNMs (Data Citation 2) annotated with phase information.

Technical Validation
As we have previously validated the accuracy of our imputation pipeline for rare and common variants 7 , we focused our validation efforts on the DNM set.

Family relationship confirmation
There were reported 3,056 mother-child and 2,473 father-child connections among our 14,688 sequenced Icelanders (over 20X average coverage) as recorded in the genealogical database, of those we confirmed 2,981 and 2,386 based on haplotype sharing, respectively. This resulted in 1,548 trios and 225 threegeneration families after restricting the analysis to probands passing the quality control.

Three-generation set
We calibrated the cutoffs of our DNM calls via consistent segregation of alleles in three-generation families. We used consistent transmission as the dependent variable in a generalized additive model (GAM). We used the fitted GAM model to derive prediction for all the DNM candidates. A histogram of the prediction for all DNM candidates is shown in Fig. 2. There is a clean separation between variants that we predicted to have been transmitted to the next generation and those that are inconsistent with three-generation haplotype sharing.

Concordance among monozygotic twins
Only 8 of the 91 monozygotic twin probands were used for the DNM calling calibration. Thus, DNM concordance among the monozygotic twins serves as an independent validation metric of the DNM calls. We calculated the DNM discordance between each monozygotic twin proband and its monozygotic twin. We treated all genotype calls of the monozygotic twin as missing if the depth was less than 10 and heterozygous calls as missing if they did not meet the following requirements: at least 2 reads supporting an alternative allele and allelic balance over or equal to 0.1. We could verify the absence or presence of the genotype in the twin of the proband for 6,000 out of 6,034 DNMs and found low discordance (2.9%) for the high quality DNMs calls (Fig. 3). This shows that our DNM data set is of high quality genome-wide as there is no ascertainment of DNMs selected for validation beyond the monozygotic proband requirement.

Phase comparison
We used two complementary approaches to assign parent of origin to the DNMs, the haplotype sharing between the parents and offspring of three-generation proband, and physically linking DNM candidates with phased variants via read pair tracing. We compared the concordance and discrepancy between the phasing methods. We found a low discrepancy 1.16% between the phasing approaches ( Table 6).