The statistical power of k-mer based aggregative statistics for alignment-free detection of horizontal gene transfer

Abstract

Alignment-based database search and sequence comparison are commonly used to detect horizontal gene transfer (HGT). However, with the rapid increase of sequencing depth, hundreds of thousands of contigs are routinely assembled from metagenomics studies, which challenges alignment-based HGT analysis by overwhelming the known reference sequences. Detecting HGT by k-mer statistics thus becomes an attractive alternative. These alignment-free statistics have been demonstrated in high performance and efficiency in whole-genome and transcriptome comparisons. To adapt k-mer statistics for HGT detection, we developed two aggregative statistics ${\text{T}}_{sum}^{\text{S}}$ and ${\text{T}}_{sum}^{\text{*}}$ , which subsample metagenome contigs by their representative regions, and summarize the regional ${\text{D}}_2^S$ and ${\text{D}}_2^{\text{*}}$ metrics by their upper bounds. We systematically studied the aggregative statistics' power at different k-mer size using simulations. Our analysis showed that, in general, the power of ${\text{T}}_{sum}^{\text{S}}$ and ${\text{T}}_{sum}^{\text{*}}$ increases with sequencing coverage, and reaches a maximum power >80% at k = 6, with 5% Type-I error and the coverage ratio >0.2x. The statistical power of ${\text{T}}_{sum}^{\text{S}}$ and ${\text{T}}_{sum}^{\text{*}}$ was evaluated with realistic simulations of HGT mechanism, sequencing depth, read length, and base error. We expect these statistics to be useful distance metrics for identifying HGT in metagenomic studies.

Keywords: Alignment-free sequence comparison; Horizontal gene transfer; Statistical power; k-mer.

Fig. 1

The ${\text{T}}_{sum}$ resampling scheme. (On Seq $\text{X}$ and Seq $\text{Y}$, F₁ to F_N are the subsampled fragments, and G is gap. ${\text{X}}_{\text{i}}^{\text{S}}$ is the maxima of ${\text{D}}_2^{\text{S}}$ between the ith subsampled fragment of $\text{X}$ and all subsampled fragments of $\text{Y}$. The same for ${\text{Y}}_i^{\text{S}}$ is the maxima of ${\text{D}}_2^{\text{S}}$ between the ith subsampled fragment of $\text{Y}$ and all subsampled fragments of $\text{X}$.)

Fig. 2

Simulating the foreground sequence pairs using the Horizontal Gene Transfer (HGT) procedure with motif length L = 5.

Fig. 3

Statistical power of ${\text{T}}_{sum}^S$ and ${\text{T}}_{sum}^{\text{*}}$ when the coverage rate R = 25%, 50%, 75%. (Two statistics ${\text{T}}_{\text{sum}}^{\text{S}}$ and ${\text{T}}_{\text{sum}}^{\text{*}}$ are used in this figure. Full length is the length of whole sequence, R is subsampling coverage rate, k is k-mer's length and L is the length of motif. The probability of four bases is P(A) = P(C)=P(G) = P(T) = 1/4 in e.i.i.d model, and P(A) = P(T) = 1/6, P(G) = P(C) = 1/3 in n.i.i.d (gc-rich) model.).

Fig. 4

Statistical power of ${\text{T}}_{sum}^S$ and ${\text{T}}_{sum}^{\text{*}}$ when k = 4, 5, 6, 7, 8. (Two statistics ${\text{T}}_{\text{sum}}^{\text{S}}$ and ${\text{T}}_{\text{sum}}^{\text{*}}$ are used in this figure. E is the entire genome segment length, R is subsampling coverage rate, k is k-mer's length and L is the length of motif. The probability of four bases is P(A) = P(C)=P(G) = P(T) = 1/4 in e.i.i.d model, and P(A) = P(T) = 1/6, P(G) = P(C) = 1/3 in n.i.i.d (gc-rich) model.)

Fig. 5

Statistical power of ${\text{T}}_{sum}^S$ and ${\text{T}}_{sum}^{\text{*}}$ when E = 400, 600, 800, 1000. (Two statistics ${\text{T}}_{\text{sum}}^{\text{S}}$ and ${\text{T}}_{\text{sum}}^{\text{*}}$ are used in this figure. E is the entire genome segment length, R is subsampling coverage rate, k is k-mer's length and L is the length of motif. The probability of four bases is P(A) = P(C)=P(G) = P(T) = 1/4 in e.i.i.d model, and P(A) = P(T) = 1/6, P(G) = P(C) = 1/3 in n.i.i.d (gc-rich) model.)