SNP and indel identification and characterization

In total, 10,871,465 distinct SNPs were identified in the 44 genomes, with an average of 3.3 million SNPs per genome (Table 1). The average SNP transition to transversion ratio was 2.13, and the average SNP heterozygote to homozygote ratio was 1.56 (Table S1), both consistent with previous reports [7], [9], [21]. In addition, we identified a total of 1,350,484 distinct short insertions (range 1–76 bp) and 1,464,731 distinct short deletions (range 1–192 bp) in our sample, with an average of 207,000 short insertions and 214,000 short deletions in an individual genome (Table 2, Table S1). The estimated heterozygote to homozygote ratios were 1.4–1.7 for short indels (Table 2, Table S1), also comparable to those calculated from a previous study [11]. Multiple SNPs and/or indels that are in close proximity (i.e., less than two reference bases between two variant sequences) were grouped together as a single variant locus, termed block substitutions. In total, we identified 394,517 distinct block substitutions in our sample, with each individual genome carrying approximately 71,398 such block substitutions (Table 2, Table S1). These block substitutions can be length-conversing (the sample's allele and reference are the same length) or length-changing, and the change in sequence length caused by these substitutions ranges 0–183 bp (Figure S3).

Table 2

Summary of identified SNPs and indels/block substitutions.

Variant Type	Total No. of Variants	Average No. of Variants in individual genomes
SNP	10,871,465	3,307,678
Intergenic	6,674,155	2,054,900
Intragenic	4,197,310	1,252,778
Intron	3,473,672	1,043,427
UTR	614,753	181,267
Splicing acceptor site	6,644	1,898
Splicing donor site	1,547	398
Coding domain	79,762	18,796
Synonymous	36,538	9,612
Non-synonymous	43,224	9,184
Missense	42,549	9,082
Nonsense	616	87
Nonstop	59	15
Indels	2,815,215	421,088
Coding domain	3,606	381
Frameshift	2,344	217
Frameshift-preserving	1,262	164
Block substitutions	394,517	71,398
Coding domain	1,870	274
Synonymous	30	20
Frameshift	264	17
Missense	1,546	234
Nonsense	29	3
Nonstop	1	<1

Open in a separate window

Of all the SNPs (n=10,862,507) mapped to autosomes and X chromosome, 7,658,805 (70.5%) were already present in dbSNP (v131), and 9,134,659 (84.1%) were in the 1000 Genomes Project Phase 1 data set (released on 5/21/2011, http://www.1000genomes.org/), whereas the rest 1,390,686 (12.8%) SNPs were in neither dataset and thus considered to be novel (Fig. 1a). In contrast, a large proportion (46.2%) of the identified short indels was novel relative to dbSNP (v131) and the 1000 Genomes Project Phase 1 data set (Fig. 1b). Interestingly, short indels reported in any of the three data sets were largely unique to that data set, showing only ~20–40% overlap with any other data set (Fig. 1b). The limited overlap of indels identified between different data sets or between different individual genomes has been previously recognized [11], [12], [15], [25]. This observation may reflect the substantial challenges in making accurate calls of indels. For instance, the estimated false discovery rates for indels were fairly high by both the Complete Genomics approach (~3.0–6.5%) [8] and the 1000 Genomes project approach (~3.7–8.1%) [29]. Moreover, 18% of indel sites identified by the 1000 Genomes project were found to be inconsistent or ambiguous [29].These results highlighted the need for developing advanced experimental and analytic methods for identification and characterization of indels.

Open in a separate window

Figure 1

Summary characterizations of the identified variants.

A–B, Venn diagram showing SNPs and indels identified in the present study overlapping with those archived in the dbSNP (v131) and the 1000 Genomes Project Phase 1 data sets (released on 5/21/2011). To account for differences in placement of many indels between different data sets, indels were considered to match if they were within 25 bp distance and of the same size. Only SNPs and indels mapped to autosomes and X chromosome were analyzed. C, Genome-wide distribution of novel SNPs. Total number of novel SNPs (compared to dbSNP v131 and the 1000 Genomes Project pilot phase) were calculated in non-overlap 1-megabases (Mb) windows across the human genome and plotted in ideograms using Idiographica [41]. The diversities were illustrated by colors, with red indicating higher numbers or proportions and blue indicating lower numbers or proportions. Genomic regions in which no SNPs were identified or no reference sequences could be determined are shown in grey. D, Allele frequency spectrum of novel SNPs.

Consistent with previous reports [12], [13], high densities of variations were observed at the human leukocyte antigen (HLA) region, and a number of telomeric and sub-telomeric regions (Fig. S1). In general, we observed strong correlations in genomic distributions between distinct forms of variants (Fig. S2). Specifically, short insertions and deletions showed very strong correlation throughout the entire genome, whereas block substitutions were relatively poorly correlated with other forms of variants (Fig. S2).

The genome-wide distribution of novel SNPs (Fig. 1c) exhibited similar patterns as that of the total SNPs (Fig. S1). As expected, the vast majority (~96%) of novel SNPs represented low-frequency (minor allele frequency, MAF, ≤5%) variants (Fig. 1d), as most common SNPs in human populations were already included in dbSNP. Intriguingly, several regions contained relatively high proportions of novel SNPs (Fig. S1). For instance, 75–83% of all SNPs identified at centromeric regions on chromosomes 3, 7, 16 and X, were not in the dbSNP (Table S2). In addition, similarly high fractions of novel SNPs were also observed in these regions, when comparing our data with the 1000 Genomes Project Phase 1 data set (novel rate: 80–97%) and when comparing the 1000 Genomes Project Phase 1 data set with the dbSNP (novel rate: 84–91%) (Table S2). The high discrepancy between SNPs from different data sets that are mapped to these regions may suggest that it is difficult to accurately sequence these regions with currently available sequencing techniques.

To evaluate the accuracy of our data, we first estimated the SNP genotyping concordance between whole genome sequencing and SNP arrays for 20 subjects that had been genotyped with the Affymetrix Genome-Wide Human SNP Array 6.0. Out of the 906,600 SNPs called on the SNP arrays, 801,557 (98.5%) were successfully called by the present sequencing project and the genotype concordance rates were 99.83% and 99.45% for homozygous and heterozygous SNPs, respectively. By comparing with the recently reported 10 Korean genomes [14], the 69 genomes data set from Complete Genomics Inc (http://www.completegenomics.com/sequence-data/download-data/), and an exome sequencing data of 100 Han Chinese subjects (Lin et al., unpublished), we in silico validated 177,554 SNPs and 411,658 indels that were not archived in the dbSNP (v131) or the 1000 Genomes Project Phase 1 data sets, supporting that these novel variants are likely to be genuine polymorphisms and not just unique to our specific population. To further evaluate the accuracy of the novel variants, we tested 25 novel nonsynonymous SNPs by traditional Sanger sequencing. Twenty-four SNPs out of the 25 were validated, highlighting the high quality and reproducibility of our data.

Gene-based characterization of SNPs and indels

Of all the identified SNPs, 79,762 were mapped to CDS (coding sequences) of RefSeq transcripts, including 36,538 synonymous, 42,549 missense, 616 nonsense, and 59 nonstop (causing the loss of a stop codon) SNPs (Table 2). Consistent with previous reports [12], [13], we observed a significantly higher fraction of low-frequency variants (MAF≤0.05) in non-synonymous SNPs, when compared to synonymous or noncoding SNPs, presumably reflecting purifying selection against deleterious mutations (Fig. S3). Of all the identified indels, only a very small fraction (0.13%) is mapped to coding sequence regions and most of the coding indels were predicted to cause frameshift changes (Table 2), which is consistent with previous reports [12], [13]. We observed a noticeable preference for 3n-bp indels in coding regions compared to indels in non-coding genic regions, such as in introns and UTRs (Fig. S3). This preference has also been observed in other studies [9], [13], [30] and may reflect the selection pressure for minimally impacting variants in coding regions. In addition, we detected 1,870 block substitutions that are mapped to coding sequence regions, and almost all of the coding block substitutions were predicted to cause non-synonymous changes (Table 2). This is expected, as block substitutions involve multiple variant bases and thus are more likely to result in amino acid changes. Similar bias towards non-synonymous variants were also reported in a previous study by Rosenfeld et al [31], who focused on a specific type of block substitution, namely, double nucleotide polymorphisms (DNPs) and triple nucleotide polymorphisms (TNPs), which by definition are length-conserving block substitutions involving two or three consecutive nucleotides.

We performed GO (Gene Ontology) analyses to determine whether genes containing high density of amino acid changing variants were enriched in specific biological processes. The GO terms related to sensory perception, such as ‘sensory perception of chemical stimulus’ (p=3.47×10⁻⁸) and ‘sensory perception of smell’ (p=7.21×10⁻⁷) showed significant enrichment for genes with high density of amino acid changing variants (Table 3). In contrast, amino acid changing variants were significantly under-represented in multiple biological processes related to ‘cellular macromolecule metabolic process’ (p=1.38×10⁻¹⁹) and ‘RNA metabolic process’ (p=1.52×10⁻¹⁷) (Table 3). These results were not unexpected because these metabolic and signal transduction related biological processes are critical for cell survival and functioning, consequently, are expected to be under strong purifying selection against deleterious mutations. On the other hand, sensory perception related genes are known to be highly associated with environmental adaptation. Thus, variants may have accumulated that differ in frequencies among human populations, suggesting balanced selection, the possible relaxation of purifying selection, and/or an increased mutation rate [9], [15].

Using the Variant Annotation Tool (version 2.0.1.) [32] with the GENCODE v7 annotation [33], we examined SNPs and indels for variants predicted to result in the complete loss-of-function (LoF). LoF variants were defined as SNPs predicted to create or disrupt a stop codon, frameshift indels, and variants predicted to disrupt a splice site. On average, each person carried approximately 179 LoF variants with ~44 of these variants occurred in a homozygous state, which would result in complete “knock-out” of the annotated genes. Altogether, we found 182 unique “knocked-out” genes in our samples. While about 40% of these genes were “knocked out” in just one or two individuals, 46 genes were “knocked out” in over 30% of our samples (Fig. 2), suggesting that a number of genes are commonly “knocked-out” in general populations. GO analysis suggested that these “knocked-out” genes were significantly enriched in several biological processes, particularly those related to antigen processing and immune responses (Table 3). Intriguingly, 8 LoF variants located in 8 genes (Table S3) were found to be homozygous in all the 44 genomes, as well as in other independently sequenced human genomes [26], [30], suggesting that the reference sequences at these sites may represent rare alleles or sequencing errors.

Open in a separate window

Figure 2

Identification of “knocked-out” genes.

A, Frequency spectrum of observed “knocked-out” genes. Genes containing homozygous LoF variants were expected to be silent or knocked-out. Numbers of “knocked-out” genes were counted with respect to the frequency of “knock-out” occurrence in the 44 genomes.

Disease-associated SNPs and indels

Using Trait-o-matic (https://github.com/xwu/trait-o-matic), we found that each individual genome carried an average of 90 non-synonymous SNPs that were known to be associated with OMIM diseases/traits. Among the 44 samples, 29 subjects carried 14 heterozygous variants causing autosomal dominant diseases/traits, and 4 subjects were homozygous for 4 variants causing autosomal recessive diseases/traits (Table S4). However, the 44 subjects were all apparently healthy with no recognizable or self-reported diseases at the time of enrollment. This might be explained by the fact that the majority of these variants were known to have reduced/age-dependent penetrance (e.g., R468H mutation in MFN2 gene for Charcot-Marie-Tooth disease-2A2) and/or did not cause severe health-threatening symptoms (e.g., wet/dry ear wax). This perspective is supported by the fact that 8 of the 18 variants were present in low to intermediate frequency in general populations (Table S4). Alternatively, some of these variants might have been erroneously assigned as disease mutations. Each genome also contained approximately 276 additional non-synonymous SNPs, including about 34 potential nonsense/nonstop SNPs and 11 frameshift indels and block substitutions in genes associated with OMIM diseases/traits.

Mitochondrial and Y chromosome analyses

By comparing to the revised Cambridge Reference Sequence (rCRS) of the Human Mitochondrial DNA [34], we identified a total of 306 mitochondrial variants, including 285 SNPs, 13 short indels, and 8 block substitutions (Table S5). Approximately 54% of these SNPs, and almost all the indels and block substitutions, were novel relative to dbSNP (v131). In the 37 mitochondrial genes, we identified 54 non-synonymous SNPs (53 missense and 1 nonstop), and one frameshift block substitution.

A total of 8,673 SNPs, 1,352 indels, and 998 block substitutions were identified on the Y chromosome (Table S5). Similar to the mitochondrial genome, approximately 54% of these SNPs and the majority (~65%) of indels/block substitutions on the Y chromosome were novel relative to dbSNP (Table S5).

Copy number variants (CNVs)

Through the analysis of read depth in 2-kb sliding windows, we discovered an average of 3,749 CNVs (size range of 2–154 kb) totaling about 12 Mb in each genome. Approximately 45% of the CNVs detected in our study overlapped (≥50% of sequence overlap) with CNVs currently annotated in the Database of Genomic Variants (DGV) [35]. Similar to previous reports [9], [20], we observed a steadily decreasing number of events with increasing CNV length (Fig. S4). Out of the 493 CNVs identified in the 20 subjects with the Affymetrix Genome-Wide Human SNP Array 6.0, 255 overlapped with CNVs called by the sequencing data.

Genomic DNA preparation and sequencing

Genomic DNA was isolated from peripheral blood using the Gentra Puregene Blood kit (Qiagen, Valencia, CA) according to the recommended protocol. DNA concentration was measured by Nanodrop 1000 (Thermo Scientific, Wilmington, DE) and Quant-iT™ Pico Green dsDNA kit (Invitrogen, Carlsbad, CA). Whole genome DNA sequencing was performed by Complete Genomics, Inc (Mountain View, CA), using its paired end library preparation and sequencing-by-ligation methodology as established previously [8], [16].

The resulting mate-paired reads were initially mapped to the NCBI reference genome (build 37.1, GRCh37/hg19) using a fast algorithm, and these initial mappings were both expanded and refined by local de novo assembly, which was applied to all regions of the genome that appeared to contain variations based on initial mappings [8]. Mapped reads were assembled into a best-fit diploid sequence with two separate resultant sequences for each locus in diploid regions (exceptions: mitochondria were assembled as haploid, and for males the non-pseudo-autosomal regions were assembled as haploid) by using a custom software suite that implements both Bayesian and de Bruijn graph techniques [8]. Variants were called by independently comparing each of the diploid assemblies to the reference. Data for each genome were delivered as lists of sequence variants (SNPs, short indels, and block substitutions) relative to the NCBI reference genome accompanied with variant quality scores. Block substitutions were called where a series of nearby reference bases had been replaced with a different series of bases in an allele. Block substitutions can be length-conserving (the same number of bases as the corresponding reference sequence region) or length-altering. For sites with multiple variant bases, if at least two reference bases on both alleles were called between two variant sequences, then the sites were broken into smaller variant events. The sequence variant data from this study are available at http://tulane.edu/publichealth/bio/genetics-and-genomics-study-of-osteoporosis.cfm.

CNV analysis

CNVs were identified by using CNV-seq (http://tiger.dbs.nus.edu.sg/cnv-seq/) [39], which compares the numbers of mapped reads in a sliding window between two individuals and assesses the hypothesis of no copy number variation based on a probabilistic model. For this analysis, one of the sequenced genomes was randomly picked as the reference genome and the remaining 43 genomes were compared with this reference genome for numbers of mapped reads in 2 kb sliding windows.

Variant annotation

Based on their locations mapped to the RefSeq transcripts, the identified variants were classified into different groups, including intergenic, CDS, intron, splicing donor and acceptor sites, and 5′- /3′- untranslated regions (UTRs). Since variants may be fall into multiple categories depending on the alternative transcript isoforms examined, a hierarchy CDS>UTRs>splicing sites>intron>intergenic was applied, such that a variant was only counted once at its highest level in the hierarchy. SNPs and indels mapped to CDS and splicing sites were further assessed for their potential effects on the gene products to identify potential LoF variants that are predicted to create or disrupt a stop codon, shift the normal reading frame, and disrupt a splice site. The identification of LoF variants was performed by using the Variant Annotation Tool (version 2.0.1.) [32] with the GENCODE v7 annotation reference [33].

The program Gorilla (http://cbl-gorilla.cs.technion.ac.il/) [40] was employed to test whether genes containing high or low density of amino acid changing variants (non-synonymous SNPs and frameshift indels/block substitutions) were enriched in certain GO terms (biological processes). The density of deleterious coding variants for each Refseq gene was calculated by dividing the total number of deleterious variants in the corresponding Refseq transcripts by the sum of CDS lengths of the transcripts. A ranked list of genes was generated by sorting the Refseq genes in three steps, first based on the density of deleterious coding variants, then by the total number of deleterious variants, and finally by the sum of CDS lengths. The ranked list of genes was then uploaded into Gorilla for GO enrichment analysis.

Non-synonymous SNPs identified in these personal genomes were screened for known associations with inherited disorders/traits by using Trait-o-matic (https://github.com/xwu/trait-o-matic), which finds and cross-references the non-synonymous SNPs with records in the Online Mendelian Inheritance in Man (OMIM) database (http://www.omim.org).

Variant correlation

The human genome was divided into consecutive 300-kb windows and the numbers of different types of variants (SNP, insertion, deletion, and block substitution) were counted within each window. Pairwise Pearson's correlation coefficients between different types of variants were computed using variant counts from 100 consecutive windows through the cor function in R.

We thank Aaron Solomon, Erica Beilharz, and Zuoming Deng at Complete Genomics Inc. for their assistance in setting up and performing the whole genome sequencing, and providing consultation for data analysis.