| Sample | Collection date | Days | Seq. length | Seq. GC (%) | Seq. gaps | Seq. ambiguities | Mapped reads | % coverage (depth ≥ 4) | Mean depth (depth ≥ 4) | Mean baseQ (depth ≥ 4) | Mean mapQ (depth ≥ 4) | Mean read length (bp) | Read length SD (bp) | Median read length (bp) | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | MIP00022 | 1995-11-01 | 0 | 4411532 | 62.2 | 21505 | 203430 | 5725105 | 98.5 | 94.4 | 36.6 | 59.7 | 72.8 | 50.3 | 51.0 |
| 2 | MIP00023 | 1996-07-01 | 243 | 4411532 | 62.2 | 22777 | 203429 | 5470629 | 98.5 | 87.8 | 36.5 | 59.7 | 70.9 | 47.8 | 51.0 |
| 3 | MIP00024 | 1996-11-01 | 366 | 4411532 | 62.1 | 28904 | 203481 | 5500055 | 98.4 | 84.8 | 36.6 | 59.7 | 68.1 | 42.9 | 51.0 |
| 4 | MIP00025 | 1997-01-01 | 427 | 4411532 | 62.2 | 22593 | 203411 | 5732091 | 98.5 | 93.3 | 36.5 | 59.7 | 72.0 | 49.8 | 51.0 |
| 5 | MIP00026 | 1998-03-01 | 851 | 4411532 | 62.2 | 22872 | 203390 | 5714197 | 98.5 | 95.0 | 36.5 | 59.7 | 73.5 | 51.9 | 51.0 |
| 6 | MIP00027 | 1998-06-01 | 943 | 4411532 | 62.2 | 23082 | 203420 | 5851858 | 98.4 | 93.8 | 36.5 | 59.7 | 70.8 | 47.8 | 51.0 |
| 7 | MIP00028 | 1998-09-01 | 1035 | 4411532 | 62.2 | 22834 | 203423 | 5680594 | 98.5 | 92.6 | 36.6 | 59.7 | 72.0 | 49.6 | 51.0 |
| 8 | MIP00029 | 1999-02-01 | 1188 | 4411532 | 62.2 | 20731 | 203442 | 5928246 | 98.5 | 99.0 | 36.6 | 59.7 | 73.8 | 51.0 | 51.0 |
| 9 | MIP00030 | 1999-05-01 | 1277 | 4411532 | 62.2 | 22256 | 203432 | 5699653 | 98.5 | 94.8 | 36.5 | 59.7 | 73.5 | 51.6 | 51.0 |
Summary
Study description
- Organism: M. tuberculosis
- Reference: MTB_anc
- Sequence alignment: none
- Variant calling: freebayes (against dataset ancestor)
Sample metadata is summarised in Table 1. Figure 1 shows the structure of the dataset.
List of nucleotide variants
The sequence of the most recent common ancestor (MRCA) of the target samples has been built using ancestral sequence reconstruction (ASR) and compared to the mapping reference (Table 2). These variants represent fixed differences already present before the processes under study began.
To distinguish pre-existing variation from changes arising during the study period, nucleotide variants identified in the ancestor relative to the reference are reported separately from those identified in each sample relative to the ancestor, listed in Table 3.
Variants overview
A total of 66 different nucleotide variants were identified in the target samples, of which 11 were ‘Frameshift’, 1 were ‘In-frame indel’, 8 were ‘Intergenic’, 36 were ‘Non-synonymous’, 1 were ‘Other’, and 9 were ‘Synonymous’ (Figure 2).
Allele frequency trajectories
The correlation of the allele frequency of each nucleotide variant with time since initial sampling has been calculated using Pearson’s correlation method (Figure 3). The graph shows which variants pass the correlation threshold (0.7), the p-value threshold (0.05), or both.
Selected correlated nucleotide variants are described in detail in Figure 4 and Figure 5. The visualizations are separated to differentiate variants with larger allele frequency amplitude (differences between their maximum and minimum alternate frequencies). In both cases, the plots show the alternative allele frequency across sampling days. Plot colors correspond to Figure 3, indicating whether each variant passes the correlation threshold, the p-value threshold, or both.
Temporal signal
To estimate the evolutionary rate based on all alleles (Figure 8), including those at low frequencies, a BIONJ (modified neighbor-joining) tree has been constructed from pairwise allele frequency-weighted distances between study samples using afwdist (Figure 6). To estimate the evolutionary rate from majority alleles (Figure 9), a maximum likelihood (ML) tree has been inferred under a GTR+F+I+G4 substitution model using IQ-TREE (Figure 7). Results are summarized in Table 6.
| Method | Alleles | Sites | Rate | Rate per site | R² | P |
|---|---|---|---|---|---|---|
| BIONJ | All (weighted) | 4111304 | 4.05 [3.31, 4.78] | 9.84e-07 [8.06e-07, 1.16e-06] | 0.961 | < 0.001 |
| ML | Majority | 4111304 | 2.37 [0.443, 4.29] | 5.76e-07 [1.08e-07, 1.04e-06] | 0.547 | 0.0227 |
Evolutionary metrics
To track signatures of selection, the allele frequency-weighted rate of change per synonymous site (dS analog) and per non-synonymous site (dN analog) have been calculated for each sample relative to the provided reference sequence, together with their ratio ω (Figure 10). These metrics quantify frequency-weighted divergence across all alleles and help identify selective pressures. Frequency-weighted πN and πS of biallelic variants are reported as a metric of standing diversity (Figure 11). Rate denominators are calculated using the Nei–Gojobori (1986) method. Frequency-weighted synonymous and non-synonymous counts are reported as raw values, relative to the first sample, and as a ratio (Figure 12).
The site frequency spectrum (SFS) has been calculated to help interpret the estimates above. Because alleles may occur at varying frequencies across samples, the SFS has been calculated using allele frequency bins, stratified by synonymous single-nucleotide variants (SNV), non-synonymous SNV, and all variants combined. To account for varying allele frequency trajectories across the dataset, the SFS is reported using two approaches: “exclusive”, which captures the overall distribution of within-host allele frequencies across the dataset (Figure 13), and “inclusive”, which identifies variants with globally stable frequency behavior (Figure 14).
Variants connectivity
Hierarchical clustering of pairwise AF correlations
To detect possible interactions between mutations, pairwise Pearson’s correlation between allele frequencies have been calculated (Figure 15). The heatmap is interactive and allows zooming into specific regions.
Network communities of pairwise AF correlations
Studying how variants are connected allows us to identify behavioral patterns and relate annotated variants to those lacking annotation. In this analysis, all variants are connected to one another (Figure 16), rather than being examined in pairs as in the previous approach (Figure 15).
The network representation (Figure 16) of pairwise correlations between variant allele frequencies as edge weights. Graph communities have been identified using the Louvain algorithm. Each color represents a community with its corresponding number.
Custom annotation of variants
An overview of the presence of variants with DRUG annotation is shown in Figure 17. The variation in frequency of these variants across different samples in time is shown in Figure 18. Ancestor variants have been identified by comparing to the reference sequence. Sample variants have been identified by comparing each sample to the provided reference sequence.
These represent genetic changes that arose after that ancestor. Any allele already present in the ancestor that has not been displaced is expected to persist in the samples at detectable frequency.