PyDamage: automated ancient damage identification and estimation for contigs in ancient DNA de novo assembly

Simulated sequencing data

In order to evaluate the performance of PyDamage with respect to the GC content of the assembled genome, the sequencing depth along the genome, the amount of observed aDNA damage on the DNA fragments, and the mean length of these DNA fragments, we simulated short-read sequencing data using gargammel (Renaud et al., 2017) varying these four parameters. We chose three microbiome-associated microbial taxa with low(Methanobrevibacter smithii, 31%), medium (Tannerella forsythia, 47%), and high (Actinomyces dentalis, 72%) GC content, following Mann et al. (2018) (Fig. 1A). Using three different read length distributions (Fig. 1B), we generated short-read sequencing data from each reference genome using gargammel’s fragSim. To the resulting short-read sequences we added different amounts of aDNA damage using gargammel’s deamSim so that ten levels of damage ranging from 0% to 20% were observed, which were measured as the amount of observed C → T substitutions on the terminal base at the 5′ end of the DNA fragments (Fig. 1C). Finally, each of these 90 simulated datasets was subsampled to generate nine coverage bins ranging from 1-fold to 500-fold genome coverage by randomly drawing a coverage value from the uniform distribution defining each bin (Fig. 1D) and these were aligned to their respective reference genome using BWA aln (Li & Durbin, 2009) with the non-default parameters optimized for aDNA -n 0.01 -o 2 -l 16500 (Meyer et al., 2012).

Test contigs of different length were simulated by defining nine contig length bins ranging from 0.5 kb to 500 kb length (Fig. 1E) and randomly drawing 100 contig lengths from the respective uniform distribution defining each bin. Next, we chose the location of these test contigs by randomly selecting a contig from all contigs of sufficient length. We determined the exact location on the selected test contig from the reference genome by randomly drawing the start position from the uniform distribution defined by the length of the selected reference contig. This resulted in 900 test contigs per reference genome. Using these test contigs, we selected the aligned DNA fragments of the simulated sequencing data that overlapped the region defined by the contig and evaluated them using PyDamage. In total, we evaluated 702,900 test contigs (243,000 contigs for both M. smithii and T. forsythia, and 216,000 contigs for A. dentalis, for which no reference contig longer than 200 kb was available).

Archaeological sample

PyDamage implementation

PyDamage takes alignment files of reads (in SAM, BAM, or CRAM format) mapped against reference sequences (i.e., contigs, a MAG, a genome, or any other reference sequences of DNA). For each read mapping to each reference sequence j, using pysam (pysam developers, 2018), we count the number of apparent C → T transitions at each position which is i bases from the 5′ terminal end, i ∈ {0, 1, …, k}, denoted $N_i^j$ (by default, we set k = 35). Similarly we denote the number of observed conserved ‘C-to-C’ sites $M_i^j$, thus ${\displaystyle {\mathit{M}}^j=\left(M_0^j,\dots ,M_k^j\right)\text{and}{\mathit{N}}^j=\left(N_0^j,\dots ,N_k^j\right).}$ Finally, we calculate the proportion of C → T transitions occurring at each position, denoted $p_i^j$, in the following way: ${\displaystyle {{\hat{p}}_i}^j=\frac{N_i^j}{M_i^j+N_i^j}.}$

For D_i, the event that we observe a C → T transition i bases from the terminal end, we define two models: a null model ${\mathcal{M}}_0$ (Eq. (1)) which assumes that damage is independent of the position from the 5′ terminal end, and a damage model ${\mathcal{M}}_1$ (Eq. (2)) which assumes a decreasing probability of damage the further a the position from the 5′ terminal end. For the damage model, we re-scale the curve to the interval defined by parameters $\left[d_{pmin}^j,d_{pmax}^j\right]$. (1)${\displaystyle P_0\left(D_i\left|\right.p_0,j\right)=p_0{=}^{{\mathcal{M}}_0}{\pi }^j}$ (2)${\displaystyle P_1\left(D_i\left|\right.p_d^j,d_{pmin}^j,d_{pmax}^j,j\right)=\frac{\left(\right.\left[\right.{\left(1-p_d^j\right)}^i\times p_d^j\left]\right.-{\hat{p}}_{min}^j\left)\right.}{{\hat{p}}_{max}^j-{\hat{p}}_{min}^j}\times \left(d_{pmax}^j-d_{pmin}^j\right)+d_{pmin}^j{=}^{{\mathcal{M}}_1}{\pi }_i^j,}$ where ${\displaystyle {{\hat{p}}_{min}}^j\left(p_j^d\right)={\left(1-p_d^j\right)}^k\times p_d^j\text{and}{{\hat{p}}_{max}}^j\left(p_j^d\right)={\left(1-p_d^j\right)}^0\times p_d^j.}$

Using the curve fitting function of Scipy (Virtanen et al., 2020), with a trf (Branch, Coleman & Li, 1999) optimization and a Huber loss (Huber, 1992), we optimize the parameters of both models using $p_i^j$, by minimising the sum of squares, giving us the optimized set of parameters ${\displaystyle {\hat{\theta }}_0=\left\{{\hat{p}}_0\right\}\text{and}{\hat{\theta }}_1=\left\{{{\hat{p}}_d}^j,{{\hat{d}}^j}_{pmin},{{\hat{d}}^j}_{pmax}\right\}}$ for ${\mathcal{M}}_0$ and ${\mathcal{M}}_1$ respectively. Under ${\mathcal{M}}_0$ and ${\mathcal{M}}_1$ we have the following likelihood functions ${\displaystyle {\mathcal{L}}_0\left(\right.{\hat{\theta }}_0|{\mathit{M}}^j,{\mathit{N}}^j\left)\right.=\prod _{i=0}^k\left(\genfrac{}{}{0ex}{}{{M_i}^j+{N_i}^j}{{N_i}^j}\right){\left({}^{{\mathcal{M}}_0}{\hat{\pi }}^j\right)}^{{N_i}^j}{\left(1{-}^{{\mathcal{M}}_0}{\hat{\pi }}^j\right)}^{{M_i}^j},}$ ${\displaystyle {\mathcal{L}}_1\left(\right.{\hat{\theta }}_1\left|\right.{\mathit{M}}^j,{\mathit{N}}^j\left)\right.=\prod _{i=0}^k\left(\genfrac{}{}{0ex}{}{M_i^j+N_i^j}{N_i^j}\right){\left({}^{{\mathcal{M}}_1}{\hat{\pi }}_i^j\right)}^{N_i^j}{\left(1{-}^{{\mathcal{M}}_1}{{\hat{\pi }}_i}^{1,j}\right)}^{M_i^j},}$ where ${{}^{{\mathcal{M}}_0}\hat{\pi }}^j$ and ${{}^{{\mathcal{M}}_1}{\hat{\pi }}_i}^j$ are calculated using Eqs. (1) and (2). Note that if $d_{pmax}^j=d_{pmin}^j=p_0$, then ${}^{{\mathcal{M}}_0}{\pi }^j{=}^{{\mathcal{M}}_1}{\pi }_i^j$ for i = 0, …, k. Hence to compare the goodness-of-fit for models ${\mathcal{M}}_0$ and ${\mathcal{M}}_1$ for each reference, we calculate a likelihood-ratio test-statistic of the form ${\displaystyle {\lambda }_j=-2\text{ln}\left[\frac{{\mathcal{L}}_0\left(\right.{\hat{\theta }}_0\left|\right.{\mathit{M}}^j,{\mathit{N}}^j\left)\right.}{{\mathcal{L}}_1\left(\right.{\hat{\theta }}_1\left|\right.{\mathit{M}}^j,{\mathit{N}}^j\left)\right.}\right],}$ from which we compute a p-value using the fact that ${\lambda }_j\sim {\chi }_2^2$, asymptotically (Neyman & Pearson, 1933). Finally, we adjust the p-values for multiple testing of all references, using the StatsModels (Seabold & Perktold, 2010) implementation of the Benjamini–Hochberg procedure (Benjamini & Hochberg, 1995).

Statistical analysis and model selection

To test the performance of PyDamage in recognizing metagenome-assembled contigs with ancient DNA damage, we used the simulated short-read sequencing data aligned against simulated contigs of different lengths. Our method correctly identified contigs as not significantly damaged for simulations with no damage in 100% of cases. However, our model only correctly identified contigs as significantly damaged in 87.71% of cases where the contigs were simulated to have damage. To assess the performance of our method, and to determine the simulation parameters that most affected model accuracy, we analysed the simulated data using logistic regression via the glm function as implemented in the stats package using R (R Core Team, 2018). We included as potential explanatory variables the median read length, the simulated coverage, the simulated contig length, the simulated level of damage, and the GC content of each of the reference contigs, yielding 32 candidate logistic regression models.

We separated the data into two data sets: half of our data was used as ‘fit data’, data for performing model fit and parameter estimation, and the remaining half was reserved as ‘test data’, data that is used to assess model accuracy on data not used in fitting the model (n = 206, 831 in both cases). Unfortunately, with so many observations in our model, classical model selection methods such as AIC and ANOVA tend to overfit (Babyak, 2004). Similarly, we also performed ten-fold down-sampling of the data for each model such that we had equal numbers of damaged and undamaged simulations so as not to bias the predictive model. Hence, for each of the fitted 32 logistic regression models (with ϵ = 1 × 10⁻¹⁴ and maximum iterations 10³) we instead report the mean F₁ and Nagelkerke’s R² values for each candidate model.

Of the 32 candidate models, four models had both F₁ and R² values greater than 0.6 (see Table 1). Each these four models contained at least the following predictor variables: contig length, mean coverage, and the simulated level of damage. However, the full model with GC content and read length as additional predictor variables had similar F₁ and R² values, and so we consider all four models (see Fig. 2). Because it is possible that there is correlation between some of our predictor variables (i.e., increased levels of simulated damage could lead to a reduced median read length), we then performed a Relative Weights Analysis (RWA) to further estimate predictor variable importance in an uncorrelated setting (Chan, 2020). In essence, RWA calculates the proportion of the overall R² for the model that can be attributed to each variable. We performed RWA on both the full model and our best performing model. We found that the median read length and GC content accounted for only 0.31% and 2.75% of the R² value in the full model respectively. However, we found that contig length, mean coverage and the simulated level of damage all accounted for approximately one third of the R² value in our best performing model, indicating that these are the predictor variables of importance.

Our final logistic regression model identified mean coverage, the level of damage, and the contig length as significant predictor variables for model accuracy. Each of these variables had positive coefficients, meaning that an increase in damage, genome coverage, or contig length all lead to improved model accuracy. Each variable contributed about one third weight to the R² value in the model, indicating roughly equal importance in the accuracy of PyDamage. We integrated the best logistic regression model in PyDamage, with the StatsModels (Seabold & Perktold, 2010) implementation of GLM to provide an estimation of PyDamage ancient contig prediction accuracy given the amount of damage, coverage, and length for each reference (Fig. 3), and found these predictions to adequately match the observed model accuracy for our simulated data set (Fig. 4).

Application of PyDamage to archeological samples

To test PyDamage on empirical data, we assembled metagenomic data from the paleofeces sample ZSM028 with the metaSPAdes de novo assembler. We obtained a total of 359,807 contigs, with an N50 of 429 bp. Such assemblies, consisting of a large number of relatively short contigs, are typical for de novo assembled aDNA datasets (Wibowo et al., 2021). After removing sequences shorter than 1,000 bp, 17,103 contigs were left. PyDamage (revision 099fd34) was able to perform a damage estimation for 99.75% of these contigs(17,061 contigs). Because the ZSM028 sequencing library was not treated with uracil-DNA-glycosylase (Rohland et al., 2015), nor amplified with a damage suppressing DNA polymerase, we expect a relatively shallow DNA damage decay curve, and thus filtered for this using the $p_d^j$ parameter. We chose a prediction accuracy threshold of 0.67 after locating the knee point on Fig. 5 with the kneedle method (Satopaa et al., 2011). After filtering PyDamage results with a q-value ≤0.05, $p_d^j\le 0.6$, and prediction accuracy ≥0.67, 1,944 contigs remain. The 5′ damage for these contigs ranges from 4.0% to 45.1% with a mean of 14.3% (Fig. 6). Their coverage spans 6.1X to 1, 579.8X with a mean of 65.6X, while their length ranges from 1,002 bp to 90,306 bp with a mean of 5,212 bp and an N50 of 10,805 bp.

The Kraken2 taxonomic profile of the microbial contigs identified by PyDamage identified as ancient (Fig. 7) is consistent with bacteria known to be members of the human gut microbiome, including Prevotella (239 contigs), Treponema (166 contigs), Bacteroides (103 contigs), Lachnospiraceae (119 contigs) Blautia (36 contigs), Ruminococcus (25 contigs), Phocaeicola (18 contigs) and Romboutsia (16 contigs) (Schnorr et al., 2016; Pasolli et al., 2019; Singh et al., 2017), as well as taxonomic groups known to be involved in initial decomposition, such as Clostridium (145 contigs) (Hyde et al., 2017; Harrison et al., 2020; Dash & Das, 2020). In addition, eukaryotic contigs were assigned to humans (18 contigs), and to the plant families Fabaceae (18 contigs) and Solanaceae (18 contigs), two families of economically important crops in the Americas that include beans, tomatoes, chile peppers, and tobacco. The remaining contigs were almost entirely assigned to higher taxonomic levels within the important gut microbiome phyla Bacteriodetes, Firmicutes, Proteobacteria, and Spirochaetes, as well as to the Streptophyta phylum of vascular plants. Collectively, these five phyla accounted for 1,283 of to 1,494 contigs that could be taxonomically assigned.

Functional annotation of the authenticated ancient contigs using Prokka was successful for 1,901 of 1,944 contigs. Among these, multiple genes of functional interest were identified, including contigs annotated as encoding the multidrug resistance proteins MdtA, MdtB, and MdtC, which convey, among other functions, bile salt resistance (Nagakubo et al., 2002) (Table 2). Kraken2 taxonomic profiling of these three contigs yields a taxonomic assignation to the gut spirochaete Treponema succinifaciens, a species absent in the gut microbiome of industrialized populations, but which is found globally in societies practicing traditional forms of subsistence (Obregon-Tito et al., 2015; Schnorr et al., 2014). Other authenticated contigs contained genes associated with resistance to the natural antimicrobial compounds fosmidomycin, colistin, daunorubicin/doxorubicin, tetracycline, polymyxin, and linearmycin. A growing body of evidence supports an ancient origin for resistance to most classes of natural antibiotics (D’Costa et al., 2011; Warinner et al., 2014; Christaki, Marcou & Tofarides, 2020; Wibowo et al., 2021).