The species-level, open-reference 16S rRNA NGS reads taxonomy assignment pipeline

1. Raw sequences reads in FASTA format were BLASTN-searched against a combined set of 16S rRNA reference sequences. It consists of MOMD (version 0.1), the HOMD (version 15.2 http://www.homd.org/index.php?name=seqDownload&file&type=R ), HOMD 16S rRNA RefSeq Extended Version 1.1 (EXT), GreenGene Gold (GG) (http://greengenes.lbl.gov/Download/Sequence_Data/Fasta_data_files/gold_strains_gg16S_aligned.fasta.gz) , and the NCBI 16S rRNA reference sequence set (https://ftp.ncbi.nlm.nih.gov/blast/db/16S_ribosomal_RNA.tar.gz). These sequences were screened and combined to remove short sequences (<1000nt), chimera, duplicated and sub-sequences, as well as sequences with poor taxonomy annotation (e.g., without species information). This process resulted in 1,015 from HOMD V15.22, 495 from EXT, 3,940 from GG and 18,044 from NCBI, a total of 25,120 sequences. Altogether these sequence represent a total of 15,601 oral and non-oral microbial species.

The NCBI BLASTN version 2.7.1+ (Zhang et al, 2000) was used with the default parameters. Reads with ≥ 98% sequence identity to the matched reference and ≥ 90% alignment length (i.e., ≥ 90% of the read length that was aligned to the reference and was used to calculate the sequence percent identity) were classified based on the taxonomy of the reference sequence with highest sequence identity. If a read matched with reference sequences representing more than one species with equal percent identity and alignment length, it was subject to chimera checking with USEARCH program version v8.1.1861 (Edgar 2010). Non-chimeric reads with multi-species best hits were considered valid and were assigned with a unique species notation (e.g., spp) denoting unresolvable multiple species.

2. Unassigned reads (i.e., reads with < 98% identity or < 90% alignment length) were pooled together and reads < 200 bases were removed. The remaining reads were subject to the de novo operational taxonomy unit (OTU) calling and chimera checking using the USEARCH program version v8.1.1861 (Edgar 2010). The de novo OTU calling and chimera checking was done using 98% as the sequence identity cutoff, i.e., the species-level OTU. The output of this step produced species-level de novo clustered OTUs with 98% identity. Representative reads from each of the OTUs/species were then BLASTN-searched against the same reference sequence set again to determine the closest species for these potential novel species. These potential novel species were pooled together with the reads that were signed to specie-level in the previous step, for down-stream analyses.

Reference:
Edgar RC. Search and clustering orders of magnitude faster than BLAST. Bioinformatics. 2010 Oct 1;26(19):2460-1. doi: 10.1093/bioinformatics/btq461. Epub 2010 Aug 12. PubMed PMID: 20709691.

3. Designations used in the taxonomy:

	1) Taxonomy levels are indicated by these prefixes:
	
	   k__: domain/kingdom
	   p__: phylum
	   c__: class
	   o__: order
	   f__: family
	   g__: genus  
	   s__: species
	
	   Example: 
	
	   k__Bacteria;p__Firmicutes;c__Clostridia;o__Clostridiales;f__Lachnospiraceae;g__Blautia;s__faecis
		
	2) Unique level identified – known species:
	   
	   k__Bacteria;p__Firmicutes;c__Clostridia;o__Clostridiales;f__Lachnospiraceae;g__Roseburia;s__hominis
	
	   The above example shows some reads match to a single species (all levels are unique)
	
	3) Non-unique level identified – known species:

	   k__Bacteria;p__Firmicutes;c__Clostridia;o__Clostridiales;f__Lachnospiraceae;g__Roseburia;s__multispecies_spp123_3
	   
	   The above example “s__multispecies_spp123_3” indicates certain reads equally match to 3 species of the 
	   genus Roseburia; the “spp123” is a temporally assigned species ID.
	
	   k__Bacteria;p__Firmicutes;c__Clostridia;o__Clostridiales;f__Lachnospiraceae;g__multigenus;s__multispecies_spp234_5
	   
	   The above example indicates certain reads match equally to 5 different species, which belong to multiple genera.; 
	   the “spp234” is a temporally assigned species ID.
	
	4) Unique level identified – unknown species, potential novel species:
	   
	   k__Bacteria;p__Firmicutes;c__Clostridia;o__Clostridiales;f__Lachnospiraceae;g__Roseburia;s__ hominis_nov_97%
	   
	   The above example indicates that some reads have no match to any of the reference sequences with 
	   sequence identity ≥ 98% and percent coverage (alignment length)  ≥ 98% as well. However this groups 
	   of reads (actually the representative read from a de novo  OTU) has 96% percent identity to 
	   Roseburia hominis, thus this is a potential novel species, closest to Roseburia hominis. 
	   (But they are not the same species).
	
	5) Multiple level identified – unknown species, potential novel species:
	   k__Bacteria;p__Firmicutes;c__Clostridia;o__Clostridiales;f__Lachnospiraceae;g__Roseburia;s__ multispecies_sppn123_3_nov_96%
	
	   The above example indicates that some reads have no match to any of the reference sequences 
	   with sequence identity ≥ 98% and percent coverage (alignment length)  ≥ 98% as well. 
	   However this groups of reads (actually the representative read from a de novo  OTU) 
	   has 96% percent identity equally to 3 species in Roseburia. Thus this is no single 
	   closest species, instead this group of reads match equally to multiple species at 96%. 
	   Since they have passed chimera check so they represent a novel species. “sppn123” is a 
	   temporary ID for this potential novel species.