Unmappable Reads from the Human Microbiome Can Be Assembled and Annotated

We sequenced a total of 37 billion molecules from the 1,351 samples of cfDNA, of which 95% of reads passed quality control. Of these, an average of 0.45% did not align to the reference human genome (GRCh38) (Fig. 1A, Left), in line with our expectations of the nonhuman DNA sources in the body. Of these putatively nonhuman reads, only approximately 1% could be identified in a curated microbiome database of almost 8,000 species of known bacteria, viruses, fungi, and eukaryotic pathogens (Fig. 1B, Left). This miniscule fraction of reads encompasses the known microbiome. Less than 1,800 known species (800 known genera) are observed across all samples. The rarefaction curve of species prevalence quickly plateaus, and the species abundance distribution has only a slight positive skew (SI Appendix, Fig. S1). These both indicate that the number of known species we measure has saturated (20), and deeper or broader sequencing of cfDNA from humans is unlikely to substantially increase the richness of known species.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

(A and B) Violin plots showing the distributions of the percentage of reads associated with various stages of the assembly pipeline. From left: (A) reads that did not align to the reference human sequence; reads that were not removed during postassembly cleanups of additional human-like or low-complexity sequences; and (B) nonhuman reads that aligned to a curated microbiome database; nonhuman reads that aligned to the cohort assemblies. The white dot is the median value, and the thicker black bar is the interquartile range. (C) The total size of novel, divergent, and known contigs in the cohort assemblies. (D) Number of reads aligning to novel, divergent, and known contigs in the cohort assemblies. (E) Schematic of the iterative assembly process and novel contig selection. Between each step of assembly, further human and low-complexity contigs are blacklisted, with reads pooled from multiple samples/patients used in the next assembly. Multiple filters are performed on the long (>1 kbp) contigs to select for sequences that have not been previously observed and are unlikely to be contaminants.

We performed de novo assembly on the remaining nonhuman reads to uncover new species in the dark matter of cfDNA. The construction of assemblies used an iterative approach (Fig. 1E, Top). Nonhuman reads were assembled on a per-sample basis, and reads that aligned to low-complexity or human-derived contigs were removed. This process of assembly and cleaning was repeated for remaining reads, pooled first by patient and then by cohort, and resulted in a total assembly of 40 Mbp. Over 25 megabases of low complexity or residual human-derived contigs were removed (15% of reads) (Fig. 1A, Right), many of which were identified as human microsatellites or primate BAC/FOSMID clones (Dataset S1). The iterative assembly process captures more reads in each stage as the number of reads pooled together increases (SI Appendix, Fig. S2). The cohort assemblies constructed 7,190 contigs larger than 1 kbp and 131 larger than 10 kbp. Compared with the proportion of reads that map to known organisms, an order-of-magnitude greater fraction of reads map to the cohort assemblies (Fig. 1B, Right).

To select for “novel” contigs likely to originate from uncharacterized genomes, a series of filtering steps were applied (Fig. 1E, Bottom). The first two filters enrich for contigs that have a high gene content (i.e., predicted genes span at least 60% of bases) and low homology at both the nucleotide and protein level to any previously known sequence (SI Appendix, Figs. S3 and S4; i.e., BLAST alignments span less than 20% of bases and an average gene identity of less than 60%). Application of these filters results in the selection of 4,354 contigs over 1 kbp. Two additional steps removed contigs that have homologies in expanded microbiome sequence databases. First, we aligned all assemblies from phases II and III of the HMP and blacklisted the 47 contigs that had matches. As many of the HMP assemblies are of previously known organisms or had been deposited in the National Center for Biotechnology Information (NCBI) nucleotide (nt) database, the BLAST coverage filter encompasses most of the HMP-filtered contigs (SI Appendix, Fig. S5). The second filter used the web service Onecodex (21) to remove contigs with homologies in their database.

To control for potential contaminants from the extraction columns, we prepared six sequencing libraries using the same protocol used in the plasma-based cfDNA samples, but instead of plasma, we used either water or DNA extracted from a human cell line. Filtered reads were aligned to the assembled contigs, resulting in blacklisting an extra 114 contigs due to their presence in control samples. An additional step to check for contaminants was performed using the nonhost reads from 300 cfDNA samples obtained from nonhuman primate plasma. No contigs were observed at a level above 1 read per kilobase in more than 158 samples, with 75% of contigs observed at this level in less than 27 samples. The highly variable and nonubiquitous expression of these contigs in primate samples indicates that these are not common contaminants from the laboratory or kits.

After all of the filters, a total of 4,187 contigs remain, which can be further reduced to 3,761 novel candidates after merging contigs with significant overlaps. For later comparisons, a further 773 contigs are classified as “known” (>80% BLAST coverage and >1 kbp) and a further 598 as “divergent” (>1 kbp and neither known nor novel). The majority of assembled bases are novel, with 11 Mbp assembled compared with 3.5 Mbp of known or divergent contigs (Fig. 1C). The number of reads aligning is similar for all classes of contig (Fig. 1D). The lower coverage of novel contigs may have kept them hidden from previous analyses, and only by pooling many samples together were these sequences found (22). Biases in read coverage across the bacterial contigs may also be used to estimate their growth dynamics (23).

Novel Contigs Are Not Artifacts or Contaminants

As validation that the novel contigs are not merely contaminants or assembly artifacts, we performed a series of additional tests of their existence. The fragmented nature of cfDNA prevents direct PCR amplification of large regions of the contigs, so instead we (i) screened external datasets and samples for novel sequences at the read level, (ii) used bioinformatic approaches to assess assembly quality, and (iii) independently measured the presence of novel contig sequences by PCR in additional samples.

Plasma aliquots from a lung transplant recipient in our cohort were independently prepared and sequenced in another laboratory using an alternative library preparation protocol (24). We downloaded the reads from the NCBI sequence read archive and processed them identically in the pipeline followed by aligning them to the database of novel contigs. The two sets of different preparations exhibit concordant rankings of the high-abundance contigs (SI Appendix, Fig. S6). We also downloaded data generated by another group (25) from pregnancy samples that were completely independently collected, extracted, and sequenced and compared the reads to our pregnancy cohort (SI Appendix, Fig. S7). Although these samples had only slightly more than a tenth as many nonhuman reads as our samples, there were still reads which mapped to the novel contigs discovered in our experiments. These external datasets provide additional evidence that the novel contigs are not local laboratory contaminants and that their presence is robust to different sample preparation protocols.

The quality of assembly was assessed using information from the iterative assembly approach. Alignments between contigs from different levels created a graph structure, allowing us to investigate the additive structure derived from pooling more samples together. Under ideal circumstances, contigs assembled at lower stages should be wholly contained at higher stages and they should not end up orphaned (i.e., having no homology to a later stage’s contig). Only 3,301 of the 25,765 sample contigs (over 500 bp, nonhuman, not low-complexity) are orphaned, compared with the patient contigs. Over 3,000 of these are from single-end read samples and are ignored by the assembler at later stages, as the longer paired-end reads tend to be more informative. The patient-orphaned contigs mostly had high BLAST coverages (10,768 of 15,321 had coverage >20%) and encompassed the following genera most frequently: Lactobacillus (5,975), Streptococcus (1,449), Saccharomyces (1,315), Leuconostoc (1,051), and Lactococcus (610). Likewise to the sample-to-patient case, the vast majority of the remaining contigs (>85%) derive from patients who only have single-end sequencing. If we look in the other direction, i.e., at the lower contigs contained within the novel contigs, the majority of them derive from a single patient’s assembly (with some reads recruited from others), and only 387 were not assembled at any lower stages. The increased complexity of assembling due to more pooling does not appear to reduce the power to assemble what is present or to produce strong evidence of chimeric sequences. The consistency of the novel microbiome was further analyzed by calculating the similarity between longitudinal samples from the same patients. Times within approximately 3 mo tend to be highly correlated before the observed contigs start to diverge.

An additional check of the pipeline with a synthetic control dataset was performed to test if relevant sequences are prematurely removed. This was comprised of 8,068 curated genomes from the NCBI of viruses, bacteria, and fungi from which high-quality synthetic 100-bp reads at 25-bp intervals were created. Only 28 genomes (Dataset S2) had more than 50% of their reads removed in the pipeline steps leading up to assembly. These include various plasmids and enterobacteria phages (expected to be filtered in UniVec), simian virus 40, a few candidate bacterial species, and two variants of the human papillomavirus. Less than 10% of genomes had more than 10% of their reads removed, and those that were removed are sequences known to be in the cleaning or subtraction databases, indicating that the subtraction steps are unlikely to remove relevant nonhuman sequences.

We performed PCR testing to verify the presence of novel contigs in 16 previously unsequenced samples derived from the plasma of eight heart transplant recipients. As the fragmented nature of cfDNA prevents the amplification of large fragments, we designed primers against numerous short (60–100 bp) templates. Sequencing libraries known to express the templates targeted were used as positive controls (fifth column in SI Appendix, Fig. S8). Eight primers for a novel TTV (LT_node_191) and nine primers from unrooted contigs (HT220, HT552, HT76, PR78) were tested (Dataset S3). One of the eight patients had evidence at both time points for four templates from the TTVs, and three additional patients had evidence for one TTV-derived template (Fig. 2). None of the unrooted templates were seen in these samples. This was not surprising, as the median of the maximum number of reads of the unrooted contigs in the samples is 21; in contrast, the TTVs have a median value of 180 reads observed. Only through sequencing many foreign molecules from hundreds of samples were the unrooted contigs (and other novel ones) able to be found.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Gel showing the presence of three segments of LT_node_191 (a novel torque teno virus) in the cfDNA deriving from the heart transplant recipients (eight patients, two samples per patient). Demultiplexed gels and no-template controls are in SI Appendix, Fig. S8.

Provisional Taxonomic Assignment of Novel Contigs

Among the 3,761 novel contigs we found 14 with ribosomal proteins, none of which resembled the 16S ribosomal subunit. This, coupled with the fact that many of the sequences are viral-like, led us to conclude that a taxonomic classification based on 16S rRNA sequences is unfeasible. Thus, the taxonomic position was estimated based on choosing the lowest common ancestor of the majority of genes with homology, based on the gene’s taxonomic assignment from alignment. For known contigs, this method agrees with the best nucleotide alignment of a contig at the species level over 85% of the time (94% accurate at genus level), indicating the suitability of this approach for taxonomic classification. Contigs with both bacterial and phage genes had their taxonomic placement refined by ignoring the phage genes if that made the assignment more specific. To visualize the taxonomic diversity, a “solar system” was constructed for all novel contigs (Fig. 3). Individual dots are contigs, with the rings representing the taxonomic levels; their radial position within the ring is proportional to the average gene identity (inner is high, outer is low), and the angular position is dependent on the superkingdom or phylum to which the contig was assigned. The 844 yellow contigs placed on the outermost ring cannot confidently be placed in any superkingdom based on gene homology because for almost all of them the predicted genes have no known homology. Over two thirds of sequences appear to be bacterial; the majority of these are most similar to proteobacteria. Most longer contigs (>5 kb) are bacterial or prophage-like and can only be placed at the phylum or superkingdom levels.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Solar system plot of all novel contigs (>1 kbp). Each ring represents a different taxonomic level, with the intraring radius representing the average gene homology. Points are randomly assigned within sectors based on their superkingdom or phylum, as represented by the gray dotted spokes. The yellow contigs outside the last ring could not be assigned to any superkingdom; the radial scattering is to help illustrate the density of contigs. Colors correspond to superkingdom (counterclockwise from 3 o’clock: bacteria, viruses, eukaryota, archaea), saturation to the phylum.

Many novel sequences cannot be classified based on their gene-level homologies. There is a clear trend in the reduction of both the proportion of genes identified and the degree of similarity to known genes as we compare the known, divergent, and novel contigs (Fig. 4 and SI Appendix, Fig. S9 A–C). This lack of homology is the reason there are 840 unrooted contigs unable to be placed into any existing superkingdom (the yellow points in the outer solar system of Fig. 3). The majority of these (n = 739) have four or fewer predicted genes. Only four of the unrooted contigs have any identified genes, and most of those that are identified are of unknown function (full table in Dataset S4).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

Examples of contigs and their gene assignments for known, divergent, and novel contigs. Width of the box is equivalent to 60 kbp. Genes are colored by superkingdom (salmon, bacteria; lime green, viruses; dodger blue, eukaryota; fuchsia, archaea; black, no homology; gray, no gene); saturation is proportional to the gene identity.

As alignment-based methods failed to identify any genes from the remaining unrooted contigs, an estimated superkingdom assignment was made based on the codon use of genes (26). Comparing all genes with and without homology shows that these classes possess significantly different codon uses (P < 1e-3 after Bonferroni correction, Mann–Whitney U test; SI Appendix, Fig. S10), which is not merely an artifact of differential GC bias (SI Appendix, Fig. S11). An adaptive boosting classifier trained on known genes allowed us to bin all except 16 of the unrooted contigs into single or multiple superkingdoms (SI Appendix, Fig. S12). The vast majority (n = 640) are bacteria, and, similarly to the rooted contigs, many others (n = 87) have viral and bacterial genes and are treated as phage or prophages. There are hints that horizontal gene transfer across superkingdoms may be prevalent in biological dark matter, with 106 unrooted contigs having diverse sets of genes. Despite this substantial taxonomic distance from known species, we were able to consistently assign the novel contigs to taxonomic locations on the tree of life.

Human Microbiome Has Been Substantially Undersampled in Previous Studies

The novel sequences indicate that the microbiome’s species richness is much higher than is estimated from known species. If we count all of the taxa the contigs are assigned to, the number of taxa of known and divergent contigs is ∼20% of the total number of species observed in the curated microbiome database. Conservatively counting contigs assigned to higher levels only once (e.g., only count “proteobacteria” once), we observe over 1,000 novel taxa. If the discovery rate is similar to known and divergent contigs, this gives a range of 1,000–5,000 new species (i.e., a 50–250% increase over known species).

The distributions of average gene identity within each rank are clearly distinct for novel, divergent, and known contigs (SI Appendix, Fig. S13), which is reflected in their solar system plots (SI Appendix, Figs. S14 and S15). Although the divergent and known contigs are fewer in number, they are much more precisely placed, both in their taxonomic level and their average gene identity (reflected in the points being placed closer to the rank rings in the solar system). As more sequences are identified or classified, we expect many of the novel contigs in the outer orbitals to be attracted toward the center. The novel contigs are not restricted to individual samples or patients and most appear ubiquitously across the cohort populations. Each novel contig is observed in a median of 51 patients, and the median number of novel contigs observed in each patient is 924 (also see histograms in SI Appendix, Figs. S16 and S17). Even though elements of the human microbiome have been well-studied by deep sequencing from particular body sites, it is evident that other niches accessed by plasma contribute substantial novel diversity, which we report here for the first time.

Human Microbiome Contains an Unexpected Diversity of Novel Phages and Viruses

From the 2,917 placed novel contigs, 276 (9%) correspond to novel viral sequences, which are predominantly either phages or torque teno viruses (TTVs). Distinguishing between a phage and its bacterial host is difficult with short sequences, as they both are prone to incorporate each other’s genes. Indeed, of the 523 contigs containing phage genes, 333 also have bacterial genes. Nonetheless, identifying these is important, as the contigs with the most predicted genes are all phage or prophage candidates. Half of the genes have no homology for the top 15 such contigs, with some (HT_node_2, HT_node_16, cluster_12) having over three quarters of genes without matches. Of the identified genes, their closest homologies are commonly hypothetical proteins. These two facts conspire to make functional annotation unfeasible at this time. The nonhypothetical genes present on the phage contigs include DNA primases, phage-tail proteins, terminases, capsid proteins, and DNA polymerases I and III. Many of their bacterial genes are from species within the proteobacteria phylum, with some from actinobacteria (HT_node_11, LT_node_6) or containing genes from a flavobacterium (bacteroides phylum) phage. Additional sequences were searched for in the recently published global virome from the Joint Genome Institute (JGI) (4). Only 21 of the unmerged contigs had matches over 50% of their lengths, and these derived from just six of JGI’s scaffolds. Three of these are from Coloradan soil samples, one each from a wastewater bioreactor, a freshwater lake, and the upper troposphere. Our contigs matching to these are mostly best assigned as a cyanobacteria phage or unclassified Siphoviridae. The phage-like contigs were highly prevalent across patients and exhibited only minor clustering by cohort (SI Appendix, Fig. S18). These phages are likely to be associated with bacteria in the background microbial flora. Given that these contigs are observed ubiquitously in our data (77–287 samples) and seen in the JGI data only in a few exotic environments, we hypothesize that they are bona fide members of the human virome, and their presence in the JGI data is due to coincident discovery of related environmental species.

In contrast, the novel TTVs show strong clustering among cohorts (SI Appendix, Fig. S19) on the basis of immune system status. TTVs encompass anellovirus families and are known to be enriched in immunocompromised patients (17). Their detection here represents another validation of the taxonomic assignment based on gene homology. A phylogenetic tree was built using these contigs along with reference sequences from previously characterized anelloviruses (Fig. 5). Both reference and de novo assembled contigs can be clustered into three distinct classes: a class of nonhuman infecting reference anelloviruses (green), a set of novel anelloviruses that have only 35–48% sequence similarity with any reference sequence (yellow), and the rest, which we identify as new species or subspecies of existing TTVs. This work has almost doubled the total number of anelloviruses found in humans, including a new potential genus of human-infecting anelloviruses. Therefore, assembly-based metagenomic methods can uncover vast amounts of prevalent diversity even in known viral families with numerous reference sequences.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 5.

Phylogenetic tree of reference and de novo assembled (bold) torque teno viruses (TTVs). Yellow sequences are a class of contigs divergent from known TTVs; green are known TTVs in animals. On the right are blue bars (with numbers) indicating the percentage identity of the nearest (nonself) reference sequence.