Kidd et al.

BMC Genomics 2014, 15:262

http://www.biomedcentral.com/1471-2164/15/262

RESEARCH ARTICLE Open Access

Exome capture from saliva produces high quality

genomic and metagenomic data
Jeffrey M Kidd1,2†, Thomas J Sharpton3,4†, Dean Bobo5, Paul J Norman6, Alicia R Martin1, Meredith L Carpenter1,
Martin Sikora1, Christopher R Gignoux7, Neda Nemat-Gorgani6, Alexandra Adams1, Moraima Guadalupe8,
Xiaosen Guo9, Qiang Feng9, Yingrui Li9, Xiao Liu9, Peter Parham6, Eileen G Hoal10, Marcus W Feldman11,
Katherine S Pollard3,12, Jeffrey D Wall12, Carlos D Bustamante1 and Brenna M Henn1,5*

Abstract
Background: Targeted capture of genomic regions reduces sequencing cost while generating higher coverage by
allowing biomedical researchers to focus on specific loci of interest, such as exons. Targeted capture also has the
potential to facilitate the generation of genomic data from DNA collected via saliva or buccal cells. DNA samples
derived from these cell types tend to have a lower human DNA yield, may be degraded from age and/or have
contamination from bacteria or other ambient oral microbiota. However, thousands of samples have been previously
collected from these cell types, and saliva collection has the advantage that it is a non-invasive and appropriate for a
wide variety of research.
Results: We demonstrate successful enrichment and sequencing of 15 South African KhoeSan exomes and 2 full
genomes with samples initially derived from saliva. The expanded exome dataset enables us to characterize genetic
diversity free from ascertainment bias for multiple KhoeSan populations, including new exome data from six HGDP
Namibian San, revealing substantial population structure across the Kalahari Desert region. Additionally, we discover
and independently verify thirty-one previously unknown KIR alleles using methods we developed to accurately map
and call the highly polymorphic HLA and KIR loci from exome capture data. Finally, we show that exome capture of
saliva-derived DNA yields sufficient non-human sequences to characterize oral microbial communities, including
detection of bacteria linked to oral disease (e.g. Prevotella melaninogenica). For comparison, two samples were
sequenced using standard full genome library preparation without exome capture and we found no systematic
bias of metagenomic information between exome-captured and non-captured data.
Conclusions: DNA from human saliva samples, collected and extracted using standard procedures, can be used to
successfully sequence high quality human exomes, and metagenomic data can be derived from non-human reads.
We find that individuals from the Kalahari carry a higher oral pathogenic microbial load than samples surveyed in
the Human Microbiome Project. Additionally, rare variants present in the exomes suggest strong population structure
across different KhoeSan populations.
Keywords: Exomes, KhoeSan, Genetic diversity, Metagenomics, Microbiome

* Correspondence: [email protected]
†
Equal contributors
1
Department of Genetics, Stanford University, Stanford, CA 94305, USA
5
Department of Ecology and Evolution, Stony Brook University, Life Sciences
Bldg, Room 640, Stony Brook, NY 11794, USA
Full list of author information is available at the end of the article

© 2014 Kidd et al.; licensee BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative
Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and
reproduction in any medium, provided the original work is properly credited. The Creative Commons Public Domain
Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article,
unless otherwise stated.
Kidd et al. BMC Genomics 2014, 15:262 Page 2 of 17
http://www.biomedcentral.com/1471-2164/15/262

Background disequilbrium (LD) of any population surveyed and thus

Sampling of saliva or via buccal cell extractions is a the largest effective population size [2]. However, in order
widely employed, non-invasive method of collecting hu- to test hypotheses regarding population sub-structure, nat-
man DNA for both biomedical and ancestry experiments. ural selection and biomedically relevant variants in Africa,
DNA extracted from saliva fluid has been used on single it is essential to have both large sample sizes and genomic
nucleotide polymorphism chip arrays, methylation arrays, data that are un-biased with regard to ascertainment
targeted resequencing, exome, and whole genome sequen- schemes.
cing [1-7]. However, the low total yield of DNA from a
single sample and the presence of many non-human DNA Results
fragments make next-generation sequencing of saliva sam- Fifteen human saliva samples were selected for exome
ples impractical for some applications. Targeted enrich- sequencing. Samples were split into two batches (“Pilot 1”
ment strategies, such as hybridization methods designed and “Pilot 2”), representing samples enriched using the
to capture the exons of annotated genes (the ‘exome’) Agilent SureSelect 50 Mb human All-Exon design and se-
prior to sequencing, offer a way to circumvent some of quenced with the Illumina GAII machine and a replication
the limitations posed by saliva-derived DNA samples. We batch enriched using the Agilent SureSelect 44 Mb human
demonstrate the successful sequencing of multiple human All-Exon design and sequenced using Illumina HiSeq. We
exomes from saliva-derived samples using commercially included a familial quartet with two daughters (Family 1),
available reagents for exome capture. an extended pedigree of first cousins and half-siblings
Exome sequencing and other capture methods permit (Family 2), and eight purportedly unrelated individuals
the high-coverage sequencing of a small portion of the (Additional file 1: Figure S1). Family 1 displayed com-
genome. This approach represents a trade off between plex ancestry from KhoeSan, European and both eastern
depth of coverage vs. breadth of the genome that is in- and western African populations (see [2]). Family 2 and the
terrogated, and has the potential to revolutionize gen- un-related individuals self-reported their ancestry as being
omic medicine [8,9]. In addition to direct applications to from only KhoeSan populations (Nama- or N|u-speakers).
human disease, exome sequencing of a modest number We obtained 3-25 ug total DNA from each saliva sample.
of individuals can reveal important aspects of human Each aliquot was processed using the Agilent SureSelectXT
evolution [10-12]. The capability to apply these approaches library preparation kit followed by enrichment with the
to DNA derived from saliva, which is more easily obtained SureSelect 44 Mb or SureSelect 50 Mb human All-Exon
and less invasive than blood or other tissue collection, will capture probes. Using standard Illumina post capture bar-
greatly facilitate the detailed examination of genetic vari- codes, libraries were sequenced on either an Illumina GAII
ants that may be associated with specific traits or have ex- or HiSeq machine. Aliquots from two samples (SA1000
perienced adaptive evolution [13,14]. and SA1025) were also sequenced without exome capture,
We focus on a unique set of DNA samples from the using the Illumina TruSeq library preparation kit (SA1000)
≠Khomani KhoeSan of South Africa to illustrate the util- and the Illumina Nextera library preparation kit (SA1025).
ity of exome sequencing via saliva. African genetic diver- The whole genome sequence (WGS) libraries were then
sity remains poorly understood, in part because many sequenced on two lanes of an Illumina HiSeq.
regions of the continent lack adequate healthcare infra-
structure, which can make blood collection impractical. Sequencing statistics
The indigenous KhoeSan peoples of southern Africa are An average of 76.5 million 75 bp paired reads and 84.3
a collection of hunter-gatherer and pastoralist groups million 100 bp paired reads were obtained for each indi-
who speak “click languages”, classified into three distinct vidual in the Pilot 1 GAII and Pilot 2 HiSeq exome ex-
language families. The genetic diversity of these, and periments (Table 1). Across all samples, 86.8%-98.1%
related populations, remains under-ascertained. The of the reads mapped to the human genome reference
genome of one Tuu-speaking San (“!Gubi”) has been (GRCh37) (Figure 1). On average, ~70-75% of non-
fully sequenced and found to contain over 700,000 duplicate, mapped reads fell in the specified target re-
novel polymorphisms [15]. Gronau et al. showed that this gions. This on-target percentage is similar to previous
San genome was highly divergent among known genomes, on-target percentages (70-87%) for standard blood or
even compared to other African individuals [16]. They esti- cell line-derived human DNA with Agilent SureSelect
mated the population divergence between western African exon designs [17,18].
individuals and the San to be about 110,000-130,000 years Two samples (SA006 and SA035) displayed a high per-
ago, over twice as old as the divergence between western centage of duplicate reads (54% and 78%) (Additional
Africans and Eurasians. Additionally, single nucleotide file 1: Figure S2, Table 1). To understand whether SA006
polymorphism (SNP) array data demonstrated that the and SA035 had high duplicate rates due to low human
≠Khomani San population had the lowest levels of linkage DNA input or whether there were other issues with
 http://www.biomedcentral.com/1471-2164/15/262
Kidd et al. BMC Genomics 2014, 15:262
Table 1 Summary statistics for KhoeSan exomes
Total Unmapped % Un-mapped % PCR % Mapped on Median target % of variants Autosomal Autosomal Non-ref.
readsa reads reads duplicates target coverageb coveredc SNV singletons concordanced
Pilot 1 SA006 69,272,282 9,122,731 13.2% 54.2% 63.5% 12 94.9% 25,225 657 0.9897
SA008 113,888,276 2,143,408 1.9% 19.8% 78.2% 73 99.5% 26,408 955 0.9947
SA011 78,006,472 1,664,959 2.1% 33.7% 77.4% 40 99.0% 26,365 67 NA
SA012 67,209,032 1,353,187 2.0% 20.5% 75.7% 42 99.3% 26,722 86 NA
SA035 85,142,498 5,812,851 6.8% 78.0% 79.4% 10 92.1% 24,692 1,726 0.9884
SA051 76,076,464 3,102,819 4.1% 27.8% 76.5% 37 98.8% 27,674 1,239 NA
SA052 60,375,472 1,247,951 2.1% 12.9% 78.2% 41 98.8% 27,779 755 0.9968
SA054 62,358,148 1,959,032 3.1% 27.9% 73.9% 31 99.3% 28,024 817 0.9956
Pilot 1 mean 76,541,081 3,300,867 4.4% 34.4% 75.4% 35.75 97.7% 26,611 788e 0.9930
Pilot 2 SA1000 77,069,730 8,387,491 10.9% 9.5% 57.3% 44 98.4% 27,921 2,483 0.9915
SA1001 85,479,934 3,551,500 4.2% 11.4% 74.2% 67 98.7% 27,694 2,318 0.9939
SA1002 92,542,846 4,674,919 5.1% 15.5% 70.1% 65 98.8% 27,886 3,286 0.9941
SA1006 83,545,692 4,002,665 4.8% 18.1% 74.5% 59 98.4% 27,446 2,442 0.9927
SA1010 87,939,484 4,445,502 5.1% 14.5% 71.0% 62 98.6% 27,295 1,782 0.9935
SA1011 82,377,158 7,810,714 9.5% 11.6% 49.2% 40 98.5% 27,484 2,717 0.9887
SA1025 81,405,650 2,498,412 3.1% 10.0% 87.8% 63 99.3% 28,696 2,676 0.9934
Pilot 2 mean 84,337,213 5,053,029 6.1% 12.9% 69.2% 57.14 98.7% 27,775 2,529 0.9925
a
Total number of DNA fragments including: mapped, unmapped and duplicate reads.
b
Limited to non-duplicate reads on autosomes, as calculated by GATK Unified Genotype.
c
Limited to XX autosomal SNPs identified at the 99% VQSR threshold.
d
Concordance at heterozygous and homozygous non-reference positions as compared to Illumina OmniExpress or 550K.v2 SNP arrays.
e
Fewer average singletons as a result of including closely related individuals in Pilot 1. See Additional file 1: Table S1 for individual data.

Page 3 of 17
Kidd et al. BMC Genomics 2014, 15:262 Page 4 of 17
http://www.biomedcentral.com/1471-2164/15/262

Figure 1 Schematic of mapping and calling pipelines. Each box summarizes the data and data format used for each step of the human
exome and microbiome mapping/calling pipelines. The pipeline begins with next-generation sequencing raw reads obtained from exome
sequencing of saliva-derived DNA and ends in finalized exome variant calls and microbiome taxonomic abundances. Arrows indicate analysis
methods used to process the human and saliva microbiome data (see Methods).

read data, we examined the distribution of mapping effective coverage and ~80% of reads with mapping
quality for all uniquely mapped reads for each sample. qualities ≥30. This difference is unlikely to be due to di-
These two samples had the lowest numbers of mapped vergence from the reference because we observed no sys-
reads and the lowest proportion of reads with mapping tematic differences in mapping quality metrics between
qualities ≥ 30 (35.2% and 68.6%, respectively, Additional the European- and Bantu- admixed Family 1 and the
file 1: Figure S3). The remaining Pilot 1 samples had higher KhoeSan Family 2. Due to lower mapping rates, SA006
Kidd et al. BMC Genomics 2014, 15:262 Page 5 of 17
http://www.biomedcentral.com/1471-2164/15/262

and SA035 displayed overall lower mapped coverage than excluding the two daughters in Family 1 (Table 1). We
the other samples. However, 90% of target sites were cov- computed genotype concordance for 12 individuals
ered at a depth of at least 10x for all individuals except (sufficient DNA was not available for SA011, SA012,
SA006 and SA035 (Additional file 1: Figure S2). The SA051) based on data from the Illumina OmniExpress
average percent of unmapped reads was higher for saliva- or 550 K.v2 SNP arrays [2]. Non-reference (NR) con-
derived exomes compared to six HGDP San samples se- cordance, that is concordance only at heterozygous or
quenced using DNA obtained from cell lines (Additional non-reference homozygous genotypes, was calculated
file 1: Table S1). However, the primary difference in se- using GATK [24,25] and concordance exceeded 98%
quencing efficiency between saliva- and cell-line derived for all individuals genotyped.
DNA results from differences in the mean rate of dupli-
cate reads: Pilot 1, 34.4%; Pilot 2, 12.9%; HGDP, 9.8%. Novelty compared to 1000 genomes project
Pilot 1 likely has a higher duplicate rate due to lower We compared 13 KhoeSan exomes from our study
DNA quality (see below). (excluding children SA011 and SA012), to exomes se-
quenced as part of the 1000 Genomes Project (1000G)
Read quality [23], HGDP Namibian San, and San Nimblegen exomes
We hypothesized that the difference in mapping quality from Schuster et al. [15]. We chose three populations
among samples could be due to different levels of DNA of African ancestry for comparison: ASW, African-
damage. To test this hypothesis, we analyzed the distri- Americans from the Southwestern United States; LWK,
bution of mismatches along the reads by comparing each Luhya from Kenya; YRI, Yoruba from Nigeria; and GBR,
read to the human reference sequence after mapping. If from Great Britain to represent European ancestry. Since
the genomic DNA had been degraded before shearing, the 1000G dataset contained many more individuals than
for example due to variable storage conditions, one our KhoeSan dataset, these populations were randomly
would expect an increase in mismatches at the ends of down-sampled to 13 individuals for comparison. We note
the reads; specifically, an excess of thymines at the 5′ that these disparate datasets were processed using differ-
end of the read and an excess of cytosines at the 3′ end ent pipelines, in some cases involving multiple-sample
of the read, similar to what is seen in ancient DNA calling and imputation with a large number of other
[19-21]. However, for SA006 and SA035 we observe an exomes, with varying degrees of coverage and sample re-
increased rate for all types of substitutions at the begin- latedness. Between 28,000-29,000 variants appear to be
ning of the reads, with the highest rates for those to- common to all 5 populations (i.e. between 38%-53% are
wards the purines G and A (Figure 2). Reads from shared in each three-way comparison) (Figure 3). The
SA006 in particular show a pronounced increase, for ex- South African ≠Khomani San appear comparable to the
ample a ~10 fold increase in the rate of A - > G substitu- Yoruba and Luhya populations in terms of the number of
tions in the first position compared to positions further private SNPs yet share slightly more variants with GBR
along the read. This pattern is absent from all other than either other sub-Saharan African population. This
samples in Pilot 1, with the exception of SA051, which reflects the degree of recent European admixture in
also shows a slight increase at the first base (Additional the ≠Khomani. In order to compare two KhoeSan popu-
file 1: Figure S4). We also observe an overall increase in lations, we used 6 Namibian HGDP San exomes se-
substitution rates towards the end of the reads, which is quenced on the same Agilent Platform to >70x coverage
shared across all samples and consistent with the in- [Martin et al, in preparation] and included African-
creased rate in sequencing error with increasing number Americans in an attempt to control for recent gene flow
of sequencing cycles. The pattern of mismatch rates from Bantu-speaking and European groups (Figure 3C).
does not support a hypothesis of simple degradation. The Namibian and South African populations share
only 4,692 unique variants that are not also found in the
Genotype and variant statistics ASW. This may reflect the small sample size, or that the
Variants were called using the Genome Analysis Tool KhoeSan populations in the north vs. south Kalahari re-
Kit (GATK) and selected using the Variant Quality Score main highly differentiated [26,27]. Given the high con-
Recalibration (VQSR) procedure with cutoffs set such cordance between our exome sequence data and the
that 99% of variants also found in the 1000 Genomes Illumina SNP array, we believe that the high genetic diver-
Omni2.5 and HapMap3 SNP training set were retained sity of the South African exomes is not an artifact caused
[22-24]. We identified 82,093 variants, with a transition/ by high false positive rates.
transversion ratio of 3.14. On average, within the target
regions, each individual had a genotype call at 98% of Population differentiation of the KhoeSan
sites variable in the 15 sample dataset (Table 1). Single- We performed principal component analysis (PCA) on
ton counts varied from 657 to 3,286 autosomal sites, the unrelated ≠Khomani KhoeSan (Additional file 1:
Kidd et al. BMC Genomics 2014, 15:262 Page 6 of 17
http://www.biomedcentral.com/1471-2164/15/262

Figure 2 (See legend on next page.)

Kidd et al. BMC Genomics 2014, 15:262 Page 7 of 17
http://www.biomedcentral.com/1471-2164/15/262

(See figure on previous page.)

Figure 2 Assessment of base substitutions from mapped reads. Each mapped read was compared to the genome reference sequence
to assess patterns consistent with DNA degradation. At each of the 75 positions along a read, we plot the frequency of substitution types,
for both the forward (left) and reverse (right) reads from each read-pair. Analysis was limited to 1 million reads from chromosome 1; all raw
reads are plotted. Three individuals with varying levels of substitution errors are shown: (A) SA006 with overall higher substitution rate and
an excess of purines at the start of the first read, (B) SA035 with a slightly elevated substitution rate and excess of purines at the start of
the first read, and (C) SA054 with a low substitution rate and no bias at the beginning of the first read. The additional five Pilot 1 individuals
tended to resemble SA054 (Additional file 1: Figure S4). Removal of reads with any soft-clipping substantively reduced the mis-incorporation
rate for SA006 and SA035.

Figure S1), Schuster et al. [15] Namibian San, and HGDP HLA and KIR
Namibian San exomes, along with several different popu- The HLA and KIR loci include some of the most poly-
lations from the 1000 Genomes Project using smartpca morphic genes in the human genome and are function-
from the EIGENSOFT software package [28] (Figure 4, ally involved in the immune system and reproduction
Additional file 1: Figure S6). PCA clearly differentiates [29,30]. Contributing to HLA and KIR polymorphism are
the populations included in this study. PC1 separates inter-locus recombination and gene duplication, factors
the African populations from the European population rendering these loci difficult to analyze with genomic-
(GBR). PC2 separates the populations of western African scale data, but among the most stringent for assessing its
ancestry (LWK and YRI) from the southern African pop- validity. We analyzed the three highly polymorphic HLA
ulations (SSAN, NSAN, and HGDP San). PC3 sepa- class I genes, HLA-A, -B and -C (6p21), and the KIR locus
rates the northern from the southern Kalahari KhoeSan (19q13.4), which has variable content of four to thirteen
populations, suggesting there is substantial substructure polymorphic genes. Despite using a highly conservative
among these groups. PC4 separates a single SSAN indi- strategy to remove read-pairs that did not map exclusively
vidual from both the HGDP San and the rest of the to one of the targeted loci, genotypes were obtained for
≠Khomani KhoeSan individuals. This individual belongs 4,070 HLA class I and KIR SNPs for the fifteen indi-
to Family 2, so the PCA was revised to include the viduals studied (Tables 2 and 3, Additional file 1: Table S2,
remaining two Family 2 individuals (1st cousin/half- Additional file 1: Table S3). Sufficient read-depth (at least
sibling relationships, Additional file 1: Figure S1) in order 20 for homozygous positions and 10 for heterozygous
to assess whether SA051 was an outlier, however PC4 still positions) was obtained for determination of all the HLA
strongly separated all Family 2 individuals from other class I and KIR alleles present, with exception of HLA-A
≠Khomani (Additional file 1: Figure S6). and -B from individual SA006. Fourteen of the individuals

Figure 3 Novelty compared to 1000 genomes project. We compared the number of nonreference variants in the South African KhoeSan
[SSAN] with (A) 1000 Genomes Yoruba samples [YRI], (B) eastern African Luhya [LWK], and (C) Namibian San from HGDP and African-Americans
[ASW]. Sites that were included in this analysis required the presence of genotype information for at least 95% of the individuals in the joint dataset.
The exome data from ASW and LWK was derived from 1000 Genomes Project, Phase 1 – March 17, 2012 release, from which 13 individuals were
randomly sampled. The Venn diagram illustrates the number of shared and unique nonreference variants among populations. The Vennerable package
in R was utilized for plotting purposes.
Kidd et al. BMC Genomics 2014, 15:262 Page 8 of 17
http://www.biomedcentral.com/1471-2164/15/262

A B C

0.3
ASW ASW
0.10

GBR GBR

0.2
LWK LWK

0.2
YRI YRI
0.05

SSAN SSAN
NSAN NSAN

0.0
0.1
HGDP HGDP
0.00

PC3 [1.27%]
PC2 [3.74%]

PC4 [1.12%]
0.0

−0.2
−0.1
−0.10

ASW

−0.4
GBR

−0.2
LWK
YRI
SSAN

−0.6
−0.20

−0.3
NSAN
HGDP

−0.2 −0.1 0.0 0.1 −0.20 −0.15 −0.10 −0.05 0.00 0.05 0.10 −0.3 −0.2 −0.1 0.0 0.1 0.2 0.3

PC1 [5.76%] PC2 [3.74%] PC3 [1.27%]

Figure 4 Principal components analysis of 61,000 exonic SNPs in the ≠Khomani San and other African populations. Exomes from 1000
Genomes Phase 1, Schuster et al. [15], and HGDP San were combined with the ≠Khomani San (related samples from Families 1 and 2 were removed).
5.76% of the variance is explained by PC1, 3.74% by PC2, 1.27% by PC3, and 1.12% by PC4. PC1 and PC2 separate Africans from Europeans, and
western Africans from southern Africans, respectively (A). The three KhoeSan populations drive PC3 and PC4 (B and C), supporting prior descriptions of
strong differentiation among Kalahari KhoeSan groups [27], and indicating even sub-structure within the ≠Khomani San samples.

were genotyped using standard methods for HLA class I We hypothesized that these unmapped reads might rep-
and eight for all of the KIR. When comparison of the resent non-human DNA carried through the saliva ex-
methods was restricted to the individuals with high traction. Although we obtained useful results, with high
genome-wide coverage, all but ten SNPs were con- concordance to SNP genotyping arrays, such microbial
cordant with standard genotyping from these samples contamination may contribute to lower effective cover-
(Tables 2 and 3), validating the sensitivity and specifi- age levels. We therefore subjected these unmapped reads
city of our analytical approach (Methods). Moreover, to an independent quality control procedure and used
all 39 discordant SNP genotypes from the two low- a fragment recruitment approach described by Rusch
coverage individuals (SA006 and SA0035) occurred in et al. [31] to identify homologs of non-human reference
clusters of low read depth where only one or neither genomes among a combined pool of 24,139,131 high-
allele was represented. Thus, following stringent filter- quality unmapped reads (Figure 1). To estimate the num-
ing there were no false positive genotypes for any of the ber of species that are detected, we applied a recruitment
HLA class I or KIR SNPs. In total, there were 36 dis- threshold based on the 95% average nucleotide iden-
tinct HLA class I and 91 KIR alleles present, including tity threshold that is commonly used to define microbial
thirty-one previously unknown KIR alleles that were species [32].
discovered by analysis of the exome-sequencing data Across all 15 sequenced exomes, we identified 1,835,400
and independently verified by standard cloning, sequencing high-quality reads (7.6%) that map to the genomes of 1,153
and family study. non-human species. The distribution of the number of re-
cruited reads per genome indicates that a small number of
Saliva metagenomes genomes recruit a large number of reads with most ge-
Although exome capture proved an efficient method of nomes recruiting an insignificant fraction of the reads.
sequencing primarily human DNA, each sample also For example, after normalizing the number of reads re-
contained more than a million unmapped reads (Table 1). cruited per genome by reference genome size, the 100

Table 2 HLA and KIR validation

Exome sequencing Standard genotyping Standard genotyping (excluding SA006 & SA035)
SNPs vs Allelesb Present SNPs vs Concordance SNPs vs Concordance
HG19a HG19 rate (%) HG19 rate (%)
Known Novel
KIR (13 genes) 1469 91 31 955 99.99 670 99.99
HLA class 1 A 690 16 690 99.98 619 100.00
HLA class 1 B 925 12 925 99.99 745 99.99
HLA class 2 C 986 8 986 100.00 814 100.00
a
Non-reference single nucleotide polymorphisms.
b
Unique coding sequences.
Kidd et al. BMC Genomics 2014, 15:262 Page 9 of 17
http://www.biomedcentral.com/1471-2164/15/262

Table 3 HLA and KIR validation for SA006 and SA035 are surprisingly consistent with estimates from traditional
Non-reference SNPS shotgun metagenomic sequencing of saliva communities.
SA006a SA035 Some of the abundant KhoeSan saliva microbiota are
HG19 ref Present Correct Present Correct known contributors to oral disease. For example, Prevo-
A*03:01 37 12 19 19
tella melaninogencia (recruits 5.9% of unmapped reads
after correcting for genome length) is associated with
B*07:02 71 58 34 34
rapidly progressing periodontitis lesions [33]. Similarly,
C*07:02 40 39 50 50 Streptococcus parasanguinis (6.3%) is a primary colonizer
a
It was not possible to obtain the HLA-A and -B genotypes from exome data of human teeth and contributes to dental plaque forma-
of SA006.
tion [34]. Granulicatella elegans (2.7%), an oral com-
mensal associated with infective endocarditis [35], is also
most abundant genomes recruit 98.3% of the reads. found in high abundance among the KhoeSan. We also
Generally, the genomes that recruit the most reads are specifically ascertained the presence of several biomedi-
well-described oral commensal microbiota (Table 4), cally important organisms, some of which may exist at
such as Neisseria subflava, Rothia mucilaginosa, Neisseria relatively low abundance. For example, the Porphyromo-
flavescens, Veillonella dispar, and Prevotella veroralis. The nas gingivalis genome, which represents organisms im-
recruitment of reads across the length of these genomes plicated in periodontal disease and has been linked to
suggests that their detection is not an artifact of a genomic rheumatoid arthritis [36] and heart disease [37], recruits
subsequence that shares similarity with the human gen- a relatively large fraction of reads from all individuals
ome (e.g. Additional file 1: Figure S7.) We verified this by (1.68%). Conversely the Campylobacter rectus genome,
comparing the per genome relative abundance distribution which is also associated with periodontitis [38], recruits
estimated through analysis of these exome-captured meta- a relatively small fraction of reads (0.24%). Only 8 reads
genomes to the corresponding distribution estimated (2.3 × 10-4% of genome length-corrected recruitments)
through analysis of non-capture metagenomes for two were recruited with high fidelity to the genome of Myco-
of the samples subjected to additional sequencing with- bacterium tuberculosis, the causative agent of tubercu-
out exome capture (SA1000 and SA1025). Specifically, losis, a disease that is common in the Northern Cape
we find a significant, positive correlation (Spearman’s region of South Africa [39]. These reads map with equally
rho > 0.65; p-value < 2.2e-16) between the relative abun- high fidelity to the genomes of other Actinobacteria, sug-
dance estimates calculated with the two sequencing gesting that they may be homologs of ancient and highly
approaches for both samples (Figure 5), indicating conserved Actinobacteria sequences and are not necessar-
that analysis of exome-captured metagenomes produces ily representatives of the M. tuberculosis genome. Robust
saliva community structure and abundance estimates that detection of M. tuberculosis from saliva-derived exome
capture sequence data requires additional experimentation
and validation.
Table 4 KhoeSan saliva microbiome abundance by read
The predominant saliva microbiota differ in their rela-
threshold
tive abundance across the KhoeSan (Figure 6). To assess
No coverage1 75% Read coverage2
whether population structure based on saliva micro-
Genome 50% Identity 80% Identity 95% Identity biome diversity exists among the KhoeSan, we clustered
Neisseria subflava 0.062 0.077 0.13 individuals based on their phylum-, genus-, or species-
Rothia mucilaginosa 0.051 0.062 0.076 level relative abundances. We find only moderate sup-
Neisseria flavescens 0.031 0.039 0.063 port for the existence of discrete clusters among the
Streptococcus 0.038 0.046 0.063 KhoeSan, with a maximum average silhouette width of
parasanguinis 0.46 (genus-level clustering). Following [40], this sug-
Prevotella melaninogenica 0.047 0.054 0.059 gests that saliva microbiome diversity varies among the
Veillonella dispar 0.041 0.049 0.056
KhoeSan along a gradient. We subjected the microbiome
abundances of the KhoeSan samples to Principal Com-
Prevotella veroralis 0.027 0.031 0.029
ponents Analysis (PCA) to identify those taxonomic
Streptococcus salivarius 0.016 0.019 0.028 groups that most strongly differentiate the samples along
Granulicatella elegans 0.02 0.025 0.027 this gradient (i.e., maximum PCA loadings). At the
Fusobacterium 0.015 0.012 0.024 phylum-level, KhoeSan saliva samples are principally
periodonticum separated by their relative composition of Proteobac-
1
The genome-length-corrected relative abundance calculated using a 50% teria, Firmicutes and Bacteroidetes (Additional file 1:
identity fragment recruitment threshold.
2
The genome-length-corrected relative abundance calculated using fragment
Figure S9). The relative abundance of Neisseria, Strepto-
recruitment thresholds of 80% or 95% identity across at least 75% of the read. coccus, Prevotella, and Veillonella primarily differentiate
Kidd et al. BMC Genomics 2014, 15:262 Page 10 of 17
http://www.biomedcentral.com/1471-2164/15/262

Figure 5 Comparison of saliva microbiome frequencies from full genome and exome-capture sequencing. Estimates of the relative
abundance of saliva microbiota obtained via exome capture (x-axis) strongly correlate with those obtained from shotgun metagenomes
produced from the same sample (y-axis). The above dot plots illustrate this result for two KhoeSan individuals involved in our study: A) SA1000
and B) SA1025. Each dot represents a genome. A linear model representing the relationship between exome-capture and non-capture estimates
of relative abundance is shown in blue; the variance in the predictions from the model are shaded in grey. A Spearman correlation test indicates
that this relationship is very strong (rho > 0.65; p < 2.2e-16).

1.00 1.00 Khoesan North American

A B C Streptococcus
Prevotella
Rothia
Veillonella
Neisseria
Porphyromonas
0.75 0.75 Fusobacterium
Species Actinomyces
Genus Granulicatella
Relative Abundance

Relative Abundance

Neisseria Rothia mucilaginosa Campylobacter

Streptococcus Prevotella melaninogenica Propionibacterium
Prevotella Streptococcus parasanguinis Corynebacterium
Veillonella Neisseria subflava Gemella
0.50 0.50 Veillonella dispar Atopobium
Fusobacterium Selenomonas
Rothia Veillonella sp. 3_1_44
Streptococcus oralis Leptotrichia
Actinomyces Eubacterium
Porphyromonas Porphyromonas gingivalis Filifactor
Other Streptococcus mitis Haemophilus
Other Capnocytophaga
0.25 0.25 Treponema
Oribacterium
Dialister
Mobiluncus
Aggregatibacter
Peptostreptococcus
Kingella
Bacteroides
0.00 0.00 Catonella
SA054
SA008
SA052
SA1000
SA051
SA012
SA011
SA1002
SA1001
SA006
SA1006
SA1010
SA1011
SA035
SA1025
SA008
SA052
SA051
SA1006
SA012
SA011
SA1002
SA006
SA054
SA1000
SA1010
SA1001
SA1025
SA035
SA1011

1.00 0.75 0.50 0.25 0.00

Normalized Relative Abundance

Figure 6 Differences in taxon ranks between South African samples and human microbiome project. A, B) Oral microbiome structure
varies among the KhoeSan. Each of the above stacked bar plots illustrates the relative abundance (y-axis) of the most abundant oral microbiota at
the A) genus, and B) species levels for each of the 15 KhoeSan individuals (x-axis). Relative abundance was measured as the fraction of high-quality
reads that were recruited to a microbial genome of a particular taxonomic rank using conservative recruitment settings (Methods). Only the nine
most abundant groups for each taxonomic level are illustrated for visualization purposes, with the remaining taxa being grouped into the ‘Other’
category. C) KhoeSan (red) and healthy North American (blue) saliva microbiomes differ in their community structure. In this bar plot, the normalized
relative abundance, which is a taxon’s median relative abundance detected within a population divided by the maximum relative abundance detected
within a population, is shown for bacterial genera that are detected in either of the two populations. Genera are ordered by their median
relative abundance across the KhoeSan. Notable differences between the populations are those where the taxon is abundant in the KhoeSan
and effectively undetected in the North Americans, especially Rothia.
Kidd et al. BMC Genomics 2014, 15:262 Page 11 of 17
http://www.biomedcentral.com/1471-2164/15/262

samples at the genus-level (Figure 6A), while Rothia muci- remain highly sensitive to coverage and calling pipe-
laginosa, Neisseria subflava, Veillonella spp., and Strepto- lines thus making direct cross-study comparisons diffi-
coccus mitis are some of the most variable species cult. However, for common SNPs, we show that the
amongst the KhoeSan (Figure 6B). KhoeSan strongly differentiate from all other human
populations in structure analyses; the KhoeSan and
North American versus south African oral microbiomes Europeans fall at opposite ends of the 1st principal compo-
We then compared the diversity of the KhoeSan oral nent, while western and eastern Africans fall at intermedi-
microbiome to the diversity observed in a recent and ate points on this axis. Furthermore, we find substantial
extensive survey of healthy North Americans in the sub-structure among the South African and Namibian
Human Microbiome Project (N = 294) [41]. This prior KhoeSan, despite recent gene flow from Bantu-speaking
HMP work was conducted through analysis of small groups and Europeans into the ≠Khomani, !Kung and Tuu
subunit ribosomal RNA (i.e., 16S rRNA) gene sequences populations.
that were taxonomically annotated to the genus level.
We used these sequences to calculate genus-level, gen- Sample quality
ome length-normalized relative abundances for each Two of our samples had demonstrably lower mapping
North American microbiome. We used the taxonomy as- quality and coverage, SA006 and SA035. We consider
sociated with each genome in our fragment recruitment three possibilities for these characteristics. First, it is dif-
database to calculate genus-level, length-normalized rela- ficult to identify the proportion of human DNA versus
tive abundances for each KhoeSan microbiome. Compar- microbial or other non-human DNA in a saliva aliquot.
ing each population’s median relative abundance for each If these two samples had by chance a lower volume of
genus, we find that most taxa exist at similar abundance human DNA input for the exome capture reaction, then
levels in the two populations (Spearman’s rho = 0.91, there would be fewer opportunities for human DNA to
p-value < 2.2e-16). However, there are five genera that bind to the specific probes and the library would likely
are present in relatively high abundance (Bonferroni- result in a higher number of duplicate read pairs. SA006
corrected Wilcoxon rank sum test p < 0.01) in the KhoeSan and SA035 do display an increased duplicate rate (54%,
and effectively undetected among the North Americans 78% respectively), but SA008 also displays high duplicate
given the level of discovery in the HMP (Figure 6C): rate with minimal effect on mapping quality. Addition-
Rothia, Granulicatella, Haemophilus, Eubacterium, and ally, poorer mapping quality might be expected if the
Filifactor. Most notable among these is Rothia, which is microbial reads map to the human genome, perhaps due
the third most abundant genus in the KhoeSan and con- to near sequence identity between some portion of the
tains Rothia mucilaginosa, a known oral opportunistic human and microbial genomes [44].
pathogen that has been linked to systemic diseases [42,43]. A second possibility is that the total amount and qual-
ity of the human DNA input initially may have been suf-
Discussion ficient, but the presence of non-human substances such
Population history as residual tobacco or bacterial DNA may have acted as
The extremely high genetic diversity in the KhoeSan, esti- inhibitors, preventing normal binding to human probes.
mated from genome-wide SNP arrays and the “Bushman” Third, the DNA in these two samples may have been
genome, has renewed interest in understanding the popu- more degraded than the other six Pilot 1 samples. How-
lation history of southern Africans [2,15,26,27]. Com- ever, although we do observe an increase in substitutions
paratively few genomic sequences are publicly available at the start of the reads for SA006 and SA0035, we find
(6 individuals total) from the KhoeSan, and ascertain- no evidence of an ancient DNA degradation pattern in
ment bias on many of the standard SNP arrays may the post-capture sequence data. While the listed possi-
strongly skew estimates of genetic diversity in these pop- bilities appear unlikely, it is possible other patterns of
ulations. We have generated 15 exomes and 2 genomes degradation occur, in relatively young DNA extractions,
from the South African ≠Khomani San greatly expanding which have not been reported in the literature.
the number of genomic sequences available. Estimates of
genetic diversity from these South African individuals are Oral microbiome from exome sequencing
comparable to genetic diversity from the Yoruba from Approximately 5.1% of the sequence data generated did
Nigeria or Luhya from Kenya (Figure 3). While we do not not map to the human genome. Using a phylogenetically
find a higher number of private SNPs in the KhoeSan, this diverse set of reference genomes and a fragment recruit-
may be biased due to endogamy among the ≠Khomani San ment approach, we identified those unmapped reads
and differences in coverage or SNP calling/imputation that are homologs of regions in non-human genomes.
pipelines between 1000 Genomes and our procedure We find that most of the reads map to genomes of well-
(Figure 1). Heterozygosity and singleton identification described commensal microorganisms of the human
Kidd et al. BMC Genomics 2014, 15:262 Page 12 of 17
http://www.biomedcentral.com/1471-2164/15/262

mouth, suggesting that this sequencing platform produces the effect of lifestyle on microbiome composition as
relevant information about the human oral microbiome. most studies focus on individuals living contemporary
We also find that analysis of exome-capture metagenomes Western lifestyles. Similar to studies conducted in Western
produces microbiome diversity estimates consistent with populations [48,49], we find that the KhoeSan salivary
those obtained from non-exome-capture metagenomes, microbiome is dominated by a small number of taxa, with
indicating that this platform can be used to reliably quan- the Firmicutes or Proteobacteria predominating, and ex-
tify microbiome diversity and abundance. We note that hibits high diversity within and between individuals. These
other capture technologies or probe designs may result in observations suggest that the general structure of the
fewer off-target reads, and a corresponding reduction in KhoeSan salivary microbiome is generally similar to that
the ability to analyze the microbiome [45,46]. Additionally, found in Western individuals.
different saliva collection kits or the use of pre-collection However, when evaluating differences in the relative
mouth washes may effect the yield of microbial-derived abundance of genera associated with the KhoeSan and a
sequences. population of healthy Americans, we identified several
The large fraction of non-human sequences that do abundant taxa in the KhoeSan that were at very low
not map to our reference genomes are likely low quality abundance or undetected among the Americans. These
and degraded sequences or are reads from organisms differences in microbiome structure may be due to dif-
that are outside of the bounds of the phylogenetic diver- ferences in (1) the evolutionary history of the popula-
sity sampled in our reference database, such as viral ge- tions, (2) demographics, or (3) host environment or
nomes. The size of this fraction may be exacerbated by lifestyle, including diet and access to health care. Given
the relatively conservative alignment thresholds applied that we find many known pathogens among the most
during our analysis. Our ability to detect oral commen- abundant members of the KhoeSan microbiome and that
sals indicates that this human exome sequencing plat- many of the differentially detected genera contain known
form provides the added benefit of being able to assay oral pathogens (e.g., Rothia, Granulicatella, Filifactor),
biogeographic patters of oral microbiome diversity. Given we speculate that the relatively limited access to dental
that many of the non-human reads can be mapped with care, antibiotics and/or absence of water fluoridation
high stringency to genomes of known pathogens, we among the KhoeSan is driving most of the observed dif-
hypothesize that this sequencing platform may be useful ferences between populations. However, the biology of
as a diagnostic tool for the detection of disease and that several of the differentially abundant genera is not well
the data obtained may be used for inferring cryptic phe- understood, especially in the context of the commensal
notypes of the sampled individuals (e.g., periodontitis oral microbiome (e.g., Mobiluncus), or is principally lim-
status). Future studies that focus on the sensitivity and ited to the pathogenic members of the genus; such gen-
specificity of pathogen detection will be required to test era may contain species that played an important role in
this hypothesis. the coevolution between the KhoeSan and their salivary
As a cautionary note, one genome that recruits a microbiome. This may include pathogenic organisms,
substantial number of reads (9.4% of total reads) is such as Aggregatibacter actinomycetemcomitans, the causa-
Beggiatoa sp. PS. Beggiatoa have been found in sulphur tive agent of adolescent periodontal disease, which is com-
springs, sewage contaminated water, and hydrothermal mon in those of African descent [50] and a member of a
vents [47]; to date, no one has described the presence relatively abundant genus in the KhoeSan. Further study of
of Beggiatoa in the human mouth. We found that the the microbiomes associated with the KhoeSan and other
Beggiatoa-recruited reads map to short, unassembled con- diverse human populations (e.g., [51]), the microbiomic
tigs that exhibit significant similarity to clone libraries of differences between these populations (e.g., [52,53]), espe-
the human genome. Thus, we suspect that our detection cially across a variety of host physiological conditions, and
of Beggiatoa is the result of low quality human reads that the biology of commensal microbiota that are under-
fail to align to the human genome reference sequence but represented in Western populations is needed to compre-
do align to regions of the Beggiatoa genome. This observa- hensively differentiate the sources of variations observed
tion highlights the importance of considering the effect of between populations and to understand the coevolution
human genome contamination when using fragment re- between humans and their microbiome.
cruitment to study the human microbiome.
Conclusions
KhoeSan microbiome diversity We have demonstrated the ability to obtain high quality
Understanding KhoeSan microbiome diversity and struc- exome sequence data from saliva-derived human DNA.
ture provides insight into the co-evolution of the human We show that even samples with low human DNA pres-
microbiome, given the ancient divergence of KhoeSan ence can be successfully captured using exome in-solution
from other African populations. It additionally clarifies target probes. Additionally, after examining some of the
Kidd et al. BMC Genomics 2014, 15:262 Page 13 of 17
http://www.biomedcentral.com/1471-2164/15/262

most diverse human loci, we find that exon-capture is able 500 ml SureSelect wash buffer #1 for 15 min. at room
to enrich and facilitate high-resolution analysis of highly temperature, and three times with 500 ml SureSelect wash
polymorphic HLA and KIR genes from DNA extracted buffer #2 for 10 min at 65°C. DNA was eluted with 50 ml
from human saliva. We also demonstrated that exon- SureSelect elution buffer for 10 min at room temperature
captured DNA sequencing of saliva reveals insight into the and neutralized with 50 ml of SureSelect neutralization
structure and diversity of the oral microbiome. buffer. The captured product was purified with SPRI beads
and amplified by PCR. The quality and concentration of
Methods the sequencing libraries was verified by the Bioanalyzer
Samples High Sensitivity DNA kit (Agilent). Indexed samples were
Sampling of the ≠Khomani KhoeSan in Upington, South pooled in an equimolar ratio and sequenced on the Illu-
Africa and neighboring villages occurred in 2006. Insti- mina HiSeq2000 according to standard protocols. A similar
tution Review Board (IRB) approval was obtained from procedure was followed for the Pilot 2 samples with the
Stanford University. Individuals who were still living in SureSelect 44 Mb All-Exon capture probe library.
2011 were re-consented under a modified protocol (IRB
approved from Stanford University and Stellenbosch Read mapping and SNP calling
University, South Africa). ≠Khomani N|u-speaking in- Illumina sequencing reads were mapped to the human
dividuals, local community leaders, traditional leaders, genome reference sequence (GRCh37) following a stand-
non-profit organizations and a legal counselor were all ard pipeline informed by the best-practices as described
consulted regarding the aims of this research, prior to col- by the 1000 Genomes project [24,54] (Figure 1). Pilot 1
lection of DNA. All individuals consented orally to partici- reads were trimmed to be 75 bp in length; Pilot 2 reads
pation, with a second, local native speaker witnessing and were 101 bp in length. Reads were mapped and paired
were re-consented with written consent. DNA via saliva using bwa version 0.6.2 [55]. Unmapped reads were identi-
(Oragene® kits) and ethnographic information regarding fied at this stage and processed via the metagenomic pipe-
self-identified ancestry (N|u, Nama, or ‘Coloured’), lan- line. Duplicate read pairs were identified using Picard
guage and parental place of birth were collected for all (http://picard.sourceforge.net/). Base qualities were empir-
participants. ically recalibrated and indel realignment was performed
jointly across all samples using the Genome Analysis Tool
Exome capture Kit (GATK) v1.6 [25]. BAM files containing only uniquely
Library preparation and exome enrichment were per- mapped reads with duplicates removed were analyzed by
formed as described in the Agilent SureSelectXT Target the program SAMStat [56]. Fraction of reads on target was
Enrichment System for Illumina Paired-End Sequencing determined using snpEff.
Library (Version 1.1.1, January 2011). First, purified DNA Sequencing reads from the samples described in
from saliva samples was concentrated to a volume com- Schuster et al. [15] were obtained from the short read
patible with the library preparation protocol. 3 μg of con- archive and remapped to the GRCh37 assembly. The
centrated genomic DNA was fragmented to a median size exome capture data from Schuster et al. was single end
of 200 bp using the Covaris-S2 instrument with the fol- sequences obtained from the 454 pyrosequencing tech-
lowing settings: duty cycle 10%, intensity 5, cycles per nology. Reads were mapped using the bwasw option in
burst 200, and mode frequency sweeping for 180 s at 4°C. bwa version 0.5.9. Processing was performed as described
The fragmentation efficiency was evaluated on the Agilent above, with the exception of omitting the ‘homopolymer’
Bioanalyzer using DNA1000 chips. After end-repair and recalibration covariate and skipping the indel realign-
A-tailing, sequencing adapters were ligated onto the ment step which is not supported for 454 reads.
DNA fragments, followed by size-selection using SPRI
beads (Agencourt AmPure XP) and PCR amplification. Read substitution bias
The amplification product was purified with SPRI beads For Pilot 1, rates of nucleotide substitutions at each pos-
and the quantity and quality was assessed using the ition along the reads were determined by comparing the
Bioanalyzer DNA1000 chip. Five hundred nanograms mapped reads to their aligned human genome reference
of the adapter-ligated DNA library were concentrated sequence. We analyzed the first 1 million reads mapped
to 3.4 ml, mixed with hybridization buffer and DNA to chr1 for each sample, using only reads without any
blocker mix, and added to the SureSelect 50 Mb All- alignment indels or clipping (with a CIGAR string of
Exon capture probe library. The mixture was incubated for ‘75 M’ in the BAM file). For each read, we retrieved the
24 hours at 65°C in a thermal cycler. The hybridization corresponding aligned reference sequence using its
mixture was added to streptavidin-coated M-280 Dyna- mapped chromosomal position in the BAM file. The
beads (Invitrogen) and incubated for 30 min at room rates for each nucleotide substitution type were then cal-
temperature, with mixing. The beads were washed with culated as the ratio of the total number of observed
Kidd et al. BMC Genomics 2014, 15:262 Page 14 of 17
http://www.biomedcentral.com/1471-2164/15/262

changes of that type and the total number of reads, for to reference sequences matching individual HLA-A, -B
each position along the reads. Because reads mapping to and -C genotypes. HLA-A, -B and -C reference alleles
the reverse strand of the reference are reverse comple- were determined using bead-based sequence specific
mented in the BAM files, we performed the analysis sep- oligonucleotide probe hybridization and were described
arately for forward and reverse strand mapping reads. in [2]. The “-phase” function of SAMtools was used to
Reverse mapping reads therefore show the complemen- attribute phase for local alignments where possible due
tary substitution patterns at the 3′ end to the forward to the close proximity of exons and/or presence of
mapping reads at the 5′ end. highly heterozygous sequence (e.g. exons 2 and 3 of
HLA class I). Post-filtered read depth was used to de-
Population differentiation termine presence or absence of the variable-content
To perform principal component analysis we used SNP KIR genes. The KIR genes present and their alleles were
genotypes for individuals from several populations and determined for comparison of eight of the individuals
the EIGENSOFT software [28]. We used 11 KhoeSan in- using pyrosequencing methods as previously described
dividuals from our dataset (excluding SA011 and SA012 [64]. Individual SNP genotypes were confirmed visually
from Family 1 and SA052 and SA054 from Family 2), 4 from independent alignments of the filtered reads,
Namibian KhoeSan individuals from Schuster et al., 6 which were created using MIRA 3 [65,66]. All newly-
Namibian San (Ju|’hoansi) from the Human Genome discovered variants were confirmed for sequence and
Diversity Project (Martin et al, in prep. SRP036155) [57], phase using standard Sanger sequencing plus one or
and 13 individuals from each of the ASW, GBR, LWK, and more of pyrosequencing, DNA cloning or segregation
YRI populations from the 1000 Genomes Project [23]. in families.
Closely related individuals were excluded from all data-
sets. Sample ‘ABT’ was excluded from Schuster et al.’s Metagenomic pipeline
dataset since it clustered with the Bantu-speaking pop- We searched for genetic signatures of non-human or-
ulations in their analyses. Individuals selected from the ganisms by adopting the fragment recruitment approach
1000 Genomes Project all had more than 20x coverage outlined by Rusch et al. [31] (Figure 1). We first
for at least 70% of exome targets. To account for differ- trimmed reads and removed low-quality (i.e., reads that
ences in coverage and target regions, variants included meet any of the following conditions: mean quality score
in this analysis had genotype information for at least less than 25, length less than 50 bp, presence of ambigu-
95% of the individuals for a given analysis. VCFtools ous bases) and exact duplicate reads from the set of
[58] was used to count the number of shared and pri- those that did not map to the human genome using
vate SNPs between populations. prinseq [67]. We then compared the remaining high-
quality reads that did not map to the human genome to
HLA/KIR calling 1,285 genomes (Additional file 2: Table S4) obtained
To analyze the whole-exome data, all read-pairs that from the Joint Genomes Institute’s Integrated Microbial
mapped within hg19 coordinates, chr6:28702021-33392022, Genomes (IMG) database [44]. In the case of species
chr19:55228188-55383188 and chr19_gl000209_random, that have multiple genome-sequenced individuals, we
were extracted using SAMtools 0.1.18 [59] and split randomly selected a single individual genome to repre-
into separate fastq files for each individual. Read-pairs sent the species group. Each read was aligned to each
having more than five bases of quality score ≤3 were genome using blast (blastall -p blastn -z 16300000000 -e
removed (FASTX Toolkit 0.0.13 [http://hannonlab. 0.01 -m 8) and the resulting alignment summary statistics
cshl.edu/fastx_toolkit/]). The analysis pipeline was de- were used to infer each read’s taxonomy [45]. We explored
signed to detect all known and any novel HLA class I several classification thresholds, including alignment
and KIR SNP variants. Using Bowtie (version 0.12.7) e-value, alignment percent identity, and the ratio between
[60] read-pairs were harvested by mapping with low- the alignment length and the read length (i.e., coverage).
stringency to a given HLA or KIR gene (positive filter). We adopted several levels of threshold stringency to recruit
To ensure specificity, pairs that mapped to any homolo- reads to genomes for the purposes of inferring taxonomic
gous gene or pseudogene were removed (negative filter). diversity. Our thresholds were similar to those used in
The remaining reads were then aligned to a final refer- Rusch et al. [31], with modifications to account for the
ence sequence and the SNP variants ascertained using short length of our sequences.
SAMtools/bcf. Data used to generate filters and reference In the lenient case (i.e., distant homology), a read was
sequences was obtained from the ImmunoPolymorphism recruited to a genome if the two sequences shared a
Database and a set of fully-sequenced KIR haplotypes local alignment having at least 50% sequence identity.
[61-63]. To accommodate the high divergence of Using these parameters we identified 5,060,454 un-
HLA exons 2 and 3, the final alignments were made mapped sequences (20.9% of total unmapped reads) that
Kidd et al. BMC Genomics 2014, 15:262 Page 15 of 17
http://www.biomedcentral.com/1471-2164/15/262

exhibit significant similarity to the collection of reference JX523654 (3DP1*014), JX523634 (2DL4*024),
genomes. In the stringent case (i.e., recent homology), a JX523637 (2DL4*027), GQ890695 (3DL1*070),
read was recruited to a genome if the alignment covered at GQ890697 (3DL1*071), GU323347 (3DL2*052),
least 75% of the read and the sequences had at least 80% GU323348 (3DL2*053), GU323349 (3DL2*054),
identity. Applying these thresholds found that 16.8% of the JX523649 (3DL2*063)
reads (N = 4,064,899) can be recruited by non-human
genomes.
Additional files
To conduct species-level binning, we applied the
aforementioned coverage thresholds, but required that Additional file 1: Table S1. Summary Statistics for HGDP San Exomes.
the read and target genome share at least 95% identity. Table S2: HLA alleles. Table S3: KIR alleles. Figure S1: Pedigree structure
In all cases of classification, we applied an e-value for sequenced individuals. Figure S2. Cumulative coverage across the
Agilent target regions for Pilot 1 (A) and Pilot 2 (B) samples. Figure S3:
threshold of 10-3. We inferred a read’s taxonomy by Mapping quality for all reads. Figure S4: Assessment of base substitutions
transferring the taxonomic annotation of the genome se- from mapped reads. Figure S5: Venn Diagram comparing ≠ Khomani San
quence that produced the best alignment score while with Namibian exome samples. Figure S6: PCA with two relatives included.
Figure S7: Distribution of mapped reads along the N. subflava genome.
also passing the classification thresholds. If a read could Figure S8: The phylogenetic distribution of three non-human exome
not be placed into a species group based on the refer- capture sequences that map with high fidelity to Mycobacterium
ence genomes, it was discarded from the subsequent tuberculosis. Figure S9: The phylum-level structure of the oral microbiome
structure varies among the KhoeSan.
diversity analyses. The IMG taxonomic annotations as-
Additional file 2: Table S4. Species included in the microbial genome
sociated with the reference database genomes were used database.
to assign species-level binned reads into genera and
phyla.
Competing interests
To quantify genus-level saliva microbiome abundances MG is employed by Agilent Technologies. BMH and CRG hold stock in
among healthy Americans, we downloaded high-quality, 23andMe, Inc. CDB is on the scientific advisory board of Ancestry.com.
taxonomically annotated V35 16S rRNA Roche amplicon 23andMe, Inc. and Ancestry.com had no role in the study design, data
collection and analysis, decision to publish, or preparation of the manuscript.
sequences associated with 294 saliva samples from the
Human Micorbiome Project (HMP) Data Analysis and
Authors’ contributions
Coordination Center (http://www.hmpdacc.org/). A prior JMK, TJS, DB, PJN, ARM, MLC, MS and BMH analyzed data; CRG, EGH, and
study used the Ribosomal Database Project classifier (v2.2) BMH collected and processed DNA samples; NN, AA, MG, XG, QF, YL, and XL
with the default 032010 training set and taxonomy to an- generated genomic data; TJS and KSP performed metagenomic analysis; PP,
MWF, KSP, CDB and JDW contributed to study design; JMK, TJS and BMH
notate these sequences [49]. Genus-level taxonomic assign- wrote the manuscript with input from all authors. All authors read and
ments were extracted for each sequence having a bootstrap approved the final manuscript.
statistic greater than 80%.
Acknowledgements
Availability of supporting data We extend our gratitude to Blanca Herrera for assistance with genome
sequencing. We thank Joanna Mountain, Julie Granka, Marlo Möller, Cedric
VCF files are available at http://ecoevo.stonybrook.edu/ Werely for their help with sample collection. Finally, we express our appreciation
hennlab/data/. Raw read data can be downloaded from to the ≠ Khomani San community for participation in our research projects.
the short-read archive (SRP038015 for saliva derived B.M.H. and C.D.B. are supported by NIH grant R01HG003229. C.R.G. is supported
by the UCSF Dissertation Year Fellowship and NIH grants T32GM007175 and
exomes and genomes, and SRP036155 for HGDP San T32HG000044. T.J.S and K.S.P. are supported by NSF grant DMS-1069303, the San
exomes). SNP variants have been deposited in dbSNP Simeon Fund, and institutional funding from Gladstone Institutes. P.J.N, N.N-G
(SS 974432427-SS974514519) Novel KIR alleles have and P.P. were supported by NIH grant AI17892. J.D.W. was supported by NIH
grant R01HG400409. J.M.K was supported by NIH grant 1DP5OD009154. A.R.M.
been deposited in Genbank and assigned Immuno Poly- was supported by NIH training grant GM007790.
morphism Database nomenclature as follows:
Author details
1
Department of Genetics, Stanford University, Stanford, CA 94305, USA.
JX523651 (3DL3*057), GQ924778 (3DL3*037), 2
Departments of Human Genetics, and Computational Medicine and
GQ924779 (3DL3*038), GQ924781 (3DL3*040), Bioinformatics, University of Michigan, Ann Arbor, MI, USA. 3The J. David
HM235773 (3DL3*041), JX523631 (2DL2*012), Gladstone Institutes, University of California, San Francisco, San Francisco,
CA 94158, USA. 4Departments of Microbiology, and Statistics, Oregon State
JX523638 (2DL5B*00803), JX523639 (2DL5B*018), University, Corvallis, OR 97331, USA. 5Department of Ecology and Evolution,
JX523640 (2DS3*007), JX523642 (2DS5*012), Stony Brook University, Life Sciences Bldg, Room 640, Stony Brook, NY 11794,
HM358896 (2DS5*0502), JX523648 (2DP1*00103), USA. 6Department of Structural Biology, Stanford University, Stanford, CA
94305, USA. 7Program in Pharmaceutical Sciences and Pharmacogenomics,
JX523646 (2DP1*00202), JX523647 (2DP1*011), University of California, San Francisco, CA 94143, USA. 8Agilent Technologies,
JX523644 (2DP1*012), JX523645 (2DP1*013), Genomics Division, Cedar Creek, TX 78612, USA. 9Translational Medicine, BGI –
JX523643 (2DP1*014), JX523630 (2DL1*026 N), Shenzhen, Shenzhen, China. 10Stellenbosch University, Tygerberg, South Africa.
11
Department of Biological Sciences, Stanford University, Stanford, CA 94305,
GU323355 (2DL1*022), JX523652 (3DP1*011), USA. 12Institute for Human Genetics, and the Departments of Epidemiology and
JX523655 (3DP1*012), JX523653 (3DP1*013), Biostatistics, University of California, San Francisco, San Francisco, CA 94143, USA.
Kidd et al. BMC Genomics 2014, 15:262 Page 16 of 17
http://www.biomedcentral.com/1471-2164/15/262

Received: 18 March 2014 Accepted: 28 March 2014 19. Briggs AW, Stenzel U, Johnson PL, Green RE, Kelso J, Prüfer K, Meyer M,
Published: 4 April 2014 Krause J, Ronan MT, Lachmann M, Pääbo S: Patterns of damage in
genomic DNA sequences from a Neandertal. Proc Natl Acad Sci U S A
2007, 104(37):14616–14621.
References 20. Stoneking M, Krause J: Learning about human population history
1. Liu J, Morgan M, Hutchison K, Calhoun VD: A study of the influence of sex from ancient and modern genomes. Nat Rev Genet 2011,
on genome wide methylation. PLoS One 2010, 5(4):e10028. 12(9):603–614.
2. Henn BM, Gignoux CR, Jobin M, Granka JM, Macpherson JM, Kidd JM, 21. Ginolhac A, Rasmussen M, Gilbert MT, Willerslev E, Orlando L: mapDamage:
Rodríguez-Botigué L, Ramachandran S, Hon L, Brisbin A, Lin AA, testing for damage patterns in ancient DNA sequences. Bioinformatics
Underhill PA, Comas D, Kidd KK, Norman PJ, Parham P, Bustamante CD, 2011, 27(15):2153–2155.
Mountain JL, Feldman MW: Hunter-gatherer genomic diversity 22. International HapMap 3 Consortium, Altshuler DM, Gibbs RA, Peltonen L,
suggests a southern African origin for modern humans. Proc Natl Altshuler DM, Gibbs RA, Peltonen L, Dermitzakis E, Schaffner SF, Yu F,
Acad Sci U S A 2011, 108(13):5154–5162. Peltonen L, Dermitzakis E, Bonnen PE, Altshuler DM, Gibbs RA, de Bakker PI,
3. Kurek KC, Luks VL, Ayturk UM, Alomari AI, Fishman SJ, Spencer SA, Mulliken JB, Deloukas P, Gabriel SB, Gwilliam R, Hunt S, Inouye M, Jia X, Palotie A, Parkin M,
Bowen ME, Yamamoto GL, Kozakewich HP, Warman ML: Somatic mosaic Whittaker P, Yu F, Chang K, Hawes A, Lewis LR, Ren Y, et al: Integrating
activating mutations in PIK3CA cause CLOVES syndrome. Am J Hum Genet common and rare genetic variation in diverse human populations.
2012, 90(6):1108–1115. Nature 2010, 467(7311):52–58.
4. Deng X: SeqGene: a comprehensive software solution for mining exome- 23. 1000 Genomes Project Consortium, Abecasis GR, Auton A, Brooks LD,
and transcriptome- sequencing data. BMC Bioinforma 2011, 12:267. DePristo MA, Durbin RM, Handsaker RE, Kang HM, Marth GT, McVean GA:
5. Shearer AE, Hildebrand MS, Smith RJ: Solution-based targeted genomic An integrated map of genetic variation from 1,092 human genomes.
enrichment for precious DNA samples. BMC Biotechnol 2012, 12:20. Nature 2012, 491(7422):56–65.
6. Kitzman JO, Snyder MW, Ventura M, Lewis AP, Qiu R, Simmons LE, Gammill HS, 24. DePristo MA, Banks E, Poplin R, Garimella KV, Maguire JR, Hartl C, Philippakis AA,
Rubens CE, Santillan DA, Murray JC, Tabor HK, Bamshad MJ, Eichler EE, del Angel G, Rivas MA, Hanna M, McKenna A, Fennell TJ, Kernytsky AM,
Shendure J: Noninvasive whole-genome sequencing of a human Sivachenko AY, Cibulskis K, Gabriel SB, Altshuler D, Daly MJ: A framework for
fetus. Sci Transl Med 2012, 4(137):76. variation discovery and genotyping using next-generation DNA sequencing
7. Patel ZH, Kottyan LC, Lazaro S, Williams MS, Ledbetter DH, Tromp H, Rupert A, data. Nat Genet 2011, 43(5):491–498.
Kohram M, Wagner M, Husami A, Qian Y, Valencia CA, Zhang K, Hostetter MK, 25. McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A,
Harley JB, Kaufman KM: The struggle to find reliable results in exome Garimella K, Altshuler D, Gabriel S, Daly M, DePristo MA: The genome
sequencing data: filtering out Mendelian errors. Front Genet 2014, 5:16. analysis toolkit: a MapReduce framework for analyzing next-generation
8. Teer JK, Mullikin JC: Exome sequencing: the sweet spot before whole DNA sequencing data. Genome Res 2010, 20(9):1297–1303.
genomes. Hum Mol Genet 2010, 19(R2):R145–R151. 26. Schlebusch CM, Skoglund P, Sjödin P, Gattepaille LM, Hernandez D, Jay F,
9. Bamshad MJ, Ng SB, Bigham AW, Tabor HK, Emond MJ, Nickerson DA, Li S, De Jongh M, Singleton A, Blum MG, Soodyall H, Jakobsson M:
Shendure J: Exome sequencing as a tool for Mendelian disease gene Genomic variation in seven Khoe-San groups reveals adaptation and
discovery. Nat Rev Genet 2011, 12(11):745–755. complex African history. Science 2012, 338(6105):374–379.
10. Bustamante CD, Fledel-Alon A, Williamson S, Nielsen R, Hubisz MT, 27. Pickrell JK, Patterson N, Barbieri C, Berthold F, Gerlach L, Güldemann T,
Glanowski S, Tanenbaum DM, White TJ, Sninsky JJ, Hernandez RD, Kure B, Mpoloka SW, Nakagawa H, Naumann C, Lipson M, Loh PR, Lachance J,
Civello D, Adams MD, Cargill M, Clark AG: Natural selection on Mountain J, Bustamante CD, Berger B, Tishkoff SA, Henn BM, Stoneking M,
protein-coding genes in the human genome. Nature 2005, Reich D, Pakendorf B: The genetic prehistory of southern Africa. Nat Commun
437(7062):1153–1157. 2012, 3:1143.
11. Tennessen JA, Madeoy J, Akey JM: Signatures of positive selection 28. Patterson N, Price AL, Reich D: Population structure and eigenanalysis.
apparent in a small sample of human exomes. Genome Res 2010, PLoS Genet 2006, 2(12):e190.
20(10):1327–1334. 29. Parham P: MHC class I molecules and KIRs in human history, health and
12. Yi X, Liang Y, Huerta-Sanchez E, Jin X, Cuo ZX, Pool JE, Xu X, Jiang H, survival. Nat Rev Immunol 2005, 5(3):201–214.
Vinckenbosch N, Korneliussen TS, Zheng H, Liu T, He W, Li K, Luo R, 30. Parham P, Norman PJ, Abi-Rached L, Hilton HG, Guethlein LA: Review:
Nie X, Wu H, Zhao M, Cao H, Zou J, Shan Y, Li S, Yang Q, Asan, Ni P, immunogenetics of human placentation. Placenta 2012, 33(Suppl):S71–S80.
Tian G, Xu J, Liu X, Jiang T, Wu R, et al: Sequencing of 50 human exomes 31. Rusch DB, Halpern AL, Sutton G, Heidelberg KB, Williamson S, Yooseph S,
reveals adaptation to high altitude. Science 2010, 329(5987):75–78. Wu D, Eisen JA, Hoffman JM, Remington K, Beeson K, Tran B, Smith H,
13. Rylander-Rudqvist T, Håkansson N, Tybring G, Wolk A: Quality and quantity Baden-Tillson H, Stewart C, Thorpe J, Freeman J, Andrews-Pfannkoch C,
of saliva DNA obtained from the self-administrated oragene method–a Venter JE, Li K, Kravitz S, Heidelberg JF, Utterback T, Rogers YH, Falcón LI,
pilot study on the cohort of Swedish men. Cancer Epidemiol Biomarkers Souza V, Bonilla-Rosso G, Eguiarte LE, Karl DM, Sathyendranath S, et al:
Prev 2006, 15(9):1742–1745. The sorcerer II global ocean Sampling expedition: Northwest Atlantic
14. Hansen TV, Simonsen MK, Nielsen FC, Hundrup YA: Collection of blood, through eastern tropical Pacific. PLoS Biol 2007, 5(3):e77.
saliva, and buccal cell samples in a pilot study on the Danish nurse 32. Konstantinidis KT, Ramette A, Tiedje JM: The bacterial species definition in the
cohort: comparison of the response rate and quality of genomic DNA. genomic era. Philos Trans R Soc Lond B Biol Sci 2006, 361(1475):1929–1940.
Cancer Epidemiol Biomarkers Prev 2007, 16(10):2072–2076. 33. Yanagisawa M, Kuriyama T, Williams DW, Nakagawa K, Karasawa T:
15. Schuster SC, Miller W, Ratan A, Tomsho LP, Giardine B, Kasson LR, Harris RS, Proteinase activity of prevotella species associated with oral purulent
Petersen DC, Zhao F, Qi J, Alkan C, Kidd JM, Sun Y, Drautz DI, Bouffard P, infection. Curr Microbiol 2006, 52(5):375–378.
Muzny DM, Reid JG, Nazareth LV, Wang Q, Burhans R, Riemer C, Wittekindt NE, 34. Peng Z, Fives-Taylor P, Ruiz T, Zhou M, Sun B, Chen Q, Wu H: Identification
Moorjani P, Tindall EA, Danko CG, Teo WS, Buboltz AM, Zhang Z, Ma Q, of critical residues in Gap3 of Streptococcus parasanguinis involved in Fap1
Oosthuysen A, et al: Complete Khoisan and Bantu genomes from southern glycosylation, fimbrial formation and in vitro adhesion. BMC Microbiol
Africa. Nature 2010, 463(7283):943–947. 2008, 8:52.
16. Gronau I, Hubisz MJ, Gulko B, Danko CG, Siepel A: Bayesian inference of 35. Ohara-Nemoto Y, Kishi K, Satho M, Tajika S, Sasaki M, Namioka A, Kimura S:
ancient human demography from individual genome sequences. Infective endocarditis caused by Granulicatella elegans originating in the
Nat Genet 2011, 43(10):1031–1034. oral cavity. J Clin Microbiol 2005, 43(3):1405–1407.
17. Asan, Xu Y, Jiang H, Tyler-Smith C, Xue Y, Jiang T, Wang J, Wu M, Liu X, 36. Gibson FC 3rd, Hong C, Chou HH, Yumoto H, Chen J, Lien E, Wong J,
Tian G, Wang J, Wang J, Yang H, Zhang X: Comprehensive comparison Genco CA: Innate immune recognition of invasive bacteria accelerates
of three commercial human whole-exome capture platforms. Genome atherosclerosis in apolipoprotein E-deficient mice. Circulation 2004,
Biol 2011, 12(9):R95. 109(22):2801–2806.
18. Clark MJ, Chen R, Lam HY, Karczewski KJ, Chen R, Euskirchen G, Butte AJ, 37. Zeituni AE, Carrion J, Cutler CW: Porphyromonas gingivalis-dendritic cell
Snyder M: Performance comparison of exome DNA sequencing interactions: consequences for coronary artery disease. J Oral Microbiol
technologies. Nat Biotechnol 2011, 29(10):908–914. 2010, 2:5782.
Kidd et al. BMC Genomics 2014, 15:262 Page 17 of 17
http://www.biomedcentral.com/1471-2164/15/262

38. Ihara H, Miura T, Kato T, Ishihara K, Nakagawa T, Yamada S, Okuda K: 62. Pyo CW, Guethlein LA, Vu Q, Wang R, Abi-Rached L, Norman PJ, Marsh SG,
Detection of Campylobacter rectus in periodontitis sites by monoclonal Miller JS, Parham P, Geraghty DE: Different patterns of evolution in the
antibodies. J Periodontal Res 2003, 38(1):64–72. centromeric and telomeric regions of group A and B haplotypes of the
39. Chimusa ER, Zaitlen N, Daya M, Möller M, van Helden PD, Mulder NJ, human killer cell Ig-like receptor locus. PLoS One 2010, 5(12):e15115.
Price AL, Hoal EG: Genome-wide association study of ancestry- 63. Robinson J, Mistry K, McWilliam H, Lopez R, Marsh SG: IPD–the
specific TB risk in the South African Coloured population. Hum Mol immuno polymorphism database. Nucleic Acids Res 2010,
Genet 2013, 23:796. doi: 10.1093/hmg/ddt462. 38(Database issue):D863–D869.
40. Koren O, Knights D, Gonzalez A, Waldron L, Segata N, Knight R, Huttenhower C, 64. Norman PJ, Hollenbach JA, Nemat-Gorgani N, Guethlein LA, Hilton HG,
Ley RE: A guide to enterotypes across the human body: meta-analysis of Pando MJ, Koram KA, Riley EM, Abi-Rached L, Parham P: Co-evolution of
microbial community structures in human microbiome datasets. human leukocyte antigen (HLA) class I ligands with killer-cell
PLoS Comput Biol 2013, 9(1):e1002863. immunoglobulin-like receptors (KIR) in a genetically diverse population
41. Human Microbiome Project Consortium: A framework for human of sub-Saharan Africans. PLoS Genet 2013, 9(10):e1003938.
microbiome research. Nature 2012, 486(7402):215–21. 65. Chevreux B, Pfisterer T, Drescher B, Driesel AJ, Müller WE, Wetter T, Suhai S:
42. Hodzic E, Snyder S: A case of peritonitis due to Rothia mucilaginosa. Using the miraEST assembler for reliable and automated mRNA
Perit Dial Int 2010, 30(3):379–380. transcript assembly and SNP detection in sequenced ESTs. Genome Res
43. Pinsky RL, Piscitelli V, Patterson JE: Endocarditis caused by relatively 2004, 14(6):1147–1159.
penicillin-resistant Stomatococcus mucilaginosus. J Clin Microbiol 1989, 66. Staden R, Beal KF, Bonfield JK: The Staden package, 1998. Methods Mol Biol
27(1):215–216. 2000, 132:115–130.
44. Liu Y, Li J: Short regions of sequence identity between the genomes of 67. Schmieder R, Edwards R: Quality control and preprocessing of
bacteria and human. Curr Microbiol 2011, 62(3):770–776. metagenomic datasets. Bioinformatics 2011, 27(6):863–864.
45. Bodi K, Perera AG, Adams PS, Bintzler D, Dewar K, Grove DS, Kieleczawa J,
Lyons RH, Neubert TA, Noll AC, Singh S, Steen R, Zianni M: Comparison of doi:10.1186/1471-2164-15-262
commercially available target enrichment methods for next-generation Cite this article as: Kidd et al.: Exome capture from saliva produces high
sequencing. J Biomol Tech 2013, 24(2):73–86. quality genomic and metagenomic data. BMC Genomics 2014 15:262.
46. Parla JS, Iossifov I, Grabill I, Spector MS, Kramer M, McCombie WR: A
comparative analysis of exome capture. Genome Biol 2011, 12(9):R97.
47. Larkin JM, Strohl WR: Beggiatoa, Thiothrix, and Thioploca. Annu Rev
Microbiol 1983, 37:341–367.
48. Lazarevic V, Whiteson K, Hernandez D, François P, Schrenzel J: Study of
inter- and intra-individual variations in the salivary microbiota.
BMC Genomics 2010, 11:523.
49. Consortium HMP: Structure, function and diversity of the healthy human
microbiome. Nature 2012, 486(7402):207–214.
50. Henderson B, Ward JM, Ready D: Aggregatibacter (Actinobacillus)
actinomycetemcomitans: a triple A* periodontopathogen? Periodontol 2000
2010, 54(1):78–105.
51. Nasidze I, Li J, Schroeder R, Creasey JL, Li M, Stoneking M: High diversity of
the saliva microbiome in Batwa Pygmies. PLoS One 2011, 6(8):e23352.
52. Nasidze I, Li J, Quinque D, Tang K, Stoneking M: Global diversity in the
human salivary microbiome. Genome Res 2009, 19(4):636–643.
53. Yatsunenko T, Rey FE, Manary MJ, Trehan I, Dominguez-Bello MG, Contreras M,
Magris M, Hidalgo G, Baldassano RN, Anokhin AP, Heath AC, Warner B, Reeder J,
Kuczynski J, Caporaso JG, Lozupone CA, Lauber C, Clemente JC, Knights D,
Knight R, Gordon JI: Human gut microbiome viewed across age and
geography. Nature 2012, 486(7402):222–227.
54. Consortium GP: A map of human genome variation from population-
scale sequencing. Nature 2010, 467(7319):1061–1073.
55. Li H, Durbin R: Fast and accurate short read alignment with Burrows-
Wheeler transform. Bioinformatics 2009, 25(14):1754–1760.
56. Lassmann T, Hayashizaki Y, Daub CO: SAMStat: monitoring biases in next
generation sequencing data. Bioinformatics 2011, 27(1):130–131.
57. Cann HM, de Toma C, Cazes L, Legrand MF, Morel V, Piouffre L, Bodmer J,
Bodmer WF, Bonne-Tamir B, Cambon-Thomsen A, Chen Z, Chu J, Carcassi C,
Contu L, Du R, Excoffier L, Ferrara GB, Friedlaender JS, Groot H, Gurwitz D,
Jenkins T, Herrera RJ, Huang X, Kidd J, Kidd KK, Langaney A, Lin AA,
Mehdi SQ, Parham P, Piazza A: A human genome diversity cell line
panel. Science 2002, 296(5566):261–262.
58. Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, Submit your next manuscript to BioMed Central
Handsaker RE, Lunter G, Marth GT, Sherry ST, McVean G, Durbin R, 1000 and take full advantage of:
Genomes Project Analysis Group: The variant call format and VCFtools.
Bioinformatics 2011, 27(15):2156–2158.
• Convenient online submission
59. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G,
Abecasis G, Durbin R, 1000 Genome Project Data Processing • Thorough peer review
Subgroup: The sequence alignment/map format and SAMtools. • No space constraints or color figure charges
Bioinformatics 2009, 25(16):2078–2079.
• Immediate publication on acceptance
60. Langmead B, Trapnell C, Pop M, Salzberg SL: Ultrafast and memory-
efficient alignment of short DNA sequences to the human genome. • Inclusion in PubMed, CAS, Scopus and Google Scholar
Genome Biol 2009, 10(3):R25. • Research which is freely available for redistribution
61. Wilson MJ, Torkar M, Haude A, Milne S, Jones T, Sheer D, Beck S, Trowsdale J:
Plasticity in the organization and sequences of human KIR/ILT gene families.
Proc Natl Acad Sci U S A 2000, 97(9):4778–4783. Submit your manuscript at
www.biomedcentral.com/submit