Score distributions of gapped multiple sequence alignments down to the low-probability tail

Assessing the significance of alignment scores of optimally aligned DNA or amino acid sequences can be achieved via the knowledge of the score distribution of random sequences. But this requires obtaining the distribution in the biologically relevant high-scoring region, where the probabilities are exponentially small. For gapless local alignments of infinitely long sequences this distribution is known analytically to follow a Gumbel distribution. Distributions for gapped local alignments and global alignments of finite lengths can only be obtained numerically. To obtain result for the small-probability region, specific statistical mechanics-based rare-event algorithms can be applied. In previous studies, this was achieved for pairwise alignments. They showed that, contrary to results from previous simple sampling studies, strong deviations from the Gumbel distribution occur in case of finite sequence lengths. Here we extend the studies to multiple sequence alignments with gaps, which are much more relevant for practical applications in molecular biology. We study the distributions of scores over a large range of the support, reaching probabilities as small as ${10}^{-160}$, for global and local (sum-of-pair scores) multiple alignments. We find that even after suitable rescaling, eliminating the sequence-length dependence, the distributions for multiple alignment differ from the pairwise alignment case. Furthermore, we also show that the previously discussed Gaussian correction to the Gumbel distribution needs to be refined, also for the case of pairwise alignments.

Figure 1

Equilibration times for different temperatures: First 10 000 Monte Carlo sweeps in the parallel tempering time series averaged over 30 different realizations for random (bottom branch) and high-scoring (top branch) initial conditions, respectively. System: $N_{\mathrm{seq}}=3$ sequences of length $L=40$, local alignment. A selection of the used temperatures is shown. At lower temperatures it takes more time for equilibration.

Figure 2

The distribution for local alignments of $N_{\mathrm{seq}}=3$ sequences of length $L=40$ as obtained by the rare-event simulation. Fits of the Gumbel distribution to the whole distribution as well as constrained to the score range $[S_l=21,S_u=115]$, i.e., $P\left(S\right)\le {10}^{-10}$, are shown. Also shown is the better-suited fit of the Gumbel distribution with a Gaussian correction. The inset shows the high-probability region.

Figure 3

The score distributions for local multiple alignments for several sequence lengths and the fitted Gumbel distributions with Gaussian correction. The inset shows the high-probability region.

Figure 4

The fit parameters found for different window sizes $(S_{\text{min}},S_{\text{end}})$. Several sequence lengths $L$ for $N_{\mathrm{seq}}=3$ are shown.

Figure 5

The fit parameters found for different window sizes $(S_{\text{min}},S_{\text{end}})$. Several sequence lengths $L$ for $N_{\mathrm{seq}}=2$ are shown.

Figure 6

The distributions for various $L$ and $N_{\mathrm{seq}}$ with rescaled scores and probabilities. The distributions coincide for low probabilities and identical number of sequences $N_{\mathrm{seq}}$ but differ among pairwise and multiple alignments. Inset: High probability region in which the distributions deviate from each other (even for same $N_{\mathrm{seq}}$).

Figure 7

The alignment score distribution as obtained with the rare-event simulations for the global alignment of $N_{\mathrm{seq}}=2$ sequences of length $L=100$, here with gap costs $(\alpha =7,\beta =-1)$ for direct comparison with Ref. [13]. Shown are the fits of the $\mathrm{\Gamma }$ distribution with and without Gaussian correction as well as the fit of the $\mathrm{\Gamma }$ distribution to the high-probability region ($P\left(S\right)>{10}^{-3}, S=[R_{\min }=-40:R_{\max }=23]$) only. Inset: high-probability region.

Figure 8

Distributions obtained for the global alignment of $N_{\mathrm{seq}}=3$ sequences of length $L=40$ with the exact algorithm and the heuristics. The lines show the fit of the uncorrected $\mathrm{\Gamma }$ function against the data. The inset shows the high-probability region, in which both distributions differ significantly from each other and the fit for exact alignments deviates significantly from the data.

Figure 9

The parameter ${\lambda }_2$ as found for the distributions of global alignment for different sequence lengths $L$ and different number of sequences $N_{\mathrm{seq}}$. Drawn lines are guides for the eye.

Figure 10

The parameter $\gamma$ as found for the distributions for global alignment of different sequence lengths $L$ and different number of sequences $N_{\mathrm{seq}}$.

Figure 11

Distributions rescaled with Eq. (11) for different number of sequences $N_{\mathrm{seq}}$ and different sequence lengths $L$. The rescaled distributions coincide for low-probabilities and identical $N_{\mathrm{seq}}$ but differing $L$.

Figure 12

Distributions obtained for global alignment of $N_{\mathrm{seq}}=2$ and $N_{\mathrm{seq}}=3$ sequences of length $L=40$ (with the exact algorithm as well as with the heuristics for $N_{\mathrm{seq}}=3$). The $log\left(P\right)$ values are rescaled with Eq. (11) and the scores with $S/S_{\text{max}}$. The inherit difference between $N_{\mathrm{seq}}=2$ and $N_{\mathrm{seq}}=3$ is obviously not based the use of an heuristics. The inset show the high-probability region, in which all three obtained distributions differ. In the case of $N_{\mathrm{seq}}=3$ this difference is small compared to the exact algorithm and is due to the bad performance of the heuristics in the low-scoring region.

Figure 13

Obtained distributions for local as well as global alignment of $N_{\mathrm{seq}}=3$ sequences of length $L=60$. The inset shows the high-scoring region, where both distributions approach each other.