BACKGROUND
during the last decade, metagenomics became one of the explosively developing areas of molecular genomics.
advent of the next-generation sequencing allowed performing genomic analysis of samples obtained directly from the environment.
such an approach provides data for an extensive quantitative examination of the microbial community structure, particularly including uncultivable and previously undiscovered components.
the sphere of metagenomic analysis has extended from science to heavy industry  <cit> , agriculture  <cit>  and healthcare  <cit> .
a large amount of metagenomic data is constantly being accumulated which leads to a demand in the means of efficient analysis  <cit> .
one of the common steps in metagenomic study is calculation of pairwise dissimilarity between the samples   <cit> .
beta-diversity is a quantitative measure of the differences in composition between two microbial communities.
its value is calculated from features like taxonomic or functional composition, phylogenetic structure of the whole community, etc.
a dissimilarity matrix composed of pairwise distances between all samples is used for further cluster analysis, classification and study of influence of the experimental factors.
in large-scale studies involving tens and hundreds of metagenomes, critical requirements in beta-diversity analysis include high algorithm performance and low memory usage.
for a long time, the standard technological approach for the evaluation of beta-diversity was based on the identification of species in metagenomic samples through 16s rrna gene sequencing.
however, this method has inherent disadvantages including incompleteness of reference databases, presence of multiple copies of 16s rrna gene in the same genome, discrepancy between phylogenetic trees constructed using 16s rrna and the other genes and lack of information about the other genes and subsequently metabolic potential of the studied community.
an alternative, more informative method is whole-genome sequencing  generating millions of reads from the total dna of the genomes of all organisms inhabiting the environment.
the identification of the organisms in the short-read wgs metagenome is commonly based on the alignment or de novo assembly  <cit> .
the alignment method is a comparison between sequences of obtained reads and sequences of reference genes or genomes, and has significant drawbacks such as high computational costs and incomplete databases.
de novo assembly is usually a time-consuming task for such complex data as metagenomes that may contain many unknown or highly similar genomes of organisms with widely varying abundance.
with a rapid increase in data output produced by sequencing technologies, efficient methods for genomic analysis based on k-mer composition analysis emerged.
such algorithms work with k-mers  obtained directly from metagenomic reads, without pre-mapping or assembly.
in comparison with reference-based methods, the main advantages of k-mer based approaches are compressed representation of sequences and inclusion of the entire data volume into analysis .
among these methods, the most simple and effective for exploratory analysis of large data sets is comparison of sequences by calculation of pairwise dissimilarity between them on the basis of k-mer spectrum - a normalized vector of frequencies of occurrences of each k-mer in the metagenome.
the k-mer length is a key factor influencing specificity and efficiency of the comparison.
for different intervals of k, the respective algorithms have been designed that target different specific tasks.
for example, for very short k-mers  only “rough” estimates are possible: sequence quality check  <cit> , taxonomic separation of individual genomes  or comparison of metagenomic communities with notably different composition .
for k=15− <dig>  the computational costs associated with the processing of the whole spectrum increase significantly.
two approaches can be applied to reduce them.
this problem can be solved in two ways.
the first is to select a fraction of k-mers that describe the studied data in the most complete way   <cit> .
another way is to combine multiple k-mers into one feature using a certain principle   <cit> .
in the intermediate range , it is still computationally feasible to analyze the complete set of k-mers, and the k-mer length remains sufficiently specific  <cit>  for comparing distinct genomes.
among the metagenomic k-mer methods, the most simple and effective for exploratory analysis of large data sets is comparison of sequences by calculation of pairwise distances between them on the basis of k-mer spectra  <cit> .
in this area, researchers are actively designing fast algorithms for calculating and assessing the k-mer spectrum  <cit> .
most studies are focused on examining clustering of the samples by one or more factors   <cit> .
since the prevalent approaches for assessment of beta-diversity today are reference-based, an important question is how their results correlate with the k-mer methods.
in this study, we compared common reference-based methods  with the k-mer approach.
we explored how various characteristics of the data influence the results of k-mer spectra analysis and identified the advantages of k-mer analysis comparing with the reference-based approaches.
to evaluate the applicability of k-mer-based dissimilarity, metagenome of human gut microbiota was selected, the study of which has great biomedical importance and perspective.
although nowadays intestinal microbiota is one of the most studied among complex microbial communities, many of its components are still not fully identified   <cit> .
the application of our method revealed significant presence of one of such components - phage - that went undetected by reference-based methods.
methods
simulated metagenomes
two set of “shotgun” gut metagenomes were simulated using metasim  <cit> .
the high-diversity set included  <dig> metagenomes generated from the genomes of ten distantly related major bacterial species accounting for more than  <dig> % of all reads in chinese group: akkermansia muciniphila atcc baa- <dig>  alistipes shahii wal  <dig>  bacteroides vulgatus atcc  <dig>  bifidobacterium adolescentis atcc  <dig>  coprococcus sp.
art55/ <dig>  eubacterium eligens atcc  <dig>  faecalibacterium prausnitzii l2- <dig>  lachnospiraceae bacterium  <dig>  <dig> 56faa, prevotella copri dsm  <dig> and ruminococcus sp.
18p <dig>
the simulation included the following steps.
first, for each genome, mean and standard deviation of its relative abundance were estimated from the taxonomic composition of the chinese metagenomes.
for each metagenome, ten abundance values were randomly generated under normal distribution with these parameters and the obtained values were normalized to  <dig> million reads; a total of  <dig> genera abundance vectors were obtained .
the metagenomes were generated by mixing ten bacterial genomes at the obtained abundance levels and sampling short reads from the genomes using metasim with read length  <dig> bp.
also we performed sampling of these reads with errors .
the low-diversity simulated group included  <dig> metagenomes generated in a similar way from the genomes of ten closely related major bacterial species accounting for more than  <dig> % of all reads in the hmp group: bacteroides vulgatus atcc  <dig>  bacteroides dorei  <dig>  <dig> 36/d <dig>  bacteroides uniformis atcc  <dig>  bacteroides stercoris atcc  <dig>  bacteroides caccae atcc  <dig>  bacteroides ovatus , bacteroides xylanisolvens xb1a and bacteroides thetaiotaomicron vpi- <dig>  bacterial proportions for these simulations are listed in additional file 1: table s <dig>
for single nucleotide polymorphism  simulations, the same ten reference genomes and abundance values as in the high-diversity dataset were used.
two different models of snps introduction were used: “independent” and “phylogenetic”.
in the “independent” snp model,  <dig> “mutated” genomes were generated for each reference species by changing nucleotide letter at random positions independently with  <dig>  % substitution rate.
thus, the average amount of snps between any two of the “mutated” genomes was ∼  <dig> %.
in “phylogenetic” snp model, the procedure was performed in iterations for each reference genome:
initialize with a single genome; iteration number =  <dig>
make a copy of each of the genomes available at the step.
introduce snps to all genomes at random positions.
increment iteration number.
if the iteration number is greater than  <dig>  stop; else return to step b.
after the  <dig> iterations, 26= <dig> “mutated” genomes are obtained.
in each model, the random “mutated” genomes of corresponding bacteria were used to generate metagenomes the same way as for high-diversity simulation above.
real metagenomic datasets
two “shotgun” gut metagenomic datasets were analyzed:  <dig> metagenomes of healthy usa population  <cit>   and  <dig> metagenomes of chinese population  <cit>  including healthy and type  <dig> diabetes individuals .
for each sample, the reads were filtered by quality using fastq quality filter script from fastx-toolkit  <cit>  .
for each metagenome,  <dig> million of high-quality reads was sampled using random _records script from  <cit> .
comparison of various sampling sizes showed that the selected size of subsampling does not significantly affect the results of the measures’ comparison .
calculation of k-mer vectors and dissimilarity measures
for each metagenome, k-mer spectrum was calculated using an ad hoc java program that processes fasta files read-wise by obtaining k-mer counts for each read and adding the counts to a global array .
after processing all reads, the counts for reverse-complementary k-mers were summed and normalized to a sum of  <dig>  the length of the final feature vector  did not exceed n=22k− <dig> for odd k and n=22k−1+2k− <dig> for even k because of reverse-complement k-mers.
the obtained k-mer spectra were used to calculate pairwise dissimilarity via bray-curtis measure defined as:
bc=1−2∑i=14kmin,mi)∑i=14k+mi)
where m is the vector of k-mer frequencies normalized to a sum of  <dig> per metagenome and x, y are two different metagenomes.
bc =  <dig>  if the frequencies are equal for all k-mers between the metagenomes, and bc =  <dig>  if no common k-mers are present in the metagenomes.
beta-diversity analysis using reference-based methods
taxonomic profiling via mapping of metagenomes to a reference genome catalog and coverage analysis was performed as described previously  <cit> , with the only difference: a non-redundant set of  <dig> genomes of gut microbes was used .
the final feature vector for each metagenome included relative abundance of microbial species was normalized to a sum of  <dig> % .
dissimilarity was calculated using these vectors both with bray-curtis measure and whole-genome adaptation of the weighted unifrac metric  <cit> .
functional profiling was performed as described previously  <cit>  to yield cog   <cit>  relative abundance vectors subsequently used for dissimilarity analysis using bray-curtis measure .
an alternative method of taxonomic profiling employed metaphlan v <dig> .
<dig> ); here mapping was performed using bowtie <dig> v <dig> .
<dig> software  <cit> , up to  <dig> mismatches per read were allowed .
all reference-based methods were summarized in table  <dig>
crassphage abundance analysis
all crassphage genes   <cit>  were aligned to the reference gene catalog .
for each gene, its relative abundance was estimated as a ratio of the total length of the reads mapped to this gene to the total length of the reads mapped to the reference gene catalog.
phage relative abundance was determined as a sum of the relative abundance values of its genes .
as an additional method of metagenomic classification, diamond aligner  <cit>  was used  in combination with megan classifier  <cit> .
statistical analysis
statistical analysis was implemented in r. the code is available at:
http://download.ripcm.com/dubinkina_2015_suppl_data/ and https://github.com/zireae1/kmer_project/.
ethical approval
the sampling procedure was approved by the ethical committee for clinical research from the peking university shenzhen hospital, shenzhen second people’s hospital and medical research center of guangdong general hospital , and enrollments were approved by the institutional review boards of the two recruitment centers  .
RESULTS
to compare the k-mer based metagenomic beta-diversity measure with traditional reference-based methods we conducted a series of computational experiments on simulated and real data.
to validate the k-mer measure and find an optimal value of k, we simulated metagenomes from prevalent human gut bacterial genomes.
then the method was applied for the analysis of a group of real human gut metagenomes sequenced in two large-scale projects: china population   <cit>  and hmp   <cit> .
comparison of beta-diversity measures for simulated metagenomes
there is considerable variation between k-mer spectra for genomes of distinct bacterial species due to the differences in the gene content, amino acid coding preferences, etc.
<cit> .
supposedly, in the case of a metagenome including sufficiently covered genomes of multiple species, one should observe significant accordance between k-mer spectrum of the metagenome and its taxonomic composition.
to verify this hypothesis, we simulated several datasets with different degrees of community richness and applied the bray-curtis measure  for both taxonomic and k-mer profiles to compare the two respective dissimilarity matrices .
simulation 1: high-diversity communities
the first synthetic dataset included  <dig> metagenomes generated by randomly sampling “reads” from ten genomes of common phylogenetically distant human gut bacteria .
comparison of the two approaches showed that, as k increases, so does the correlation value between the two dissimilarity matrices based on k-mers and taxonomic composition.
with high values of k, the two matrices are highly similar .
simulation 2: sequencing errors and snps
besides the considerable number of species within microbiota, the other factors contributing to the diversity of metagenomic k-mers are presence of mutations and sequencing errors in reads.
therefore, we conducted two experiments by introducing artificial snps into genomes and, separately, random single-nucleotide changes  into the reads in order to explore their influence on the correlation of beta-diversity estimates using k-mer and taxonomic methods.
datasets from simulation  <dig> were used here.
for sequencing errors modeling, the reads for each metagenome were simulated with per-nucleotide substitution rate of  <dig> % .
introduction of such “errors” did not lead to a significant change in correlation between the two methods .
for snp modeling, bacterial genomes with  <dig> % of randomly introduced single-nucleotide substitutions  were used to generate simulated metagenomes with the same abundance proportions as in simulation  <dig>  we employed two different models of snps introduction - “independent” and “phylogenetic”.
with the former simulation being more straight-forward, “phylogenetic” approach was developed to model the accumulation of mutations in bacterial species during evolutionary process .
the results of simulations showed that, independent of the model choice, in general snps had a minor effect on k-mer spectra comparable to the effect of simulated sequencing errors: correlation between the k-mer and taxonomic methods decreased from  <dig>  to  <dig>  for “independent” model and  <dig>  for “phylogenetic” model .
noteworthy, introduction of snps had a more pronounced effect on metagenomes with highly similar taxonomic composition.
this was particularly marked when the snp rate was increased to  <dig> % .
simulation 3: low-diversity communities
the second synthetic dataset included  <dig> metagenomes generated by randomly sampling “reads” from ten genomes of common phylogenetically close human gut bacteria - belonging to the same genus - bacteroides .
the correlation between the methods was found to be lower for such homogeneous community than for a heterogeneous one .
the correlation value tends to increase with k but does not achieve the level of simulation  <dig> .
it suggests that higher values of k should be used to increase accuracy; however, the size of the feature vector increases as 4k, hence the computational time quickly becomes unacceptable.
to select the optimal k value, we evaluated the correlation between k-mer and taxonomic dissimilarity matrices together with the computational time of k-mer matrix generation for k=5− <dig> using both high- and low-diversity simulated datasets .
as the results in both simulations showed, with k =  <dig> the dissimilarity matrices are highly correlated while the computational time is still acceptable  took about  <dig> seconds to calculate the spectra) further statistical analysis - calculation of dissimilarity matrix - took about  <dig> -  <dig> minutes, see additional file 5: figure s3).
at the same time, it is the highest value practically acceptable in terms of memory usage: for k= <dig>  the spectra occupied ∼ <dig> gb of memory in r environment, but for k =  <dig> - as much as  <dig> gb.
considering these observations we selected k= <dig> for further analyses.
comparison of beta-diversity measures on real human gut metagenomes
after testing the method on simulated data, we examined two real human gut datasets from large-scale metagenomic projects: china  <cit>  and hmp  <cit> , with the former cohort representing more diverse microbial community structures than the latter  <cit> .
using this data, the pairwise dissimilarity matrix obtained by the k-mer approach with bray-curtis measure  was compared with the dissimilarity matrices obtained by each of the four methods based on taxonomic and functional reference .
to visualize the distributions of beta-diversity values, we applied two types of scatter plots.
the first type is a basic principle coordinate analysis  plot constructed using a single dissimilarity measure, with dots representing distinct metagenomes .
on the second type of plot, two dissimilarity measures are compared: each triangle corresponds to a pair of metagenomes, one measure is plotted against the other .
samples from the two studies  tended to cluster separately by functional, as well as by k-mer composition, but not by taxonomic composition .
therefore, the two cohorts were further analyzed separately.
another interesting fact was that  <dig> of the outliers  present on k-mer scatter plot were also on the periphery of cog scatter plot but not of the taxonomic scatter plot .
these samples were examined in details.
fig.
<dig> variation of metagenomes using different dissimilarity measures.
pcoa plots for different dissimilarity measures: a bc kmer, b bc cog, c wg unifrac, d bc tax , e bc metaphlan .
three samples-outliers are marked with asterisks.
f heatmap of spearman correlation coefficient between dissimilarity matrices obtained using different measures
fig.
<dig> comparison of pairwise difference measures obtained by k-mer and reference-based methods.
for each plot, y-axis represents k-mer distance, x-axis - distance by one of the reference-based methods.
distribution of dissimilarity measures is shown for a bc kmer for all reads and bc tax ; b bc kmer for all reads and bc cog; c bc kmer for reads mapped to the catalog of genomes and bc tax ; d bc kmer for reads mapped to the catalog of genes and bc cog
comparison of the five beta-diversity measures showed that the k-mer measure has a significant similarity with each of the reference-based ones .
the closest was the measure based on cog composition.
for the chinese group, the correlation values tended to be higher than for the hmp group in all comparisons .
the phylogeny-aware metric wg unifrac was among the most dissimilar .
influence of reads mappability
to assess the contribution of the unmapped reads to the results, k-mer spectra were also computed only using the reads that successfully mapped to the corresponding catalog .
this analysis led to interesting results .
first, we observed an equalization of bc tax org vs. bc kmer correlation between the two cohorts .
therefore, fraction of unmapped reads appears to be one of the major factors contributing to the difference between the cohorts.
this parameter is dependent on the representability of the reference catalog and quality of sequencing run.
second, we assessed the shift of each outlier in the direction of the central cloud of points.
quantitatively, for each outlier we calculated the bc kmer difference value: the difference between its bc kmer value and the linearly interpolated middle of the cloud obtained for the same reference-based value .
for comparison with bc tax, the bc kmer difference decreased significantly - from  <dig>  ±  <dig>  to  <dig>  ±  <dig>  .
for comparison with bc cog, the bc kmer difference changed slightly: from  <dig>  ±  <dig>  to  <dig>  ±  <dig> .
correspondingly, a group of pairs-outliers mentioned above moved into the central cloud of points in the bc tax org vs. bc kmer comparison, but did not change their visual location in bc cog vs. bc kmer comparison.
this observation is in agreement with the fact that the gene reference catalog is more complete than the genome reference catalog and the percentage of mapping to the gene catalog is higher .
presumably, the presence of pairs-outliers could be caused by k-mers from certain dominant sequences that are present in the reference base of genes but not genomes.
we investigated these outliers in details.
investigation of samples-outliers
the total human gut metagenome is a phylogenetically diverse structure including not only the sequences of bacterial genomes but also ones from bacterial mobilome , fungi, protozoa, traces of dna of dietary origin, host.
our reference genome catalog partly accounts for such non-microbial components by including the genomes of several common intestinal eukaryotes - clinically relevant yeasts candida  and protozoan blastocystis .
however, many sequences are not present in our genome catalog, particularly viral genomes.
therefore, in our analysis the potential reads of viral origin would not contribute to the taxonomic difference but would change the k-mer spectrum.
recently, a new bacteriophage was discovered - crassphage - shown to be a sole major dominant of the human gut viriome  <cit> .
moreover, its presence was estimated for the hmp metagenomes analyzed in our work: crassphage genome amounts for up to  <dig> % of the reads for this group.
obviously, such a prevalent genome should have a significant influence on k-mer spectra and thus on our comparison of the beta-diversity measures.
basing on the available data on the abundance of crassphage in hmp samples , the cohort was split into two groups - with high phage abundance  and with low phage abundance .
the whole group of extreme outliers was found to consist of the pairs where at least one of the samples belonged to the phage-enriched group .
moreover, the outliers on fig.
1a were also found to be the samples with high fraction of crassphage reads .
fig.
<dig> analysis of samples-outliers.
a distribution of pairwise dissimilarity obtained using k-mer and taxonomic composition for hmp cohort.
different colors indicate groups of dissimilarities for: all hmp pairs, pairs-outliers - where at least one of the samples belonged to the phage-enriched group; cp-filtered pairs - extreme outliers  after removal of k-mers from reads mapped to crassphage  genome; b composition of sample srs <dig> according to the combined results from two analyses
as the extreme outliers were found to be generated by pairs including at least one of the two samples , these samples were analyzed in detail.
the reads that did not map to the genome catalog  were subject to metagenomic classification using an alternative method - using diamond alignment and megan classifier algorithms .
as a result, additionally  <dig> % and  <dig> % of the reads were identified as crassphage .
to further confirm the contribution of high phage fraction to formation of outliers, we subtracted the k-mers of the crassphage reads from k-mer spectra of the samples.
indeed, such operation significantly decreased the k-mer-based dissimilarity for the respective pairs .
discussion
here we have developed an algorithm for assessing pairwise dissimilarity of “shotgun” metagenomes basing on k-mer spectrum and compared it with commonly used reference-based approaches.
the comparison was performed using various measures  on a set of simulated metagenomes as well as on real metagenomes from two large-scale human gut microbiota studies.
for simulated metagenomes, we showed that k =  <dig> is an optimal value in terms of balancing between the resolution of the method and computational time.
this value of k performed well for both high- and low-diversity simulated metagenomes; however, for low-diversity simulations the dissimilarity matrices based on k-mer method and taxonomic composition were less correlated .
this fact was likely due to the decreased diversity of k-mers and thus reduced differentiating resolution.
for real gut metagenomes with complex community structure, the k-mer approach allows to delineate the samples with a wide range of functional composition, as demonstrated on two international cohorts .
on the other hand, we observed that k-mers are less correlated with taxonomic composition than with functional  one.
we speculate that this difference could be associated with significant subspecies-level genomic diversity of gut microbes: a recent analysis of publicly available metagenomic data showed that the average gene variation between individuals across  <dig> abundant species was as high as  <dig> ±  <dig>  %  <cit> .
the k-mer frequencies as well as the gene relative abundance features are sensitive to gene content variation, while in the case of species relative abundance features it would be ignored.
besides gene presence/absence, another common form of genomic variability is snps.
we attempted to model their influence on k-mer beta-diversity.
theoretically, introduction of snps would lead to change in frequencies of k-mers and thus deteriorate correlation between k-mer and reference-based dissimilarity.
in our simulations, when  <dig> % snps were introduced to simulated datasets , the correlation between the methods dropped slightly , independent of whether the evolutionary character of snp accumulation was considered during modeling or not.
however, for real metagenomes the correlation between the methods was lower .
these results suggest the existence of other real-life effects having stronger influence on the correlation than snps .
a major advantage of the k-mer approach is that it exploits the totality of the reads - unlike reference-based methods that inherently discard the reads that failed to map to the reference catalog - and thus the information contained within them.
such a feature promotes the application of the k-mer approach as a tool for assessing the representability of the reference set for given metagenomes.
currently representative sets of reference genomes are available for microbial communities of few environments .
however, recent discoveries imply that the so called reference genomes do not capture a wide intra-species level variation even for this extensively examined community  <cit> .
the situation is even worse for less popular environments - like marine ecosystems  <cit>  or human skin  <cit> : reference catalogs for their microbiota are considerably less complete, thus rendering beta-diversity assessment difficult.
we propose to assess the representability of a reference genome catalog via examining the k-mer content of the metagenomic reads mapped to it.
on the analyzed gut metagenomic data, we observed that k-mer spectra of the mapped reads produced dissimilarity profiles that had higher correlation with those obtained with taxonomic composition than the k-mer spectra of the whole set of reads.
however, lower correlation between the two methods observed for some pairs of samples suggested the presence of dominant genomic sequences not included in the reference catalog.
detailed analysis showed that these outliers corresponded to the hmp samples enriched in crassphage, a recently discovered gut bacteriophage  <cit> ; the genome of this phage was not included in the respective reference catalog.
subtraction of the crassphage k-mers moved the outliers towards the main cloud of points but not into it completely .
presumably, such incomplete compensation can be linked to high level of genomic variability inherent to gut phages  <cit> : originally the consensus crassphage genome was obtained by combined assembly of  <dig> metagenomes from individuals not included in our groups  <cit> , so its sequence in the latter might be quite distant than all the crassphage-related k-mers in our analysis.
additionally, over  <dig> % of the reads remain unidentified by two different methods  they can correspond to genome contributing to formation of outliers.
considering the gene catalog, dedicated analysis of the reads mapped to reference genes did not lead to shift of outliers .
first, a likely reason for this is that the crassphage sequences were included in the catalog: search for crassphage genes in the reference gene catalog  identified highly similar hits for  <dig> of the  <dig> phage genes  that were detected in at least one metagenome.
second, the gut microbial gene catalog was originally constructed from the contigs assembled from total dna reads  <cit>  and is known to contain not only the bacterial genes, but viral and eukaryotic, too.
interestingly, our results also imply that the chinese cohort lacks metagenomes with such high prevalence of this phage, provoking speculations on world-wide phage phylogeography.
while no clinical associations for crassphage have been described to date, omission of phage components could be a significant miss in biomedical studies of microbiota.
there is a growing understanding that gut phages play an important role in the ecology of “phage-gut microbiota-human” system and include potential biomarkers; they are able to transfer clinically important bacterial genes - e.g.
antibiotic resistance and pathogenicity determinants  <cit> .
application of our reference-free k-mer approach can facilitate early detection of such sequences in biomedical diagnostics data and discovery of novel biomarkers.
our approach is not only applicable to metagenomes from an arbitrary environment, but is indispensable for dissimilarity and cluster analysis of communities with poorly described components.
the approach allows to detect a major presence of an unknown organism and/or virus in a metagenome.
we suggest that the approach should be introduced as a necessary method of “shotgun” metagenome composition analysis complementary to reference mapping in order to avoid biases associated with unrepresentative reference database.
although we did not find evidence for outliers caused by technical issues in the examined datasets, the approach can also be used for primary detection of metagenomes with abnormal composition caused by high abundance of host dna , dna of dietary origin  and technical artifacts .
finally, comparison of the metagenomes basing on k-mer spectrum provides more information than mapping to reference sequence catalogs.
essentially, k-mer analysis is a feature extraction procedure applied to metagenomic reads.
the produced set of features  is several orders of magnitude larger than one yielded in reference-based approaches.
therefore, it provides higher discriminative resolution that opens a promising opportunity for developing a new generation of methods for metagenomic analysis, and our method makes a step towards understanding of how to explore such high-dimensional feature space efficiently.
CONCLUSIONS
analysis of k-mer spectra for both simulated and real “shotgun” metagenomes showed that this method allows quick assessment of the pairwise dissimilarity of such datasets.
simulations show that the method is robust to variability introduced by sequencing errors and genomic mutations.
the obtained dissimilarity matrix can be used not only for cluster analysis and classification purposes, but also for early detection of major unknown components and quality control of reference-based approaches.
it is recommended that the method should be included as a complementary step in high-throughput computational pipelines for metagenomic data analysis.
additional files
additional file  <dig>
supplementary tables.
supplementary table s <dig>  bacterial abundances in simulation  <dig> .
supplementary table s <dig>  bacterial abundances in simulation  <dig> .
supplementary table s <dig>  list of genomes in taxonomic catalog for human gut.
supplementary table s <dig>  taxonomic composition for real dataset .
supplementary table s <dig>  taxonomic composition for real dataset .
supplementary table s <dig>  functional composition for real dataset .
supplementary table s <dig>  taxonomic composition for real dataset by metaphlan .
supplementary table s <dig>  mapped read counts and percentage of mapping on taxonomic and functional catalog and phage genome.
additional file  <dig>
figure s <dig>  correlation between k-mer-based dissimilarity matrix obtained using the original full metagenomic readsets and their subsampled versions.
the same amount of reads was repeatedly sampled from  <dig> randomly selected metagenomes .
additional file  <dig>
figure s <dig>  comparison of pairwise difference measures obtained by k-mer and taxonomic composition .
for each pair y-axis represents k-mer dissimilarity, x-axis represents dissimilarity by taxonomic composition.
k= <dig> a) using the original bacterial genomes; b) using the bacterial genomes with  <dig> % of snps.
additional file  <dig>
figure s <dig>  comparison of pairwise dissimilarity measures obtained by k-mer and taxonomic composition for simulated for high- and low-diversity metagenomes.
as seen, satisfactory correlation of k-mers with taxonomic composition can be obtained only at relatively high values of k.
additional file  <dig>
figure s <dig>  correlation between k-mer and taxonomic composition dissimilarity matrices, as well as k-mer dissimilarity matrix computation time with varying values of k. all computations were performed on a compute node with cpu opteron  <dig>  <dig>  ghz  and  <dig> gb ram.
additional file  <dig>
figure s <dig>  comparison of dissimilarity measures obtained by k-mer and  <dig> reference-based methods: bc metaphlan genus, bc metaphlan org and wg unifrac.
for each plot, y-axis represents k-mer dissimilarity, x-axis - dissimilarity using one of reference-based methods.
abbreviations
16s rrna16s ribosomal ribonucleic acid
bcbray-curtis dissimilarity measure
bc kmerbc based on k-mer spectrum
bc taxbc based on taxonomic composition, obtained by mapping to catalog of genomes
bc metaphlanbc based on taxonomic composition, obtained by mapping to clade-specific marker genes
bc cogbc based on functional composition, obtained by mapping to catalog of cog  genes
dnadeoxyribonucleic acid
snpsingle nucleotide polymorphism
wgswhole genome sequencing, “shotgun”
wg unifracwhole genome adaptation of unifrac measure
competing interests
the authors declare that they have no competing interests.
authors’ contributions
dga and avt initiated the project and designed the study.
viu wrote the programs for k-mer counting.
vbd wrote the programs for dissimilarity matrix calculation, generated the simulation datasets and analyzed the data.
dsi performed statistical modeling of snps effect.
dga and avt guided the comparison of dissimilarity measures obtained by k-mer and reference-based methods.
vbd and avt wrote the manuscript under the supervision of dga.
all authors read and approved the final manuscript.

 during the last decade, metagenomics became one of the explosively developing areas of molecular genomics. advent of the next-generation sequencing allowed performing genomic analysis of samples obtained directly from the environment. such an approach provides data for an extensive quantitative examination of the microbial community structure, particularly including uncultivable and previously undiscovered components. the sphere of metagenomic analysis has extended from science to heavy industry, agriculture and healthcare. a large amount of metagenomic data is constantly being accumulated which leads to a demand in the means of efficient analysis. one of the common steps in metagenomic study is calculation of pairwise dissimilarity between the samples. beta-diversity is a quantitative measure of the differences in composition between two microbial communities. its value is calculated from features like taxonomic or functional composition, phylogenetic structure of the whole community, etc. a dissimilarity matrix composed of pairwise distances between all samples is used for further cluster analysis, classification and study of influence of the experimental factors. in large-scale studies involving tens and hundreds of metagenomes, critical requirements in beta-diversity analysis include high algorithm performance and low memory usage. for a long time, the standard technological approach for the evaluation of beta-diversity was based on the identification of species in metagenomic samples through 16s rrna gene sequencing. however, this method has inherent disadvantages including incompleteness of reference databases, presence of multiple copies of 16s rrna gene in the same genome, discrepancy between phylogenetic trees constructed using 16s rrna and the other genes and lack of information about the other genes and subsequently metabolic potential of the studied community. an alternative, more informative method is whole-genome sequencing generating millions of reads from the total dna of the genomes of all organisms inhabiting the environment. the identification of the organisms in the short-read wgs metagenome is commonly based on the alignment or de novo assembly. the alignment method is a comparison between sequences of obtained reads and sequences of reference genes or genomes, and has significant drawbacks such as high computational costs and incomplete databases. de novo assembly is usually a time-consuming task for such complex data as metagenomes that may contain many unknown or highly similar genomes of organisms with widely varying abundance. with a rapid increase in data output produced by sequencing technologies, efficient methods for genomic analysis based on k-mer composition analysis emerged. such algorithms work with k-mers obtained directly from metagenomic reads, without pre-mapping or assembly. in comparison with reference-based methods, the main advantages of k-mer based approaches are compressed representation of sequences and inclusion of the entire data volume into analysis. among these methods, the most simple and effective for exploratory analysis of large data sets is comparison of sequences by calculation of pairwise dissimilarity between them on the basis of k-mer spectrum - a normalized vector of frequencies of occurrences of each k-mer in the metagenome. the k-mer length is a key factor influencing specificity and efficiency of the comparison. for different intervals of k, the respective algorithms have been designed that target different specific tasks. for example, for very short k-mers only “rough” estimates are possible: sequence quality check, taxonomic separation of individual genomes or comparison of metagenomic communities with notably different composition. for k=15− the computational costs associated with the processing of the whole spectrum increase significantly. two approaches can be applied to reduce them. this problem can be solved in two ways. the first is to select a fraction of k-mers that describe the studied data in the most complete way. another way is to combine multiple k-mers into one feature using a certain principle. in the intermediate range, it is still computationally feasible to analyze the complete set of k-mers, and the k-mer length remains sufficiently specific for comparing distinct genomes. among the metagenomic k-mer methods, the most simple and effective for exploratory analysis of large data sets is comparison of sequences by calculation of pairwise distances between them on the basis of k-mer spectra. in this area, researchers are actively designing fast algorithms for calculating and assessing the k-mer spectrum. most studies are focused on examining clustering of the samples by one or more factors. since the prevalent approaches for assessment of beta-diversity today are reference-based, an important question is how their results correlate with the k-mer methods. in this study, we compared common reference-based methods with the k-mer approach. we explored how various characteristics of the data influence the results of k-mer spectra analysis and identified the advantages of k-mer analysis comparing with the reference-based approaches. to evaluate the applicability of k-mer-based dissimilarity, metagenome of human gut microbiota was selected, the study of which has great biomedical importance and perspective. although nowadays intestinal microbiota is one of the most studied among complex microbial communities, many of its components are still not fully identified. the application of our method revealed significant presence of one of such components - phage - that went undetected by reference-based methods. methods simulated metagenomes two set of “shotgun” gut metagenomes were simulated using metasim. the high-diversity set included metagenomes generated from the genomes of ten distantly related major bacterial species accounting for more than % of all reads in chinese group: akkermansia muciniphila atcc baa- alistipes shahii wal bacteroides vulgatus atcc bifidobacterium adolescentis atcc coprococcus sp. art55/ eubacterium eligens atcc faecalibacterium prausnitzii l2- lachnospiraceae bacterium 56faa, prevotella copri dsm and ruminococcus sp. 18p the simulation included the following steps. first, for each genome, mean and standard deviation of its relative abundance were estimated from the taxonomic composition of the chinese metagenomes. for each metagenome, ten abundance values were randomly generated under normal distribution with these parameters and the obtained values were normalized to million reads; a total of genera abundance vectors were obtained. the metagenomes were generated by mixing ten bacterial genomes at the obtained abundance levels and sampling short reads from the genomes using metasim with read length bp. also we performed sampling of these reads with errors. the low-diversity simulated group included metagenomes generated in a similar way from the genomes of ten closely related major bacterial species accounting for more than % of all reads in the hmp group: bacteroides vulgatus atcc bacteroides dorei 36/d bacteroides uniformis atcc bacteroides stercoris atcc bacteroides caccae atcc bacteroides ovatus, bacteroides xylanisolvens xb1a and bacteroides thetaiotaomicron vpi- bacterial proportions for these simulations are listed in additional file 1: table s for single nucleotide polymorphism simulations, the same ten reference genomes and abundance values as in the high-diversity dataset were used. two different models of snps introduction were used: “independent” and “phylogenetic”. in the “independent” snp model, “mutated” genomes were generated for each reference species by changing nucleotide letter at random positions independently with % substitution rate. thus, the average amount of snps between any two of the “mutated” genomes was ∼ %. in “phylogenetic” snp model, the procedure was performed in iterations for each reference genome: initialize with a single genome; iteration number = make a copy of each of the genomes available at the step. introduce snps to all genomes at random positions. increment iteration number. if the iteration number is greater than stop; else return to step b. after the iterations, 26= “mutated” genomes are obtained. in each model, the random “mutated” genomes of corresponding bacteria were used to generate metagenomes the same way as for high-diversity simulation above. real metagenomic datasets two “shotgun” gut metagenomic datasets were analyzed: metagenomes of healthy usa population and metagenomes of chinese population including healthy and type diabetes individuals. for each sample, the reads were filtered by quality using fastq quality filter script from fastx-toolkit. for each metagenome, million of high-quality reads was sampled using random _records script from. comparison of various sampling sizes showed that the selected size of subsampling does not significantly affect the results of the measures’ comparison. calculation of k-mer vectors and dissimilarity measures for each metagenome, k-mer spectrum was calculated using an ad hoc java program that processes fasta files read-wise by obtaining k-mer counts for each read and adding the counts to a global array. after processing all reads, the counts for reverse-complementary k-mers were summed and normalized to a sum of the length of the final feature vector did not exceed n=22k− for odd k and n=22k−1+2k− for even k because of reverse-complement k-mers. the obtained k-mer spectra were used to calculate pairwise dissimilarity via bray-curtis measure defined as: bc=1−2∑i=14kmin,mi)∑i=14k+mi) where m is the vector of k-mer frequencies normalized to a sum of per metagenome and x, y are two different metagenomes. bc = if the frequencies are equal for all k-mers between the metagenomes, and bc = if no common k-mers are present in the metagenomes. beta-diversity analysis using reference-based methods taxonomic profiling via mapping of metagenomes to a reference genome catalog and coverage analysis was performed as described previously, with the only difference: a non-redundant set of genomes of gut microbes was used. the final feature vector for each metagenome included relative abundance of microbial species was normalized to a sum of %. dissimilarity was calculated using these vectors both with bray-curtis measure and whole-genome adaptation of the weighted unifrac metric. functional profiling was performed as described previously to yield cog relative abundance vectors subsequently used for dissimilarity analysis using bray-curtis measure. an alternative method of taxonomic profiling employed metaphlan v. ); here mapping was performed using bowtie v. software, up to mismatches per read were allowed. all reference-based methods were summarized in table crassphage abundance analysis all crassphage genes were aligned to the reference gene catalog. for each gene, its relative abundance was estimated as a ratio of the total length of the reads mapped to this gene to the total length of the reads mapped to the reference gene catalog. phage relative abundance was determined as a sum of the relative abundance values of its genes. as an additional method of metagenomic classification, diamond aligner was used in combination with megan classifier. statistical analysis statistical analysis was implemented in r. the code is available at: http://download.ripcm.com/dubinkina_2015_suppl_data/ and https://github.com/zireae1/kmer_project/. ethical approval the sampling procedure was approved by the ethical committee for clinical research from the peking university shenzhen hospital, shenzhen second people’s hospital and medical research center of guangdong general hospital, and enrollments were approved by the institutional review boards of the two recruitment centers.  to compare the k-mer based metagenomic beta-diversity measure with traditional reference-based methods we conducted a series of computational experiments on simulated and real data. to validate the k-mer measure and find an optimal value of k, we simulated metagenomes from prevalent human gut bacterial genomes. then the method was applied for the analysis of a group of real human gut metagenomes sequenced in two large-scale projects: china population and hmp. comparison of beta-diversity measures for simulated metagenomes there is considerable variation between k-mer spectra for genomes of distinct bacterial species due to the differences in the gene content, amino acid coding preferences, etc.. supposedly, in the case of a metagenome including sufficiently covered genomes of multiple species, one should observe significant accordance between k-mer spectrum of the metagenome and its taxonomic composition. to verify this hypothesis, we simulated several datasets with different degrees of community richness and applied the bray-curtis measure for both taxonomic and k-mer profiles to compare the two respective dissimilarity matrices. simulation 1: high-diversity communities the first synthetic dataset included metagenomes generated by randomly sampling “reads” from ten genomes of common phylogenetically distant human gut bacteria. comparison of the two approaches showed that, as k increases, so does the correlation value between the two dissimilarity matrices based on k-mers and taxonomic composition. with high values of k, the two matrices are highly similar. simulation 2: sequencing errors and snps besides the considerable number of species within microbiota, the other factors contributing to the diversity of metagenomic k-mers are presence of mutations and sequencing errors in reads. therefore, we conducted two experiments by introducing artificial snps into genomes and, separately, random single-nucleotide changes into the reads in order to explore their influence on the correlation of beta-diversity estimates using k-mer and taxonomic methods. datasets from simulation were used here. for sequencing errors modeling, the reads for each metagenome were simulated with per-nucleotide substitution rate of %. introduction of such “errors” did not lead to a significant change in correlation between the two methods. for snp modeling, bacterial genomes with % of randomly introduced single-nucleotide substitutions were used to generate simulated metagenomes with the same abundance proportions as in simulation we employed two different models of snps introduction - “independent” and “phylogenetic”. with the former simulation being more straight-forward, “phylogenetic” approach was developed to model the accumulation of mutations in bacterial species during evolutionary process. the results of simulations showed that, independent of the model choice, in general snps had a minor effect on k-mer spectra comparable to the effect of simulated sequencing errors: correlation between the k-mer and taxonomic methods decreased from to for “independent” model and for “phylogenetic” model. noteworthy, introduction of snps had a more pronounced effect on metagenomes with highly similar taxonomic composition. this was particularly marked when the snp rate was increased to %. simulation 3: low-diversity communities the second synthetic dataset included metagenomes generated by randomly sampling “reads” from ten genomes of common phylogenetically close human gut bacteria - belonging to the same genus - bacteroides. the correlation between the methods was found to be lower for such homogeneous community than for a heterogeneous one. the correlation value tends to increase with k but does not achieve the level of simulation. it suggests that higher values of k should be used to increase accuracy; however, the size of the feature vector increases as 4k, hence the computational time quickly becomes unacceptable. to select the optimal k value, we evaluated the correlation between k-mer and taxonomic dissimilarity matrices together with the computational time of k-mer matrix generation for k=5− using both high- and low-diversity simulated datasets. as the results in both simulations showed, with k = the dissimilarity matrices are highly correlated while the computational time is still acceptable took about seconds to calculate the spectra) further statistical analysis - calculation of dissimilarity matrix - took about - minutes, see additional file 5: figure s3). at the same time, it is the highest value practically acceptable in terms of memory usage: for k= the spectra occupied ∼ gb of memory in r environment, but for k = - as much as gb. considering these observations we selected k= for further analyses. comparison of beta-diversity measures on real human gut metagenomes after testing the method on simulated data, we examined two real human gut datasets from large-scale metagenomic projects: china and hmp, with the former cohort representing more diverse microbial community structures than the latter. using this data, the pairwise dissimilarity matrix obtained by the k-mer approach with bray-curtis measure was compared with the dissimilarity matrices obtained by each of the four methods based on taxonomic and functional reference. to visualize the distributions of beta-diversity values, we applied two types of scatter plots. the first type is a basic principle coordinate analysis plot constructed using a single dissimilarity measure, with dots representing distinct metagenomes. on the second type of plot, two dissimilarity measures are compared: each triangle corresponds to a pair of metagenomes, one measure is plotted against the other. samples from the two studies tended to cluster separately by functional, as well as by k-mer composition, but not by taxonomic composition. therefore, the two cohorts were further analyzed separately. another interesting fact was that of the outliers present on k-mer scatter plot were also on the periphery of cog scatter plot but not of the taxonomic scatter plot. these samples were examined in details. fig. variation of metagenomes using different dissimilarity measures. pcoa plots for different dissimilarity measures: a bc kmer, b bc cog, c wg unifrac, d bc tax, e bc metaphlan. three samples-outliers are marked with asterisks. f heatmap of spearman correlation coefficient between dissimilarity matrices obtained using different measures fig. comparison of pairwise difference measures obtained by k-mer and reference-based methods. for each plot, y-axis represents k-mer distance, x-axis - distance by one of the reference-based methods. distribution of dissimilarity measures is shown for a bc kmer for all reads and bc tax; b bc kmer for all reads and bc cog; c bc kmer for reads mapped to the catalog of genomes and bc tax; d bc kmer for reads mapped to the catalog of genes and bc cog comparison of the five beta-diversity measures showed that the k-mer measure has a significant similarity with each of the reference-based ones. the closest was the measure based on cog composition. for the chinese group, the correlation values tended to be higher than for the hmp group in all comparisons. the phylogeny-aware metric wg unifrac was among the most dissimilar. influence of reads mappability to assess the contribution of the unmapped reads to the results, k-mer spectra were also computed only using the reads that successfully mapped to the corresponding catalog. this analysis led to interesting results. first, we observed an equalization of bc tax org vs. bc kmer correlation between the two cohorts. therefore, fraction of unmapped reads appears to be one of the major factors contributing to the difference between the cohorts. this parameter is dependent on the representability of the reference catalog and quality of sequencing run. second, we assessed the shift of each outlier in the direction of the central cloud of points. quantitatively, for each outlier we calculated the bc kmer difference value: the difference between its bc kmer value and the linearly interpolated middle of the cloud obtained for the same reference-based value. for comparison with bc tax, the bc kmer difference decreased significantly - from ± to ±. for comparison with bc cog, the bc kmer difference changed slightly: from ± to ±. correspondingly, a group of pairs-outliers mentioned above moved into the central cloud of points in the bc tax org vs. bc kmer comparison, but did not change their visual location in bc cog vs. bc kmer comparison. this observation is in agreement with the fact that the gene reference catalog is more complete than the genome reference catalog and the percentage of mapping to the gene catalog is higher. presumably, the presence of pairs-outliers could be caused by k-mers from certain dominant sequences that are present in the reference base of genes but not genomes. we investigated these outliers in details. investigation of samples-outliers the total human gut metagenome is a phylogenetically diverse structure including not only the sequences of bacterial genomes but also ones from bacterial mobilome, fungi, protozoa, traces of dna of dietary origin, host. our reference genome catalog partly accounts for such non-microbial components by including the genomes of several common intestinal eukaryotes - clinically relevant yeasts candida and protozoan blastocystis. however, many sequences are not present in our genome catalog, particularly viral genomes. therefore, in our analysis the potential reads of viral origin would not contribute to the taxonomic difference but would change the k-mer spectrum. recently, a new bacteriophage was discovered - crassphage - shown to be a sole major dominant of the human gut viriome. moreover, its presence was estimated for the hmp metagenomes analyzed in our work: crassphage genome amounts for up to % of the reads for this group. obviously, such a prevalent genome should have a significant influence on k-mer spectra and thus on our comparison of the beta-diversity measures. basing on the available data on the abundance of crassphage in hmp samples, the cohort was split into two groups - with high phage abundance and with low phage abundance. the whole group of extreme outliers was found to consist of the pairs where at least one of the samples belonged to the phage-enriched group. moreover, the outliers on fig. 1a were also found to be the samples with high fraction of crassphage reads. fig. analysis of samples-outliers. a distribution of pairwise dissimilarity obtained using k-mer and taxonomic composition for hmp cohort. different colors indicate groups of dissimilarities for: all hmp pairs, pairs-outliers - where at least one of the samples belonged to the phage-enriched group; cp-filtered pairs - extreme outliers after removal of k-mers from reads mapped to crassphage genome; b composition of sample srs according to the combined results from two analyses as the extreme outliers were found to be generated by pairs including at least one of the two samples, these samples were analyzed in detail. the reads that did not map to the genome catalog were subject to metagenomic classification using an alternative method - using diamond alignment and megan classifier algorithms. as a result, additionally % and % of the reads were identified as crassphage. to further confirm the contribution of high phage fraction to formation of outliers, we subtracted the k-mers of the crassphage reads from k-mer spectra of the samples. indeed, such operation significantly decreased the k-mer-based dissimilarity for the respective pairs. discussion here we have developed an algorithm for assessing pairwise dissimilarity of “shotgun” metagenomes basing on k-mer spectrum and compared it with commonly used reference-based approaches. the comparison was performed using various measures on a set of simulated metagenomes as well as on real metagenomes from two large-scale human gut microbiota studies. for simulated metagenomes, we showed that k = is an optimal value in terms of balancing between the resolution of the method and computational time. this value of k performed well for both high- and low-diversity simulated metagenomes; however, for low-diversity simulations the dissimilarity matrices based on k-mer method and taxonomic composition were less correlated. this fact was likely due to the decreased diversity of k-mers and thus reduced differentiating resolution. for real gut metagenomes with complex community structure, the k-mer approach allows to delineate the samples with a wide range of functional composition, as demonstrated on two international cohorts. on the other hand, we observed that k-mers are less correlated with taxonomic composition than with functional one. we speculate that this difference could be associated with significant subspecies-level genomic diversity of gut microbes: a recent analysis of publicly available metagenomic data showed that the average gene variation between individuals across abundant species was as high as ± %. the k-mer frequencies as well as the gene relative abundance features are sensitive to gene content variation, while in the case of species relative abundance features it would be ignored. besides gene presence/absence, another common form of genomic variability is snps. we attempted to model their influence on k-mer beta-diversity. theoretically, introduction of snps would lead to change in frequencies of k-mers and thus deteriorate correlation between k-mer and reference-based dissimilarity. in our simulations, when % snps were introduced to simulated datasets, the correlation between the methods dropped slightly, independent of whether the evolutionary character of snp accumulation was considered during modeling or not. however, for real metagenomes the correlation between the methods was lower. these results suggest the existence of other real-life effects having stronger influence on the correlation than snps. a major advantage of the k-mer approach is that it exploits the totality of the reads - unlike reference-based methods that inherently discard the reads that failed to map to the reference catalog - and thus the information contained within them. such a feature promotes the application of the k-mer approach as a tool for assessing the representability of the reference set for given metagenomes. currently representative sets of reference genomes are available for microbial communities of few environments. however, recent discoveries imply that the so called reference genomes do not capture a wide intra-species level variation even for this extensively examined community. the situation is even worse for less popular environments - like marine ecosystems or human skin: reference catalogs for their microbiota are considerably less complete, thus rendering beta-diversity assessment difficult. we propose to assess the representability of a reference genome catalog via examining the k-mer content of the metagenomic reads mapped to it. on the analyzed gut metagenomic data, we observed that k-mer spectra of the mapped reads produced dissimilarity profiles that had higher correlation with those obtained with taxonomic composition than the k-mer spectra of the whole set of reads. however, lower correlation between the two methods observed for some pairs of samples suggested the presence of dominant genomic sequences not included in the reference catalog. detailed analysis showed that these outliers corresponded to the hmp samples enriched in crassphage, a recently discovered gut bacteriophage; the genome of this phage was not included in the respective reference catalog. subtraction of the crassphage k-mers moved the outliers towards the main cloud of points but not into it completely. presumably, such incomplete compensation can be linked to high level of genomic variability inherent to gut phages: originally the consensus crassphage genome was obtained by combined assembly of metagenomes from individuals not included in our groups, so its sequence in the latter might be quite distant than all the crassphage-related k-mers in our analysis. additionally, over % of the reads remain unidentified by two different methods they can correspond to genome contributing to formation of outliers. considering the gene catalog, dedicated analysis of the reads mapped to reference genes did not lead to shift of outliers. first, a likely reason for this is that the crassphage sequences were included in the catalog: search for crassphage genes in the reference gene catalog identified highly similar hits for of the phage genes that were detected in at least one metagenome. second, the gut microbial gene catalog was originally constructed from the contigs assembled from total dna reads and is known to contain not only the bacterial genes, but viral and eukaryotic, too. interestingly, our results also imply that the chinese cohort lacks metagenomes with such high prevalence of this phage, provoking speculations on world-wide phage phylogeography. while no clinical associations for crassphage have been described to date, omission of phage components could be a significant miss in biomedical studies of microbiota. there is a growing understanding that gut phages play an important role in the ecology of “phage-gut microbiota-human” system and include potential biomarkers; they are able to transfer clinically important bacterial genes - e.g. antibiotic resistance and pathogenicity determinants. application of our reference-free k-mer approach can facilitate early detection of such sequences in biomedical diagnostics data and discovery of novel biomarkers. our approach is not only applicable to metagenomes from an arbitrary environment, but is indispensable for dissimilarity and cluster analysis of communities with poorly described components. the approach allows to detect a major presence of an unknown organism and/or virus in a metagenome. we suggest that the approach should be introduced as a necessary method of “shotgun” metagenome composition analysis complementary to reference mapping in order to avoid biases associated with unrepresentative reference database. although we did not find evidence for outliers caused by technical issues in the examined datasets, the approach can also be used for primary detection of metagenomes with abnormal composition caused by high abundance of host dna, dna of dietary origin and technical artifacts. finally, comparison of the metagenomes basing on k-mer spectrum provides more information than mapping to reference sequence catalogs. essentially, k-mer analysis is a feature extraction procedure applied to metagenomic reads. the produced set of features is several orders of magnitude larger than one yielded in reference-based approaches. therefore, it provides higher discriminative resolution that opens a promising opportunity for developing a new generation of methods for metagenomic analysis, and our method makes a step towards understanding of how to explore such high-dimensional feature space efficiently.  analysis of k-mer spectra for both simulated and real “shotgun” metagenomes showed that this method allows quick assessment of the pairwise dissimilarity of such datasets. simulations show that the method is robust to variability introduced by sequencing errors and genomic mutations. the obtained dissimilarity matrix can be used not only for cluster analysis and classification purposes, but also for early detection of major unknown components and quality control of reference-based approaches. it is recommended that the method should be included as a complementary step in high-throughput computational pipelines for metagenomic data analysis. additional files additional file supplementary tables. supplementary table s bacterial abundances in simulation. supplementary table s bacterial abundances in simulation. supplementary table s list of genomes in taxonomic catalog for human gut. supplementary table s taxonomic composition for real dataset. supplementary table s taxonomic composition for real dataset. supplementary table s functional composition for real dataset. supplementary table s taxonomic composition for real dataset by metaphlan. supplementary table s mapped read counts and percentage of mapping on taxonomic and functional catalog and phage genome. additional file figure s correlation between k-mer-based dissimilarity matrix obtained using the original full metagenomic readsets and their subsampled versions. the same amount of reads was repeatedly sampled from randomly selected metagenomes. additional file figure s comparison of pairwise difference measures obtained by k-mer and taxonomic composition. for each pair y-axis represents k-mer dissimilarity, x-axis represents dissimilarity by taxonomic composition. k= a) using the original bacterial genomes; b) using the bacterial genomes with % of snps. additional file figure s comparison of pairwise dissimilarity measures obtained by k-mer and taxonomic composition for simulated for high- and low-diversity metagenomes. as seen, satisfactory correlation of k-mers with taxonomic composition can be obtained only at relatively high values of k. additional file figure s correlation between k-mer and taxonomic composition dissimilarity matrices, as well as k-mer dissimilarity matrix computation time with varying values of k. all computations were performed on a compute node with cpu opteron ghz and gb ram. additional file figure s comparison of dissimilarity measures obtained by k-mer and reference-based methods: bc metaphlan genus, bc metaphlan org and wg unifrac. for each plot, y-axis represents k-mer dissimilarity, x-axis - dissimilarity using one of reference-based methods. abbreviations 16s rrna16s ribosomal ribonucleic acid bcbray-curtis dissimilarity measure bc kmerbc based on k-mer spectrum bc taxbc based on taxonomic composition, obtained by mapping to catalog of genomes bc metaphlanbc based on taxonomic composition, obtained by mapping to clade-specific marker genes bc cogbc based on functional composition, obtained by mapping to catalog of cog genes dnadeoxyribonucleic acid snpsingle nucleotide polymorphism wgswhole genome sequencing, “shotgun” wg unifracwhole genome adaptation of unifrac measure competing interests the authors declare that they have no competing interests. authors’ contributions dga and avt initiated the project and designed the study. viu wrote the programs for k-mer counting. vbd wrote the programs for dissimilarity matrix calculation, generated the simulation datasets and analyzed the data. dsi performed statistical modeling of snps effect. dga and avt guided the comparison of dissimilarity measures obtained by k-mer and reference-based methods. vbd and avt wrote the manuscript under the supervision of dga. all authors read and approved the final manuscript. 
the k-mer spectrum-based measure was found to behave similarly to one based on mapping to a reference gene catalog, but different from one using a genome catalog. here we evaluated a metagenomic pairwise dissimilarity measure based on short k-mer spectrum using the example of human gut microbiota, a biomedically significant object of study. dissimilarity measure based on k-mer spectrum provides a wider perspective in comparison with the ones based on the alignment against reference sequence sets. we developed a method for calculating pairwise dissimilarity of “shotgun” metagenomes based on short k-mer spectra. our approach allows rapid estimation of pairwise dissimilarity between metagenomes. a rapidly increasing flow of genomic data requires the development of efficient methods for obtaining its compact representation. it helps not to miss possible outstanding features of metagenomic composition, particularly related to the presence of an unknown bacteria, virus or eukaryote, as well as to technical artifacts at the early stages of bioinformatic analysis. this difference turned out to be associated with a significant presence of viral reads in a number of metagenomes. simulations showed limited impact of bacterial genetic variability as well as sequencing errors on k-mer spectra.
 a rapidly increasing flow of genomic data requires the development of efficient methods for obtaining its compact representation. feature extraction facilitates classification, clustering and model analysis for testing and refining biological hypotheses. “shotgun” metagenome is an analytically challenging type of genomic data - containing sequences of all genes from the totality of a complex microbial community. recently, researchers started to analyze metagenomes using reference-free methods based on the analysis of oligonucleotides frequency spectrum previously applied to isolated genomes. however, little is known about their correlation with the existing approaches for metagenomic feature extraction, as well as the limits of applicability. here we evaluated a metagenomic pairwise dissimilarity measure based on short k-mer spectrum using the example of human gut microbiota, a biomedically significant object of study.  we developed a method for calculating pairwise dissimilarity of “shotgun” metagenomes based on short k-mer spectra. the method was validated on simulated metagenomes and further applied to a large collection of human gut metagenomes from the populations of the world. the k-mer spectrum-based measure was found to behave similarly to one based on mapping to a reference gene catalog, but different from one using a genome catalog. this difference turned out to be associated with a significant presence of viral reads in a number of metagenomes. simulations showed limited impact of bacterial genetic variability as well as sequencing errors on k-mer spectra. specific differences between the datasets from individual populations were identified.  our approach allows rapid estimation of pairwise dissimilarity between metagenomes. though we applied this technique to gut microbiota, it should be useful for arbitrary metagenomes, even metagenomes with novel microbiota. dissimilarity measure based on k-mer spectrum provides a wider perspective in comparison with the ones based on the alignment against reference sequence sets. it helps not to miss possible outstanding features of metagenomic composition, particularly related to the presence of an unknown bacteria, virus or eukaryote, as well as to technical artifacts at the early stages of bioinformatic analysis. our method is complementary to reference-based approaches and can be easily integrated into metagenomic analysis pipelines. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords k-mern-graml-tuplesequence signaturegut metagenomephagereference-free metagenomic analysisgenomic variabilitymapping biashttp://dx.doi.org/ /501100003443ministry of education and science of the russian federation rfmefi60414x0119alexeev dmitry g. government of russian federationgrant 074-u01ulyantsev vladimir i. russian scientific foundationproject #15-14-00066tyakht alexander v. http://dx.doi.org/ /501100003443ministry of education and science of the russian federation rfmefi60414x0119dubinkina veronika b. russian scientific foundationproject #15-14-00066ischenko dmitry s. issue-copyright-statement© the author 2016 