BACKGROUND
metagenomics is a new field that studies the microbes under different environmental conditions such as ocean, soil, human distal gut, and many others  <cit> .
using culture-independent genomic sequencing technologies, metagenomics provides a more global and less biased view of an entire microbial community than the traditional isolated genomics.
the earlier metagenomic studies were largely carried out with sanger sequencing, but recently, more studies  <cit>  were performed with the new breaking through next-generation sequencing technologies  <cit> .
today, the pyrosequencing by roche's  <dig> life science serves as a dominant sequencing platform in metagenomics.
however, it is known that the  <dig> sequencers produce artificially duplicated reads, which might lead to misleading conclusions.
exact duplicates sometimes were removed before data analyses  <cit> .
recently, in the study by gomez-alvarez et al  <cit> , nearly identical duplicates, the reads that begin at the same position but may vary in length or bear mismatches, were also classified as artifacts.
exact and nearly identical duplicates may make up 11~35% of the raw reads.
in metagenomics, the amount of reads is used as an abundance measure, so artificial duplicates will introduce overestimation of abundance of taxon, gene, and function.
duplicated reads observed in a pyrosequencing run include not only artificial duplicates but also natural duplicates - reads from the same origin that start at the same genomic position by chance.
therefore, simply removing all duplicates might also cause underestimation of abundance associated with naturally duplicated reads.
the occurrence of natural duplicates can be very low for metagenomic samples lacking dominant species  <cit> , or can be very high for other samples like transcriptomic samples .
so it is important not only to identify the duplicates, but also to distinguish whether they are artificial or natural duplicates.
exactly identical sequences can be easily found, but identification of non-exact duplicates requires sophisticated algorithms to process the massive sequence comparisons between reads.
in gomez-alvarez et al's study  <cit> , the duplicates were identified by first clustering the reads at 90% sequence identity using cd-hit program  <cit>  and then parsing the clustering results.
in this study, we first implemented a method for identification of exact and nearly identical duplicates from pyrosequencing reads.
as the original developers of cd-hit, we reengineered cd-hit into a new program that can process the duplicates more effectively than the original program.
this method can process a typical  <dig> dataset in ~ <dig> minutes; it also provides a consensus sequence for each group of duplicates.
secondly, we validated this method using the underlying raw reads from a list of genome projects utilizing pyrosequencing technology.
we compared the occurrences of the duplicates identified by our method and the natural duplicates made by independent simulations.
lastly, we studied duplicates for several metagenomic samples and estimated their natural duplicates under different situations.
RESULTS
program for identification of duplicates
we implemented a computer program called cdhit- <dig> to identify duplicated reads by reengineering our ultra-fast sequence clustering algorithm cd-hit  <cit> .
the algorithm and other details of cdhit- <dig> are introduced in the methods section.
briefly, we constrained cdhit- <dig> to find exact duplicates and nearly identical duplicates that start at the same position and are within the user-defined level of mismatches .
we allow mismatches in order to tolerate sequencing errors.
the default parameters of mismatches are based on the pyrosequencing error model derived in this study.
we provide a tool in cdhit- <dig> to build a consensus sequence for each group of duplicates.
the model used for consensus generation is also described in the methods section.
cdhit- <dig> software and web server are available at http://cd-hit.org/.
duplicated reads of genomic datasets
we tested and validated cdhit- <dig> with data from pyrosequencing-based genome projects where both the complete genomes and the underlying raw reads are available from ncbi at refseq and short read archive .
for a project with multiple sequencing runs, only the run with the most reads was selected.
we identified  <dig> such genome projects on september  <dig>  with  <dig> from gs- <dig> and  <dig> from gs-flx platform .
the details of the  <dig> projects, including their project identifiers, sra accession numbers, genome sizes, gc contents, and other calculated results, are summarized in table  <dig>
aproject ids, sra study accessions, and sra run accessions are from ncbi short read archive at http://www.ncbi.nlm.nih.gov/traces/sra/sra.cgi.
bread density is the number of reads divided by the genome length.
dσ is the standard deviation, which is based on the results of  <dig> simulations .
ethe platform provided by sra is gs_flx, and the read length  suggests gs_flx titanium.
fthe platform provided by sra is gs_ <dig>  but the read length  suggests gs_flx.
for each genome project, we applied cdhit- <dig> on the raw reads to identify the duplicates, which include both artificial and natural duplicates.
the number of natural duplicates was empirically estimated by applying cdhit- <dig> on simulated reads, which are fragments randomly cut from the complete genomes.
a simulated reads set and the experimental raw reads set in each genome project have the exactly the same number of sequences of exactly the same lengths.
we generated  <dig> sets of simulated reads for each genome project and selected  <dig> sets that are most similar to its corresponding raw reads in gc content.
we further introduced sequencing errors  to the simulated reads according to the error model derived in this study .
these processes made the simulated read sets as similar as possible to the real reads set, except that the former only contains natural duplicates.
using cdhit- <dig>  we identified the duplicates for the  <dig> sets of simulated reads of each project and calculated the average duplicate ratio and the standard deviations.
figure  <dig> shows the ratio of all duplicates and the average natural duplicates for these  <dig> projects.
the results and the standard deviations are also available in table  <dig>
aproject id is the same as in table  <dig>
berror rate is calculated for aligned reads as the number of errors  divided by the number of bases of reads.
aproject id is same as in table  <dig>
berror rate is calculated for aligned reads as the number of errors  divided by the number of bases of reads.
as illustrated in figure  <dig>  the duplicates make up to 4-44% of reads.
we observed that the ratio of natural duplicates, which ranges from 1-11%, highly correlates with the read density  with a pearson correlation factor of  <dig> .
the ratio of artificial duplicates  varies from 3-42%.
on average, artificial duplicates make up 74% of all duplicates, and this percentage varies from 28% to 98% .
artificial duplicates happen randomly without observed correlation with the sample's gc content, genome size, or platform .
here, we define the sensitivity and specificity for the evaluation of duplicate identification.
within the simulated datasets, the reads that start at the same position are considered as true duplicates.
the sensitivity of a method is the ratio of predicted true duplicates by this method to all true duplicates.
the specificity is the ratio of predicted true duplicates to all predictions by this method.
the averaged sensitivity and specificity for the  <dig> datasets are both ~ <dig> % using the default parameters of cdhit- <dig>
pyrosequencing errors
the original  <dig> publication reported an error rate at ~4%  <cit> .
but later studies yielded much higher accuracy.
for example, huse et al.
concluded an error rate at ~ <dig> % for gs <dig> system  <cit> .
quinlan et al.
provided a similar error rate at about ~ <dig> %  <cit> .
accurate estimation of the pyrosequencing error rate is very important for this study, because we use the error rate to optimize the parameters for cdhit- <dig> program to identify the duplicated reads with sequencing errors.
the error model is also used to guide the generation of sequencing errors in the simulated datasets in above analysis.
therefore, we re-evaluated the pyrosequencing error rate using the data from table  <dig>
we used megablast  <cit>  with the default parameters to align the raw reads back to the corresponding reference genomes.
only the reads with at least 90% of the length aligned were selected to calculate the error rates.
other reads were treated as from contamination material, and therefore discarded.
the error rates and fractions by error types  for all projects are shown in table  <dig> for gs- <dig> and table  <dig> for gs-flx.
we found that the error rate  for pyrosequencing is from  <dig> % to  <dig> %.
about 75% of the reads have no error; and about another 20% of the reads have ≤ 2% errors.
if the sequencing error rate is 2%, two reads may have up to 4% mismatches.
we set the default mismatch cutoff at 4% for cdhit- <dig> so that about 95% of reads can be covered.
we examined several mismatch cutoffs from 95% to 98% on the simulated datasets; the 96% cutoff gave the best sensitivity and specificity.
the mismatch cutoff parameter in cdhit- <dig> is a user-configurable parameter.
if the low-quality reads are already filtered out, a higher cutoff such as 98% may be used.
duplicated reads of metagenomic datasets
we studied the pyrosequencing reads for  <dig> metagenomic datasets  of different environments from ncbi sra or from camera metagenomic project http://camera.calit <dig> net.
we identified their duplicates with cdhit- <dig> .
the duplicates make up 5-23% of reads.
as concluded earlier in this study, the quantity of natural duplicates of metagenomic samples depends on the read density of their individual species, and therefore can vary significantly.
since the exact species composition and genome sequences are unknown for metagenomic samples, we could not calculate the amount of natural duplicates as we did for the genome projects.
so, we simulated the occurrence of natural duplicates under several hypothetical sample types.
adatasets are either from ncbi short read archive with project ids and run accession numbers at http://www.ncbi.nlm.nih.gov/traces/sra/sra.cgi or from camera with project and sample names at http://camera.calit <dig> net.
bhigh-, cmoderate-complexity microbial  environment with average genome length of  <dig> mb,  <dig> kb, and  <dig> kb
metagenomic samples are roughly grouped into low-, moderate-, and high-complexity samples, which represent samples dominated by a single near-clonal population, samples with more than one dominant species, and those lacking dominant populations respectively <cit> .
we constructed  <dig> hypothetical sample types: m- <dig> mb, m- <dig> kb, m- <dig> kb, h- <dig> mb, h- <dig> kb, and h- <dig> kb; here the name of a sample type starts with m or h  followed by the average genome size.
therefore, m- <dig> mb represents a moderate-complexity microbial sample with  <dig> mb genomes; and h- <dig> kb may represent a high-complexity small viral sample.
we assumed that each hypothetical sample contained  <dig> different genomes of certain length.
given a set of reads, in order to calculate the natural duplicates under a high-complexity hypothetical sample type, we assigned these reads to the  <dig> genomes randomly at arbitrary positions on either strand.
for a moderate-complexity type, 50% of the reads were randomly assigned to  <dig> dominant genomes, and other reads were randomly mapped to the remaining  <dig> low-abundance genomes.
the natural duplicates were identified by comparing the mapping coordinates.
we applied this method to the  <dig> metagenomic datasets to calculate the ratio of their natural duplicates under different hypothetical sample types .
from figure  <dig>  we can see that if these metagenomic samples match h- <dig> mb, m- <dig> mb, and h- <dig> kb, their natural duplicates only make up  <dig> %,  <dig> %, and  <dig> % of all duplicates on average; so it is appropriate to filter out the duplicates.
however, if a metagenomic sample matches other types, the number of anticipated natural duplicates may be similar or even larger than the artificial duplicates.
here, we want to discuss a particular type of samples - metatranscriptomic samples.
similar to metagenomics, metatranscriptomics is a field that studies the microbial gene expression via sequencing of the total rnas extracted directly from environments.
recent metatranscriptomic studies  <cit>  were performed with pyrosequencing.
since most microbial transcript sequences are only 102~ <dig> bases in length and one transcript can have many copies in a cell, so the read density of metatranscriptomic samples is several orders of magnitude higher than the read density of metagenomic samples.
it is expected that metatranscriptomic samples have high occurrence of natural duplicates.
for example, 61% of reads in the mrna samples from  <cit>  are found as duplicates by cdhit-454; it is reasonable to believe most of them are natural duplicates and therefore should be kept for abundance analysis.
CONCLUSIONS
in this study, we present an effective method to identify exact and nearly identical duplicated sequences from pyrosequencing reads.
but since the identified duplicates contain natural duplicates, it is important to estimate the proportion of natural duplicates.
in the cdhit- <dig> package, we provide a tool to estimate the number of natural duplicates under any hypothetical sample type defined by users, so users can decide whether to retain or remove duplicates in their projects.

 metagenomics is a new field that studies the microbes under different environmental conditions such as ocean, soil, human distal gut, and many others. using culture-independent genomic sequencing technologies, metagenomics provides a more global and less biased view of an entire microbial community than the traditional isolated genomics. the earlier metagenomic studies were largely carried out with sanger sequencing, but recently, more studies were performed with the new breaking through next-generation sequencing technologies. today, the pyrosequencing by roche's life science serves as a dominant sequencing platform in metagenomics. however, it is known that the sequencers produce artificially duplicated reads, which might lead to misleading conclusions. exact duplicates sometimes were removed before data analyses. recently, in the study by gomez-alvarez et al, nearly identical duplicates, the reads that begin at the same position but may vary in length or bear mismatches, were also classified as artifacts. exact and nearly identical duplicates may make up 11~35% of the raw reads. in metagenomics, the amount of reads is used as an abundance measure, so artificial duplicates will introduce overestimation of abundance of taxon, gene, and function. duplicated reads observed in a pyrosequencing run include not only artificial duplicates but also natural duplicates - reads from the same origin that start at the same genomic position by chance. therefore, simply removing all duplicates might also cause underestimation of abundance associated with naturally duplicated reads. the occurrence of natural duplicates can be very low for metagenomic samples lacking dominant species, or can be very high for other samples like transcriptomic samples. so it is important not only to identify the duplicates, but also to distinguish whether they are artificial or natural duplicates. exactly identical sequences can be easily found, but identification of non-exact duplicates requires sophisticated algorithms to process the massive sequence comparisons between reads. in gomez-alvarez et al's study, the duplicates were identified by first clustering the reads at 90% sequence identity using cd-hit program and then parsing the clustering results. in this study, we first implemented a method for identification of exact and nearly identical duplicates from pyrosequencing reads. as the original developers of cd-hit, we reengineered cd-hit into a new program that can process the duplicates more effectively than the original program. this method can process a typical dataset in ~ minutes; it also provides a consensus sequence for each group of duplicates. secondly, we validated this method using the underlying raw reads from a list of genome projects utilizing pyrosequencing technology. we compared the occurrences of the duplicates identified by our method and the natural duplicates made by independent simulations. lastly, we studied duplicates for several metagenomic samples and estimated their natural duplicates under different situations.  program for identification of duplicates we implemented a computer program called cdhit- to identify duplicated reads by reengineering our ultra-fast sequence clustering algorithm cd-hit. the algorithm and other details of cdhit- are introduced in the methods section. briefly, we constrained cdhit- to find exact duplicates and nearly identical duplicates that start at the same position and are within the user-defined level of mismatches. we allow mismatches in order to tolerate sequencing errors. the default parameters of mismatches are based on the pyrosequencing error model derived in this study. we provide a tool in cdhit- to build a consensus sequence for each group of duplicates. the model used for consensus generation is also described in the methods section. cdhit- software and web server are available at http://cd-hit.org/. duplicated reads of genomic datasets we tested and validated cdhit- with data from pyrosequencing-based genome projects where both the complete genomes and the underlying raw reads are available from ncbi at refseq and short read archive. for a project with multiple sequencing runs, only the run with the most reads was selected. we identified such genome projects on september with from gs- and from gs-flx platform. the details of the projects, including their project identifiers, sra accession numbers, genome sizes, gc contents, and other calculated results, are summarized in table aproject ids, sra study accessions, and sra run accessions are from ncbi short read archive at http://www.ncbi.nlm.nih.gov/traces/sra/sra.cgi. bread density is the number of reads divided by the genome length. dσ is the standard deviation, which is based on the results of simulations. ethe platform provided by sra is gs_flx, and the read length suggests gs_flx titanium. fthe platform provided by sra is gs_ but the read length suggests gs_flx. for each genome project, we applied cdhit- on the raw reads to identify the duplicates, which include both artificial and natural duplicates. the number of natural duplicates was empirically estimated by applying cdhit- on simulated reads, which are fragments randomly cut from the complete genomes. a simulated reads set and the experimental raw reads set in each genome project have the exactly the same number of sequences of exactly the same lengths. we generated sets of simulated reads for each genome project and selected sets that are most similar to its corresponding raw reads in gc content. we further introduced sequencing errors to the simulated reads according to the error model derived in this study. these processes made the simulated read sets as similar as possible to the real reads set, except that the former only contains natural duplicates. using cdhit- we identified the duplicates for the sets of simulated reads of each project and calculated the average duplicate ratio and the standard deviations. figure shows the ratio of all duplicates and the average natural duplicates for these projects. the results and the standard deviations are also available in table aproject id is the same as in table berror rate is calculated for aligned reads as the number of errors divided by the number of bases of reads. aproject id is same as in table berror rate is calculated for aligned reads as the number of errors divided by the number of bases of reads. as illustrated in figure the duplicates make up to 4-44% of reads. we observed that the ratio of natural duplicates, which ranges from 1-11%, highly correlates with the read density with a pearson correlation factor of. the ratio of artificial duplicates varies from 3-42%. on average, artificial duplicates make up 74% of all duplicates, and this percentage varies from 28% to 98%. artificial duplicates happen randomly without observed correlation with the sample's gc content, genome size, or platform. here, we define the sensitivity and specificity for the evaluation of duplicate identification. within the simulated datasets, the reads that start at the same position are considered as true duplicates. the sensitivity of a method is the ratio of predicted true duplicates by this method to all true duplicates. the specificity is the ratio of predicted true duplicates to all predictions by this method. the averaged sensitivity and specificity for the datasets are both ~ % using the default parameters of cdhit- pyrosequencing errors the original publication reported an error rate at ~4%. but later studies yielded much higher accuracy. for example, huse et al. concluded an error rate at ~ % for gs system. quinlan et al. provided a similar error rate at about ~ %. accurate estimation of the pyrosequencing error rate is very important for this study, because we use the error rate to optimize the parameters for cdhit- program to identify the duplicated reads with sequencing errors. the error model is also used to guide the generation of sequencing errors in the simulated datasets in above analysis. therefore, we re-evaluated the pyrosequencing error rate using the data from table we used megablast with the default parameters to align the raw reads back to the corresponding reference genomes. only the reads with at least 90% of the length aligned were selected to calculate the error rates. other reads were treated as from contamination material, and therefore discarded. the error rates and fractions by error types for all projects are shown in table for gs- and table for gs-flx. we found that the error rate for pyrosequencing is from % to %. about 75% of the reads have no error; and about another 20% of the reads have ≤ 2% errors. if the sequencing error rate is 2%, two reads may have up to 4% mismatches. we set the default mismatch cutoff at 4% for cdhit- so that about 95% of reads can be covered. we examined several mismatch cutoffs from 95% to 98% on the simulated datasets; the 96% cutoff gave the best sensitivity and specificity. the mismatch cutoff parameter in cdhit- is a user-configurable parameter. if the low-quality reads are already filtered out, a higher cutoff such as 98% may be used. duplicated reads of metagenomic datasets we studied the pyrosequencing reads for metagenomic datasets of different environments from ncbi sra or from camera metagenomic project http://camera.calit net. we identified their duplicates with cdhit-. the duplicates make up 5-23% of reads. as concluded earlier in this study, the quantity of natural duplicates of metagenomic samples depends on the read density of their individual species, and therefore can vary significantly. since the exact species composition and genome sequences are unknown for metagenomic samples, we could not calculate the amount of natural duplicates as we did for the genome projects. so, we simulated the occurrence of natural duplicates under several hypothetical sample types. adatasets are either from ncbi short read archive with project ids and run accession numbers at http://www.ncbi.nlm.nih.gov/traces/sra/sra.cgi or from camera with project and sample names at http://camera.calit net. bhigh-, cmoderate-complexity microbial environment with average genome length of mb, kb, and kb metagenomic samples are roughly grouped into low-, moderate-, and high-complexity samples, which represent samples dominated by a single near-clonal population, samples with more than one dominant species, and those lacking dominant populations respectively. we constructed hypothetical sample types: m- mb, m- kb, m- kb, h- mb, h- kb, and h- kb; here the name of a sample type starts with m or h followed by the average genome size. therefore, m- mb represents a moderate-complexity microbial sample with mb genomes; and h- kb may represent a high-complexity small viral sample. we assumed that each hypothetical sample contained different genomes of certain length. given a set of reads, in order to calculate the natural duplicates under a high-complexity hypothetical sample type, we assigned these reads to the genomes randomly at arbitrary positions on either strand. for a moderate-complexity type, 50% of the reads were randomly assigned to dominant genomes, and other reads were randomly mapped to the remaining low-abundance genomes. the natural duplicates were identified by comparing the mapping coordinates. we applied this method to the metagenomic datasets to calculate the ratio of their natural duplicates under different hypothetical sample types. from figure we can see that if these metagenomic samples match h- mb, m- mb, and h- kb, their natural duplicates only make up %, %, and % of all duplicates on average; so it is appropriate to filter out the duplicates. however, if a metagenomic sample matches other types, the number of anticipated natural duplicates may be similar or even larger than the artificial duplicates. here, we want to discuss a particular type of samples - metatranscriptomic samples. similar to metagenomics, metatranscriptomics is a field that studies the microbial gene expression via sequencing of the total rnas extracted directly from environments. recent metatranscriptomic studies were performed with pyrosequencing. since most microbial transcript sequences are only 102~ bases in length and one transcript can have many copies in a cell, so the read density of metatranscriptomic samples is several orders of magnitude higher than the read density of metagenomic samples. it is expected that metatranscriptomic samples have high occurrence of natural duplicates. for example, 61% of reads in the mrna samples from are found as duplicates by cdhit-454; it is reasonable to believe most of them are natural duplicates and therefore should be kept for abundance analysis.  in this study, we present an effective method to identify exact and nearly identical duplicated sequences from pyrosequencing reads. but since the identified duplicates contain natural duplicates, it is important to estimate the proportion of natural duplicates. in the cdhit- package, we provide a tool to estimate the number of natural duplicates under any hypothetical sample type defined by users, so users can decide whether to retain or remove duplicates in their projects. 
but for some other samples like transcriptomic samples, majority of the observed duplicates might be natural duplicates. artificial duplicates from pyrosequencing reads may lead to incorrect interpretation of the abundance of species and genes in metagenomic studies. we implemented a method for identification of exact and nearly identical duplicates from pyrosequencing reads. for high-complexity metagenomic samples lacking dominant species, natural duplicates only make up <1% of all duplicates. we compared the occurrences of the duplicates identified by our method and the natural duplicates made by independent simulations. however, since the duplicated reads observed in a pyrosequencing run also include natural duplicates, simply removing all duplicates may also cause underestimation of abundance associated with natural duplicates. we observed that the duplicates, including both artificial and natural duplicates, make up 4-44% of reads. the number of natural duplicates highly correlates with the samples' read density. it is important not only to identify the duplicates from metagenomic datasets but also to distinguish whether they are artificial or natural duplicates. we provide a tool to estimate the number of natural duplicates according to user-defined sample types, so users can decide whether to retain or remove duplicates in their projects. this method can process a typical dataset in ~ minutes; it also provides a consensus sequence for each group of duplicates. our method is available from http://cd-hit.org as a downloadable program and a web server.
 artificial duplicates from pyrosequencing reads may lead to incorrect interpretation of the abundance of species and genes in metagenomic studies. duplicated reads were filtered out in many metagenomic projects. however, since the duplicated reads observed in a pyrosequencing run also include natural duplicates, simply removing all duplicates may also cause underestimation of abundance associated with natural duplicates.  we implemented a method for identification of exact and nearly identical duplicates from pyrosequencing reads. this method performs an all-against-all sequence comparison and clusters the duplicates into groups using an algorithm modified from our previous sequence clustering method cd-hit. this method can process a typical dataset in ~ minutes; it also provides a consensus sequence for each group of duplicates. we applied this method to the underlying raw reads of genomic projects and metagenomic projects that utilized pyrosequencing technique. we compared the occurrences of the duplicates identified by our method and the natural duplicates made by independent simulations. we observed that the duplicates, including both artificial and natural duplicates, make up 4-44% of reads. the number of natural duplicates highly correlates with the samples' read density. for high-complexity metagenomic samples lacking dominant species, natural duplicates only make up <1% of all duplicates. but for some other samples like transcriptomic samples, majority of the observed duplicates might be natural duplicates.  our method is available from http://cd-hit.org as a downloadable program and a web server. it is important not only to identify the duplicates from metagenomic datasets but also to distinguish whether they are artificial or natural duplicates. we provide a tool to estimate the number of natural duplicates according to user-defined sample types, so users can decide whether to retain or remove duplicates in their projects. 