Three types of metagenomic datasets have been utilized in this study: simulated, ab initio and in vivo (Figure 1; Table 1).

For in silico simulated datasets, reads of target, and contamination species were generated by NeSSM (Jia et al., 2013). In this study, we assume only one target bacteria species is present in each sample. The target bacteria species used in this study include three model organisms: Bacillus subtilis, Staphylococcus aureus, Escherichia coli (dataset A, dataset B, dataset C, Table 1). Their reads were mixed with reads generated from the genome of 5 or 10 representative species in human oral microbial community [referred to as HOB(5/10)], which were used as possible human contaminations (Hasan et al., 2014). Gradient proportions of reads from target species were set to 5, 35, 65, 95%. We also combined Saccharomyces cerevisiae with B. subtilis and 10 oral bacteria to simulate a special condition with eukaryotic contamination (dataset D, Table 1). In each dataset, over 10 million pair-end reads with 100X coverage were generated at the length of 120 bp every 200 bp bin. All other parameters were set as default (Jia et al., 2013).

For ab initio datasets, we mixed the real sequencing data of B. subtilis strain 168 with real metagenomic sequences from human saliva samples (dataset E, Table 1), with the relative proportion of reads from target species (B. subtilis) set at 35, 65, 95% for different datasets. These samples were named AB_BS_35%, AB_BS_65%, AB_BS_95%, respectively. In order to maintain them in their natural state, we did not perform any filtration of any human contamination from saliva samples (Hasan et al., 2014), even though the lack of filtration may adversely affect the contig assembly and clustering process.

For both ab initio and in vivo dataset preparations, real DNA extraction and sequencing were needed. Their procedures are detailed in the sub-section “in vivo sample preparation, DNA extraction and sequencing”.

The in vivo datasets used in this study were metagenomic (not 16s rRNA) datasets from real community samples prepared as follows: after being activated, strain B. subtilis 168 was cultured overnight till its OD₆₀₀ value reaches between 0.6 and 0.8. All the B. subtilis were centrifuged at 12,000 rev min⁻¹ (12114 g) for the following experiments. The fresh saliva was collected from three healthy adults abstaining from drinking water or gargling for about 30 min before sample collection. Then 200 μl fresh saliva was added to the B. subtilis culture before DNA extraction, amplification and sequencing. This sample was named Real_BS (dataset F, Table 1).

Modified CTAB method (Porebski et al., 1997; Cheng et al., 2014a,b) was chosen for obtaining high MW metagenomic DNA of samples. Five milliliter lysis buffer (Cetyl Trimethyl Ammonium Bromide, 1% w/v; EDTA, 100 mM; NaCl, 1.5 mol l⁻¹; Sodium phosphate, 100 mmol l⁻¹; Tris-Cl pH 8.0, 100 mmol l⁻¹) and 20 μl Proteinase K was added into 15 ml liquid sample (B.S. culture or its mixture with saliva), followed by gentle shaking at 100rev min⁻¹. SDS was added to a final concentration of 1% and the reaction was incubated at 65°C for 30 min with intermittent shaking. After the above steps, an equal volume of saturated phenol, chloroform and isoamyl alcohol (25: 24:1) was added to the mixture and centrifuge at 12,000 rev min⁻¹ (12,114 g) for 10 min to collect supernatant that are free from protein. This step is then repeated once. Metagenomic DNA was precipitated with 0.6 volumes of isopropanol for 30 min at −20°C and pelleted by centrifugation at 12,000 rev min⁻¹ (12,114 g) for 10 min. DNA was washed twice with 70% ethanol and finally dissolved into a 200 μl of TE (1X), pH 8.0. The genomic DNA of B. subtilis 168 and the mixture of B. subtilis 168 with human saliva were extracted with Soil Genomic DNA kit(CWBIO).

Before sequencing with Illumina Miseq-2000, DNA samples were quantified using a Qubit® 2.0 Fluorometer (Invitrogen, Carlsbad, CA) and its quality was checked on a 0.8% agarose gel. 5–50 ng metagenomic DNA in high quality was used as the template for amplifying the V3-V4 diameter region of 16S rRNA genes for each individual sample, with ″5′-CCTACGGRRBGCASCAGKVRVGAAT-3^‴ used as the forward primer and ″5′-GGACTACNVGGGTWTCTAATCC-3^‴ as the reverse primer. The sequencing library was constructed using a MetaVxTM Library Preparation kit (GENEWIZ, Inc., South Plainfield, NJ, USA). Then indexed adapters were added to the ends of 16S rDNA amplicons by limited-cycle PCR. Verified by Agilent 2100 Bioanalyzier (Agilent Technologies, Palo Alto, CA, USA), and quantified by Qubit® 2.0 Fluorometer (Invitrogen, Carlsbad, CA) and real-time PCR (Applied Biosystems, Carlsbad, CA, USA), the DNA libraries were normalized for sequencing. All sequencing reactions were performed on the Illumina MiSeq platform using paired-end sequencing technology (2^*300 bp).

The taxonomical profiles were generated by the Parallel-Meta pipeline (version 2.0) (Jia et al., 2013). 16S rRNA sequences were extracted from raw sequencing data through an HMM model, and these sequences were searched against the Greengene database (http://greengenes.secondgenome.com/) to identify their species.

At the contig binning step, the total number of species identified was used to set the initial cluster number, which aims to provide better accuracy for clustering. An additional eukaryotic18S rRNA database was used as reference when processing the dataset with S. cerevisiae. Though choosing the number of clusters become difficult for unknown contaminants whose information are not recorded in 16S RNA or 18S rRNA database, this approach is practical and performs well for target identification and contamination filtration.

Two assemblers were applied to assemble contigs from community reads. One of the assemblers selected was Velvet (Jia et al., 2013), which could deal with de novo genomic assembly and short sequencing reads alignment. For Velvet, we use the velveth command to construct the dataset as preparation work, and velvet command to build the de Bruijn graph from the k-mers obtained by velveth and extract the contigs. We use k = 12, and set other parameters to auto or default. The other one was MEGAHIT (Li et al., 2015), an assembler designed specifically for complex metagenomics via succinct de Bruijn graph. It is worth mentioning that by using these two assemblers, the abundance information has been intrinsically taken into consideration.

For simulated metagenomic datasets, assembly was performed on two assemblers to compare their performance. Basic assembly statistics were extracted and compared. As MEGAHIT has been shown to be superior to Velvet on simulated data, only MEGAHIT was used to process ab initio and in vivo dataset.

Contig binning is another key step of QC-Blind. Out of all the existing and evaluated binning algorithms, CONCOCT is selected because: firstly, both sequence composition and coverage across multiple samples were considered in contig binning; and secondly, it could handle both single sample and multiple samples. These make CONCOCT a suitable choice for batch processing of possibly contaminated samples (Knights et al., 2011). For multiple species, co-assembly is a necessary prerequisite to running CONCOCT, as CONCOCT takes contigs as input to maximize the number of genomes that could potentially be resolved. Contigs are limited to lengths from 1,000 to 10,000 bp, the lower limit filters low quality contigs while the upper limit cuts fragments for more statistical weight. The number of clusters was precisely determined with 16S rRNA method by Parallel-Meta (Su et al., 2014). Contigs would be clustered into species-level groups after CONCOCT processing. Again, we have to emphasize that by using CONCOCT, the abundances (read depth) of contigs in each of the clusters are similar.

Utilizing marker genes for target species cluster identification helps us overcome the situation where we do not have the complete or partial reference genome. This is realistic for most of the contemporary sequencing tasks, in which we only know a few of the targets' marker genes (Davey et al., 2011; Sunagawa et al., 2013).

Except for a few well-studied species that have had their whole genome sequenced, we have knowledge on only a few genes for most of the other species, such as important regulators in their metabolism, proliferation, or 16S rRNA genes. Therefore, we can only utilize these genes as markers to identify contig clusters that belong to the target species. Marker genes were generated by the combination of MetaPhlAn2 (Truong et al., 2015) and manual curation. The majority of these marker genes are obtained from MetaPhlAn2's unique clade-specific marker genes identified from ~17,000 reference genomes (http://huttenhower.sph.harvard.edu/metaphlan2). For those species (newly sequenced) not included in MetaPhlAn2's marker gene sets, we have extracted marker genes manually and ensure they are clade-specific.

The clusters generated in the previous step were mapped to marker genes of target species by BLAST (e-value cutoff = 1e-20). The number of usable marker genes are species-dependent. The more unique the genes, the more specific the identification would be. Through MetaPhlAn2 selection and literature review, marker genes ftsZ, lytF, nsrR, spo0A, ygxB, yjbH, yjbI were selected for B. subtilis (Amann et al., 1995; Lindgreen et al., 2016; Vikram et al., 2016), acpP, casA, cof, dxs, fabB, fabF, leuO, tesA, uidA were chosen for E. coli (Coleman et al., 2015; Katta et al., 2015; Freedman et al., 2016; Yan et al., 2016; Yu et al., 2017).

After marker gene assignment, contigs containing marker genes for target species are identified as belonging to the target species (defined as target contigs). Based on these assignments, raw reads were mapped to contigs identified as belonging to the targets with BOWTIE2 (Langmead and Salzberg, 2012) (defined as target reads).

DS (dominant species) is defined as the species whose contigs (and reads) outnumber other species' in the cluster. In each cluster, all the contigs were mapped to their reference genome database consisted of both target bacteria and contaminations by BLAST, to identify their source species. Reads were mapped to contigs by BOWTIE2 (Langmead and Salzberg, 2012) and thus inheriting taxonomical information from contigs. Purity was defined as the proportion of contigs or reads of the DS in each cluster, after contig binning through CONCOCT (Knights et al., 2011; Formulas 1, 2).

Purity of each cluster was categorized and evaluated at three levels, 100%, 90%+, 80%+. The proportion of clusters exceeding these thresholds among all the clusters reflects the quality of the clustering.

We define clusters that contain contigs (and reads) from target species (TS) as target clusters (TC), no matter how many contigs are in the TC. The proportion of target contig or target read in each clusters, as well as their distributions among all TC are measured in Formulas (3, 4).

Target concentrations in all clusters would thus form a distribution, which is named “Target distribution.” A more biased target distribution indicates a better clustering result, and vice versa.

“Target/Contamination Dichotomization” were also measured by sensitivity and specificity, which provided quantitative information about contaminations correctly removed and target sequences successfully preserved (Formulas 5–8). Here, the contigs or reads of target species in the clusters found by marker genes alignment were considered as true positive (TP), while those of other species were considered as false positive (FP). The contigs or reads mapping on the target genome were considered as ground truth (GT), which contained those identified or unidentified contigs or reads in target clusters.

We performed functional evaluations to examine how much of the target bacteria species sequencing data has been kept after the quality control process, and whether they retain the functional genomics of the target species. The coverage in target genome was calculated at both base and gene levels (Formulas 9, 10). At the base level, the number of bases in the genome that have been mapped for one or more times was measured (mapped base, MB) and compared with the length of the whole genome (total base, TB). Statistical data were summarized after target reads alignment. At the gene level, genes covered by contigs in our target clusters were viewed as mapped genes (MG), while the genes for the all species were considered as total gene (TG). Genes preserved after the processing were identified by the intersect command of BEDTools (Quinlan and Hall, 2010). Obviously, the gene and base coverage for unprocessed simulated data are both 100%.

Target species can be completely identified at genus level, along with 88% of contaminations identified on average (Table 1 in Supplementary File 2). The setting of number of species acquired in this step could increase the purity of clustering, although the unsupervised contig binning method utilized in QC-Blind does not necessarily require this parameter (Knights et al., 2011). However, it was also found that most eukaryotic organisms (usually not contaminations) would not be identified at species level by this method.

For read assembly, results based on Velvet (Zerbino and Birney, 2008) and MEGAHIT (Li et al., 2015) were compared, showing that MEGAHIT would generate less contigs with longer N50 (e.g., for BS 5%, 603 contigs were generated with N50 = 154,200 by MEGAHIT, while 10,667 contigs were assembled by Velvet with N50 = 39,369), thus we deemed contigs from MEGAHIT outperform those from Velvet for downstream analysis (Table 2 in Supplementary File 2).

Purity of cluster evaluation indicated that most of these contigs generated by QC-Blind were dominated by single species. At the read level, 59.6% (28 of 47) of the clusters reached 90% purity on average, and 72.3% clusters reached 80% purity (34 of 47) (Figure 3A). At contig level, more than 50% of the clusters reached 90% purity. For target concentration, each dataset had a single main cluster that contained over 94% target contigs. Noticeably, the dominant clusters in three simulated datasets were of 100% purity, except that purity of Simu_BS_95.0% dominant cluster 8 was 95% (Figure 3B).

Taken together, the contig binning method could resolve single highly concentrated and pure target cluster from multiple species. Considering possible artifacts produced during read mapping on the simulated datasets, we anticipated that the method would actually perform better for real datasets.

The sensitivity and specificity values for target species read assignment of MEGAHIT-processed data were both high (Figure 3C). Sensitivity values were 92.7% on average in four samples, while specificity values of those were even higher for both target contigs and reads: 100% assignment specificity in Simu_BS_5%, Simu_BS_35%, and Simu_BS_65%, showing that the target information in target cluster can be extracted with very few contaminations remaining. However, the sensitivity and specificity evaluation of Velvet-processed data were extremely low at the dataset with 5% target reads (34.3%, compared to 93.5% in MEGAHIT), which raised question on the ability of Velvet to deal with severely contaminated data. Velvet's sensitivity at the contig level was also not optimistic.

This evaluation of sensitivity and specificity for target species read assignments showed the superiority of using MEGAHIT in QC-Blind method. Thus, in the following analyses, we adopted MEGAHIT in the QC-Blind method as default.

The information loss that we generally experience as the target information progress from raw reads, to read assembly, contig binning, and then marker gene mapping, decreases as the proportion of target species increases (Figure 4A). The greatest data loss on read level occurred in marker gene mapping. the proportion of reads loss were up to 5.31, 5.91, 5.32, 4.87% for each simulated dataset with target fraction from 5 to 95%. Samples without a dominant species (e.g., Simu_BS_5%, Simu_BS_35%) encountered difficulty on assigning all reads correctly, as there may not be enough unique reads from them to reconstruct a complete genome for species identification (Knights et al., 2011). Certain proportions of short contigs (6.47, 9.66, 10.07, 9.05%, for 4 samples, respectively) were also filtered out during contig binning (Figure 4B). Data loss at read level is almost negligible (<2%) in reads assembly and contig binning (Figure 4C).

By contrast, analysis results on the data loss issue by Velvet was unsatisfactory, indicating its inability to deal with highly contaminated data. At contig level, over 35% contigs were lost after contig binning and mapping of Velvet-processed data in all four samples. At read level, Velvet performed comparably as MEGAHIT in three samples, except that in Simu_BS_5%, 61.6% reads were lost after contig assembly by Velvet.

QC-Blind results have nearly perfect coverage of genomic information regardless of the different complexities of the samples. Through mapping these reads back to the genome, on average 93.5% of the bases could be covered for one or more times, indicating the potential of this method to reconstruct a complete genome (Figure 5A). More notably, this coverage is consistent among the four samples (94.1, 92.9, 93.2, 93.8%), which suggests that QC-Blind is able to work with samples of different complexities. Similarly, most of the annotated genes are also covered with processed reads (on average 93.8%). The Arginine biosynthesis pathway, a critical pathway for B. subtilis's metabolism (Kunst et al., 1997), was selected as an example for examination (Figure 5B); all the 20 genes that are involved in this pathway could be found in the processed reads (after quality control) of the four datasets (Figure 5C).

The above assessments of QC-Blind based on simulated data have not only demonstrated the possibility but also the high fidelity of the reference-free QC. This performance can be attributed to the fact that the vast majority of target contigs were binned into a single cluster. Certainly, the selection of marker genes was very crucial, as their uniqueness among microbial community would assure high specificity in target/contamination classification.

Hence for simulated data, the resultant near-perfect coverage proved that the additional work performed on screening is worthwhile. We have put other simulation settings and results in Supplementary File 1 including analytical results of Datasets B, C, and D (Figure 1A).

Before bacterial contamination screening, QC-Blind was able to capture the genetic information of target species in ab initio and in vivo datasets that were contaminated by large proportion of human-oriented reads (416,679,339 reads were identified from bacteria floras of human saliva in this study). This might affect the contig assembly and clustering process.

Single dominant cluster for the target species B. subtilis were successfully determined in each ab initio and in vivo dataset, with cluster 34 in AB_BS 35% containing 99.6% target reads, cluster 33 in AB_BS 65% containing 99.9% target reads, cluster 32 in AB_BS 95% containing 99.5% target reads and cluster 14, 65, 78 together containing 59.7% of target reads (Figure 6A), while a lot of contigs from contamination species with very few sequences were not classified into independent cluster in AB_BS 65% and AB_BS 95%. All three dominant clusters were identified by marker genes with high specificity (Figure 6B). However, the sensitivity of AB_BS 65% at read level dropped to 47.5% while the sensitivity of AB_BS 35% and AB_BS 95% and Real_BS remained high.

For data loss ratios, our method's performance on read datasets remained high level (Figure 7A), except that <20% contigs remained in AB_BS 65% with <30% reads after assembly (Figures 7B,C). For possible explanation for this abnormal phenomenon, we found that the N50 and average contig length of AB_BS 65% were the lowest among the three (Table 2), a large number of its contigs were filtered due to the 600 bp cutoff of CONCOCT. This stringent cutoff was set to remove low quality reads and keep specificity at high level, which if set lower, might recover the data loss (Supplementary File 3). No more than 6% of reads were lost in marker gene mapping.

For base- and gene-level coverage, the performance of QC-Blind on real datasets was consistent with that on simulated datasets, as the analytical results of AB_BS 35%, AB_BS 95% and Real_BS datasets all show coverage of over 98%, except for AB_ BS 65% (Figure 8A). The results reaffirmed the potential of QC-Blind to reconstruct genome from real sequencing data with contaminations (Figures 8B,C).

There is no existing reference-free QC tool similar to QC-Blind up till now, so it is difficult to conduct a completely fair performance comparison with QC-Blind. Here we briefly describe three reference-based methods, followed by a simple comparison study.

Firstly, Kontaminant is a k-mer-based contamination screening tool (Ramirez-Gonzalez, 2013), which has proved to be effective in host filtering for novel viral discovery. Compared to QC-Blind, Kontaminant is limited by the completeness of existing k-mer database, which lead to its inability to work on unknown contaminations. Another general QC analytical tool is FastQC, which is able to generate a comprehensive report on the quality profile of the reads (Andrews, 2010; Brown et al., 2017). With respect to contaminant identification, it could detect the overrepresented sequence through assessing per base sequence content, per base GC content and per sequence GC content, which can be used as evidence that the library is contaminated (Andrews, 2010). Although FastQC is able to detect unknown contaminations, it cannot remove contaminants well, especially in situations where there are multiple contaminating species. Compared to FastQC, QC-blind implements species-based contigs binning method (ideally one cluster is from one species), which performs better in the case that the contaminant is a mixture of multiple species. The third reference-based method is QC-Chain (Zhou et al., 2013), which has good performance in identifying and removing contaminations. However, the false positive rate remains high. Compared to this method, the purity of target clusters in QC-Blind is much higher, and marker genes can map the target clusters accurately, which could ensure higher true positive rate.

As it is still valuable to compare QC-Blind with commonly used reference-based QC tools in terms of contaminant identification and filtration, while the above-mentioned reference-based QC tools are all based on certain assumptions, here we have conducted a simple yet direct comparison of QC-Blind with reference-based Bowtie2 method (https://sourceforge.net/projects/bowtie-bio/files/bowtie2/). Since Bowtie2 is simply a set of mapping tools, it is obvious that if all reference genomes are provided, Bowtie2 would yield the best specificity for target identification and read assignment.

For a direct comparison with reference-based method, we have utilized reference genome for strain B. subtilis 168, as well as the most abundant oral microbes (HOB10, as detailed in Table 1) for the assessment. The Bowtie2 mapping tools with default parameters was used. Results have shown that based on simulated dataset A (Table 1), when using reference genome for strain B. subtilis 168, almost perfect quantities of reads could be assigned to target successfully, which is not surprising. However, when using HOB10 as reference genome for contaminations, 0.3% of error could still occur. On in vivo datasets, it was observed that 88.6 ± 1.5% read-level assignment specificity could be achieved by using reference genome for strain B. subtilis 168. Such results on in vivo datasets are close to the results of QC-Blind (Figure 6). Therefore, it is still true that reference-based method is superior to reference-free method in specificity. However, the gap between the reference-based method and the reference-free method is much closer when QC-Blind is adopted, especially on in vivo sequencing data.

On a computational server with Intel Xeon CPU E7-4820 v2 (each core 2.00 GHz, 8 core, 16 processors—we used only 1 processor) and 512GB RAM, QC-Blind ran in <12 h on datasets with 40 million paired-end reads, with varying time depending on the sequencing quality and contig number. The greatest proportion of time was consumed on contig assembly and contig binning. Due to the time complexity of the clustering algorithm, the running time increased significantly when working with larger number of contigs, which result from the use of lower contig length cutoff values. Thus, a reasonable cutoff, an improved clustering algorithm, as well as the utilization of multiple processors could be taken into consideration to reach higher efficiency.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.