{'fraction_non_alphanumeric': 0.017719491440584642, 'fraction_numerical': 0.027680448493342676, 'mean_word_length': 4.759296627270107, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Improved assemblies and comparison of two ancient Yersinia pestis genomes

Yersinia pestis is the causative agent of the bubonic plague, a disease responsible for several dramatic historical pandemics. Progress in ancient DNA (aDNA) sequencing rendered possible the sequencing of whole genomes of important human pathogens, including the ancient Yersinia pestis strains responsible for outbreaks of the bubonic plague in London in the 14th century and in Marseille in the 18th century among others. However, aDNA sequencing data are still characterized by short reads and non-uniform coverage, so assembling ancient pathogen genomes remains challenging and prevents in many cases a detailed study of genome rearrangements. It has recently been shown that comparative scaffolding approaches can improve the assembly of ancient Yersinia pestis genomes at a chromosome level. In the present work, we address the last step of genome assembly, the gap-filling stage. We describe an optimization-based method AGapEs (Ancestral Gap Estimation) to fill in inter-contig gaps using a combination of a template obtained from related extant genomes and aDNA reads. We show how this approach can be used to refine comparative scaffolding by selecting contig adjacencies supported by a mix of unassembled aDNA reads and comparative signal. We apply our method to two data sets from the London and Marseilles outbreaks of the bubonic plague. We obtain highly improved genome assemblies for both the London strain and Marseille strain genomes, comprised of respectively five and six scaffolds, with 95% of the assemblies supported by ancient reads. We analyze the genome evolution between both ancient genomes in terms of genome rearrangements, and observe a high level of synteny conservation between these two strains.

fragmented assemblies of extant strains to analyze the evolution of pathogen genome organization.

48
Without long-read sequencing data, comparative scaffolding based on the comparison of the contigs of a genome of 49 interest with related assembled genomes has proven to be a useful approach to improve the assembly of fragmented 50 genomes, especially bacterial genomes 8, 25, 39, 41 . Among such methods, FPSAC 39 was introduced to improve ancient 51 genome assemblies within a phylogenetic context. It was applied to aDNA contigs from the Yersinia pestis strain 52 responsible for the medieval London bubonic plague outbreak -that was shown to be ancestral to several extant 53 Yersinia pestis strains 7 -and resulted in an improvement of the initial contig assembly from thousands of contigs to 54 a chromosome-scale scaffolding. Moreover, taking advantage of the high sequence conservation in Yersinia pestis 55 genomes, the inter-contigs gaps of the ancient Yersinia pestis strain were filled with putative sequences reconstructed 56 from multiple sequence alignments of conserved extant gaps. This gap-filling step shed an interesting light on genomic 57 features hidden within the assembly gaps, in particular IS and their correlation with rearrangement breakpoints reuse, 58 but also allowed the potential reconstruction of regions that were not recovered or were absent from the aDNA material. 59 However, the scaffolding of adjacencies and gap sequences obtained in 39 , that accounted for roughly 20% of the genome 60 size, were inferred through computational methods within a parsimony framework, that can be sensitive to convergent 61 evolution that cannot be ruled out for genomes with a high rate of genome rearrangements such as Yersinia pestis 11 . 62 In the present work, we address this issue by using the large set of aDNA reads that are unassembled after the contig 63 assembly stage, to confirm the scaffolding of contigs as well as sequences for inter-contigs gaps. We introduce the 64 method AGapEs (Ancestral Gap Estimation) which attempts to fill the inter-contig gap between two adjacent ancient 2 Materials and Methods 78 We first describe the input to our analysis, namely ancient sequencing data, ancient and extant assemblies and annotations 79 of IS, before outlining the general pipeline we used to improve the assembly of the ancient genomes.  Figure S3). We rely on seven extant Yersinia pestis and four Yersinia pseudotuberculosis as 84 reference and outgroup genomes (see Table S1). The phylogeny of the considered strains is depicted in Figure S1 and is 85 taken from 6, 7 .

86
Contig assembly and preprocessing. We de novo assembled aDNA reads into contigs using Minia 10 for both 87 aDNA data sets (London outbreak and Marseille outbreak). Minia is a conservative assembler based on an efficient 88 implementation of the de Bruijn graph methodology. In general, Minia produces shorter contigs than competing 89 assemblers, as it avoids assembly decisions in case of ambiguity in the sequence data. We will refer to the Minia 90 assemblies as de novo assemblies in the following. To allow the comparison with extant genomes, contigs above a 91 minimum length threshold were aligned with the extant genomes to define families of homologous synteny blocks 92 (called markers from now) as described in 39 . Marker families were then filtered to retain only one-to-one orthologous 93 families, i.e. families that contain one and exactly one marker in each considered extant and ancient genome.

94
Insertion sequence annotation. Insertion Sequences (IS) are strongly related to rearrangements in Yersinia pestis 95 evolution, and their annotation in the considered extant genomes is crucial. In order to annotate IS, we designed our own 96 annotation pipeline. Because IS elements in the original Genbank files were rather disparately annotated, we relied on 97 automated annotations from the Basys annotation server 49 . Basys identified 11 families of IS transposase proteins (see 98   Table S2). For each of these families, we produced a multiple alignment of their annotated sequences using muscle 15 99 which was subsequently used to train Hidden Markov Model (HMM) profiles. Using hmmer 14 , we then annotated 100 those regions as associated to IS elements that showed significant correlation to any of the HMM profiles. We eventually 101 combined the Genbank annotations with these derived annotations. The number of these IS annotations per reference 102 genome ranges from 151 in Yersinia pestis KIM10+ to 293 in Yersinia pestis Antiqua (see Table S1). The length of the 103 annotations ranges from 60bp to 2,417bp; some short annotations deviate from the expected length for IS, however, in 104 order to avoid filtering any true annotations, we include them all as potential IS coordinates in the downstream analysis.

105
Ancestral marker adjacencies. Each marker can be defined by a pair of marker extremities. An adjacency consists 106 of two markers extremities that are contiguous along a genome, i.e. are not separated by a sequence containing another 107 marker. For extant genomes, extant adjacencies can be observed directly, while for an ancestral genome of interest, 108 we infer potential ancestral adjacencies using the Dollo parsimony principle 39 : two ancient marker extremities are 109 potentially adjacent if there exist two extant genomes whose evolutionary path contains the most common recent 110 ancestor of the London and Marseille strains and where the two corresponding extant marker extremities are contiguous 111 (see Figure S5 for an example). Hence every potential ancestral adjacency is supported by a set of extant adjacencies. A 112 gap is the sequence between the two marker extremities defining an adjacency. Therefore each putative ancestral gap is 113 likewise supported by a set of extant gap sequences. 114 We say that two potential ancestral adjacencies are conflicting if they share a common marker extremity. An 115 IS-annotated adjacency is supported by at least one extant adjacency whose gap contains an IS annotation. An adjacency 116 that is neither conflicting nor IS-annotated is said to be simple. introduce is a template-based method to assess the validity of a potential ancestral adjacency. The general principle is to 119 associate to every ancestral gap a template sequence obtained from the supporting extant gaps sequences. We can then 120 map aDNA reads onto this template and assemble the mapped aDNA reads into a sequence that minimizes the edit 121 distance to the template sequence. The rationale for this template-based approach is that, due to the low coverage of 122 the aDNA reads and their short length, existing gap-filling methods fail to fill a large number of ancestral gaps. For example, the method gap2Seq 42 , a recent efficient gap-closing algorithm based on finding a path of given length in a de We describe now the AGapEs algorithm. Assume we are given a template sequence t for a gap in an adjacency a = {m 1 , m 2 } between two marker extremities. We define R = m 1 + t + m 2 as the concatenated nucleotide sequence of 127 the oriented markers and the respective template. We first align the aDNA reads onto R, using BWA 29 , where we only 128 consider mappings whose start and/or end position is in t (i. e. either fully included in t or overlapping the junction 129 between a marker and the gap template). Next, we construct a graph G where vertices are mappings m 2 M t and there 130 is an edge between two vertices (i. e. mappings) if the two mapping coordinates (segments of R) overlap. For each 131 such edge/overlap, we define s as the non-overlapping suffix of the mappings with the highest end coordinate. We 132 can then associate a weight to each edge given by the edit distance between s and the subsequence R s of R it aligns to.

133
A sequence of overlapping reads that minimizes the distance to t can then be found by searching for a shortest path 134 between the vertex labeled with the smallest start position (i. e. the first mapping covering the junction between m 1 and 135 t) and the vertex labeled with the largest start position (i. e. the last mapping covering the junction between t and m 2 ).

136
See Figure S8 for an illustration. bases in the sequence are covered, but not both by the same read. In these cases, uncovered regions and breakpoints 141 need to be identified in the mapping beforehand to identify start and end vertex of the shortest path. We can then obtain 142 a partial gap filling, precisely for the regions covered by mapped reads.

143
The key element of the approach described above lies in defining the template sequence or set of alternative template 144 sequences associated to each ancestral gap. We follow the general approach described in 39 , that computes a multiple 145 sequence alignment of the supporting extant sequence gaps and applies the Fitch-Hartigan parsimony algorithm 16 to 146 each alignment column to reconstruct a most parsimonious ancestral sequence. If the multiple sequence alignment of 147 extant gaps shows little variation, as is the case for most gaps in our data sets, then a single template sequence can be 148 considered, as we expect that minor variations compared to the true ancestral sequence (substitutions, small indels) will 149 be corrected during the local assembly process outlined above. Alternatively, if larger variations are observed, such 150 as larger indels or a contradicting pattern of presence/absence of an IS in the supporting extant gaps, then alternative 151 templates can be considered, under the hypothesis that the true variant can be recovered from the mapped aDNA reads.

152
Hence in the following analysis, we separate all potential ancestral gaps into groups of simple, conflicting and 153 IS-annotated gaps. For simple and conflicting gaps without IS annotation, we can follow the process described above 154 directly. For IS-annotated gaps, we reduce the described large variations in the multiple alignment by further dividing 155 its supporting extant gaps into sets of IS-annotated and non-IS-annotated sequences respectively. Building the multiple 156 alignment on each of these sets separately allows us to define two alternative templates that can be used as a basis to fill 157 the gap. Ideally, differences in read coverage or breakpoints naturally identified by AGapEs then point to one of the 158 alternative templates for each IS-annotated gap. Further, for each template that is only partly covered by mapped reads, 159 we will correct the covered parts according to the reads using AGapEs and use the template sequence otherwise.

160
The implementation of AGapEs is available at http://github.com/nluhmann/AGapEs, the data underly- are 21 gaps whose lengths difference falls into the length range of potential annotated IS elements, thus raising the 176 question of the presence of an IS within these adjacencies in the ancestral genome. We note a small number of 5 177 potential ancestral adjacencies with strikingly large extant gap length differences. All of these gaps accumulate more 178 than one IS annotation in some extant genomes. Most problematic is a gap with length difference of more than 100.000 179 bp. As this gap is not well conserved in general (apart from the inserted sequences), it is difficult to obtain a good 180 template sequence based on a very fragmented multiple alignment. We will get back later to this special gap.

181
Ancestral gaps filling. We applied AGapEs to all potential ancestral gaps. We assume a gap to be filled, if we find a 182 sequence of reads that covers the whole ancestral gap. As we test two alternative templates for an IS-annotated gap, we mapped aDNA reads, we correct the covered regions as described above and use the template sequence of the uncovered 188 regions to complete filling the gap. Figure 1 summarizes the gap-filling results (see also Table S5).

189
A high number of gaps is supported by sufficient read coverage that enables us to fill the gap with a sequence of 190 overlapping aDNA reads. Especially considering partially covered gaps improves the length of the genome that is 191 supported by reads. Note that we also find covering reads for all gaps of length 0, spanning the breakpoint between 192 directly adjacent markers. 193 We further computed the edit distance between each reconstructed gap sequence and its previous gap template. For to remove all unsupported conflicting adjacencies. See Figure S12 for the read coverage of discarded adjacencies. Note 209 that filling these gaps only partially does not provide much information, as uncovered regions can be either breakpoints 210 or correspond to regions of the ancestral genome that were not sequenced.

211
The set of ancestral adjacencies can then be ordered into five Contiguous Ancestral Regions (CARs). We converted 212 the reconstructed sequences of markers back to genome sequences by filling the gaps with the read sequences if possible 213 and resorting to the template sequence otherwise.

214
As mentioned earlier, we observe one gap with highly differing extant gap lengths and very little conservation in the 215 reconstruction. The multiple alignment based on extant gap sequence is very fragmented and the mapping of reads onto 216 this template is poor: the gap contains 211 uncovered regions of 9, 319 bp in total. See Figure S13 for an overview over 217 the read coverage for this gap in the de novo assembly. As the reconstructed sequence has a high edit distance after 218 partial gap filling, we remove this gap sequence completely at this point to avoid dubious and non-robust reconstructed 219 ancestral sequences.

228
In the improved assembly, 49.88% of the sequence is based on markers and hence directly adopted from the 229 initial assembly. Together with the gaps that have been filled by read sequences, we can say that in total 95.25% are 230 reconstructed using only the available aDNA reads.

243
We can see in Figure 4 that with the combined set of reads, we can fill nearly all simple gaps by read sequences. In identified the same number of potential ancestral IS as for the London strain. 247 We identified two conflicting components in this set of potential adjacencies (see Figure S14). Both of them align in  Figure S15.

255
This results into 6 CARs for the ancestral genome. Again, we used BWA 29 to align reads from all five samples again 256 to the assembly to assess the amount of uncovered regions in the reconstructed sequences. In total, only 54,672bp in 257 this mapping are not covered by any read and the length of the uncovered regions is rather short (see Figure S17).   CO92 as a reference. In order to assess the influence of the reference sequence in the assembly of the ancient genome, 279 we compare our pipeline using this initial assembly to our results based on the de novo assembly. 280 We compared the two sets of CARs obtained from both initial assemblies by aligning the resulting genome sequences 281 using MUMmer 28 . We observe no rearrangements between both resulting sets of CARs (see also Figure S16), showing 282 that, in terms of large-scale genome organization, the final result does not depend on the initial contig assembly.

283
Assembly validation We compare our results to assemblies obtained with several other assembly pipelines. We 284 used the iMetAMOS pipeline 26 to determine the best de novo assembly for both data sets testing different assemblers.

285
The winning assembly computed by SPAdes 3 for both data sets as well as the minia assemblies on both datasets were 286 subsequently used as input for two comparative scaffolding programs, Ragout 25 and MeDuSa 8 , to obtain a scaffolding 287 We will distinguish the results according to the scaffolding tool used in the following.

289
As shown in Table 1, for both datasets, Ragout is reconstructing the smallest number of contigs, however the 290 scaffolds still contain a high number of unfilled gaps that cannot be closed by gap2Seq. See Table S10  SPAdes and minia assemblies for the presence of 10 Yersinia pestis species-specific IS elements found in the ISFinder 300 database 45 and using all potential IS gaps as references.

301
While ISseeker is not able to annotate IS elements in the Minia assembly, we see in Table 2 30 annotations that are 302 found in the SPAdes assembly. Seven of these are not annotated in the AGapEs reconstruction, and they all concern 303 gaps that are only partially covered by reads. However a manual check of these gaps determined the presence of the 304 respective IS element in five gaps, indicating that ISseeker was not able to correctly annotate these elements in these 305 cases.

307
In this paper, we present a method to fill the gaps between contigs assembled from aDNA reads that combines 308 comparative scaffolding using related extant genomes and direct aDNA sequencing data, and we apply it to two ancient 309 Yersinia pestis strains isolated from the remains of victims of the second plague pandemic.

310
The comparison of the two assemblies for the London strain illustrates that relying on a shorter initial de novo contig 311 assembly does not impact significantly the final result. The results we obtain with the Marseille data set illustrates that 312 Applied to the same data set for the London strain, the method FPSAC 39 was able to obtain a single scaffold 318 based on parsimonious optimization. Comparing our resulting assembly to this single scaffold, we can identify two 319 breakpoints between both assemblies, hence both methods do not entirely support the same scaffold structure for the 320 London strain. These disagreements should be seen as weak points in both assemblies, as they are not reconstructed by 321 different scaffolding objectives and would need to be confirmed more confidently by additional sequencing data. 322 We see a clear connection between conflicts in the set of potential adjacencies and the presence of IS elements in 323 the corresponding gaps. Solving these conflicts based on aDNA read data provides a useful way to identify ancestral