{'fraction_non_alphanumeric': 0.02338337182448037, 'fraction_numerical': 0.02594028373474101, 'mean_word_length': 5.31155648099948, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Accuracy through Subsampling of Protein EvolutioN : Analyzing and reconstructing protein divergence using an ensemble approach

Mapping the history of gene duplications which gave rise to a protein family encoded in a genome (a set of paralogs) can be critical to understanding how those proteins function in their host cells today. However, since each member of a family is recapitulated in the genomes of related species (a set of orthologs), selection of sequences to be included in the history reconstruction is non-trivial. Reconstruction is extremely sensitive to the choice of sequences, which is deeply problematic given no mechanism exists for assessing the accuracy of individual reconstructions. Here, we capitalize on the variability of phylogenetic tree reconstruction to selected input sequences, by subsampling from the available ortholog sequences of a protein family to create an ensemble of trees, which explores the space of plausible tree topologies. We hypothesize that the most consistent topological features across an ensemble are more likely to be true and propose a tree reconstruction algorithm (ASPEN) based on this hypothesis. We simulate 600 protein families over known phylogenies, with varying branch lengths, and compare the accuracy of ASPEN reconstructions to those of traditional phylogeny inference methods. We find that ASPEN trees are more accurate than trees reconstructed traditionally. Additionally, we develop an observable metric calculated form subsampling, reconstruction Precision, for assessing the likely accuracy of a traditional, single-alignment all-sequence reconstruction of the divergence history for a set of paralogs. Together these findings suggest that an ensemble of imperfect reconstructions can provide more accurate insight than any individual reconstruction.

Protein families grow in size and diversity through duplication of genes encoding existing family members followed by 2 functional divergence of the duplicates [1,2]. Immediately following a gene duplication event the affected genome 3 contains two identical copies of the duplicated gene. Because the genes are redundant, relaxed purifying selection 4 allows mutations to accumulate rapidly. Since the added energy cost of expressing identical products from redundant 5 loci confers a selective disadvantage, mutations resulting in loss of functionality by one of the copies are typically 6 favored by selection. However, the rapid accumulation of mutations can also result in partial or complete functional divergence reconstructions for well-studied protein domain families [20][21][22] have been relied upon extensively by the 19 scientific community. Because such reconstructions are created from single sequence alignments, they ignore the great 20 deal of uncertainty in topology reconstruction under equally valid alignment representations of input sequence data. 21 Divergence topology reconstruction is extremely sensitive to the input alignment. For example, the same 22 sequences aligned by different algorithms [23][24][25][26][27] or using different guide trees [28] yield different topology 23 reconstructions. So does reversing input sequences prior to alignment [29,30], or removing less than 0.1% of columns 24 from an alignment containing over 600,000 columns [31]. For paralog divergence topologies, another source of 25 uncertainty likely to influence reconstruction is the set of orthologs selected to represent each paralog. Because 26 duplications usually predate numerous speciation events, they tend to correspond to deep internal nodes -nodes with 27 many descendant leaves -in full phylogenies of protein families. MRCAs of ortholgs descend from duplications 28 ( Figure 1), meaning every ortholog descended from each MRCA is also descended from the duplication. Deep internal 29 nodes tend to be most sensitive to perturbations of the input alignment [32]. Unfortunately, since the true history of 30 protein divergence is hidden from us in time, we have no way of knowing which divergence topologies are more 31 accurate, given the equal validity of input alignments. 32 Although traditional tree reconstruction produces phylogenies -topologies parametrized with branch lengths 33 reflecting extent of divergence -we disregard the branch lengths here to focus on the topologies alone. In traditional 34 inference topologies and branch lengths are inferred jointly, alternating between topology modifications and branch 35 length optimization in the case of statistical (Maximum Likelihood and Bayesian) methods. Because the likelihood 36 function is evaluated many times for each proposed topology, and topology space is almost unfathomably large, 37 statistical methods can suffer extremely long run times on large sequence collections. However, if accurate candidate 38 topologies can be identified by other means, the computational cost of optimizing branch lengths for individual 39 topologies is nearly trivial, while optimization for multiple topologies is embarrassingly parallel. Our approach 40 permits separating topology reconstruction from branch length optimization. 41 Furthermore, we focus on reconstructing only the hardest topology nodes -the deep internal nodes corresponding 42 to protein or domain paralog divergence. We treat MRCAs of ortholog sets as leaves in our reconstructions and 43 disregard ortholog divergence, which overwhelmingly recapitulates the species tree. Species divergence is 44 reconstructed more accurately by other approaches [31,32]. Instead, we capitalize on the variance in reconstructed 45 topologies under changes in ortholog representation of paralogs to separate topological features we believe to be 46 supported by phylogenetic signal from spurious ones we believe to result from noise. We hypothesize that features 47 observed more frequently under ortholog resampling are more likely to reflect signal and, therefore, be more accurate, 48 than less frequently observed ones. We explore the relationship between accuracy and variability in reconstructing 49 paralog divergence topologies and propose a metric for assessing the likely accuracy of a single-alignment 50 reconstruction for a given protein family. We then present ASPEN, a topology reconstruction algorithm that 51 integrates over the uncertainty of single alignment reconstructions to build and rank trees according to observations 52 across reconstructions from many equally valid alignments. ASPEN produces more accurate topologies than 53 Figure 2. Analysis framework for comparing reconstruction Accuracy and Precision. (A) Sequence evolution was simulated over synthetic phylogenies. Synthetic phylogenies were pruned to MRCAs of ortholog sets and branch lengths were discarded to obtain true paralog divergence topologies. Simulated sequences were aligned, phylogenies were inferred from those alignments, and "all-sequence" reconstructions of paralog divergence topologies were extracted. (B) Sequences were repeatedly sampled from each ortholog set in a family and phylogeny inference and topology extraction were done to produce a "subsampled topology". Repeating this N times yields an ensemble of topologies. (textbfC) We define Accuracy as the similarity between the all-sequence reconstruction and Precision as the comparison between subsampled topologies and the all-sequence topology.
traditional reconstructions from single, all-sequence alignments.

54
Experimental framework for reconstruction analysis 55 We generated test sequence data by simulating evolution of protein families instead of using natural protein sequences 56 for two previously noted reasons [25]. First, simulating evolution over known phylogenies allowed us make a 57 quantitative assessment of reconstruction accuracy compared to the "true" divergence topology. Second, it allowed us 58 to explore a range of divergence conditions by systematically varying branch lengths of input phylogenies, while 59 controlling for other factors such as overall sequence length and the distribution of secondary structure elements and 60 disordered loops. Assembling a comparable biological data set would have been impossible. 61 We simulated families containing 15 paralogs, each represented by 66 orthologs. In order to make the assessment 62 statistically robust, we generated 600 families across a range of post-duplication branch lengths. An alignment of 63 human tyrosine kinase domains (median length 269 a.a.) was used as template for all simulations (see Methods for 64 simulation details). We then used all combinations of three multiple sequence alignment algorithms (MAFFT's 65 L-INS-i protocol [33], ClustalOmega [34], and Muscle [35]) and two phylogeny inference algorithms (FastTree2 [36] 66 and RAxML [37]) to reconstruct phylogenies for the 600 simulated families. We compared the reconstructed paralog 67 divergence topologies, excluding speciation nodes by pruning orthologs' MRCAs to leaves, to the true divergence 68 topology over which evolution was simulated ( Figure 2A). We quantified topology differences with the 69 Robinson-Foulds symmetric distance metric [38], modified to account for the occasionally non-monophyletic  Figure S1 ). Based on these results, we selected the combination of L-INS-i and FastTree2 for all remaining analysis. 75 Subsampling reveals an observable measure of accuracy 76 Given the known sensitivity of reconstruction to input alignment, we explored reconstruction variance resulting from 77 differences in ortholog representation of paralogs using the framework outlined in Figure 2. We gathered the sets of 78 ortholog sequences representing each paralog in a simulated family ( Fig. 2A) and performed a resampling experiment 79 (Fig. 2B): 50 times we randomly sampled 60 of 66 sequences (91%) from each ortholog set and performed traditional 80 reconstruction, using L-INS-i and FastTree2 with each collection of subsampled input sequences. We retained a large 81 fraction of sequences to minimize both the input variation and the loss of phylogenetic signal. To quantify 82 reconstruction uncertainty, we measured the similarity (1 the average of RF ⇤ ) between topologies reconstructed 83 from most of the sequences to the "all-sequence" topology ( Fig. 2C). Since this quantity is a measure of how close 84 the estimates are to each other, we refer to it as Precision.  reconstruction Accuracy (unknowable for natural proteins) and, by proxy, the overall "complexity" of reconstruction 89 for that family. Importantly, this also suggests that our 600 synthetic protein families span a range of complexities, 90 allowing us to observe the performance of reconstruction as a function of complexity, via its proxy -Precision.

91
Using variability to distinguish phylogenetic signal from noise 92 Although we observed high reconstruction Precision for many families, only four of 600 families had identical paralog 93 divergence topologies reconstructed from every subsampled alignment (Precision=1). Even among families with the 94 highest Precision, and under dense subsampling, reconstruction variability was pervasive. On the other hand, 95 Salichos and Rokas [32] argued that pairwise RF distances smaller than 1 (the average RF distance among randomly 96 generated topologies) indicates consistent phylogenetic signal among the topologies being compared. Most of our 600 97 families had Precision (1 hRF ⇤ i) significantly greater than 0, but less than 1. Thus we sought to go a step further 98 and test our central hypothesis: not only does intermediate Precision indicate consistent signal, but more frequently 99 recapitulated features are more likely to be accurate, and this fact can be used to reconstruct more accurate 100 topologies. In order to test this we first needed a way to extract frequently recapitulated features, and then a way to 101 identify topologies most consistent with those features. Next we describe our method, ASPEN, which accomplishes 102 both tasks. 103

4/15
Reconstructing topologies from ensemble sampling 104 We created a method we call ASPEN, for Accuracy through Subsampling of Protein EvolutioN, to construct and 105 score topologies according to their consistency with topological features frequently represented in an ensemble of 106 subsampled reconstructions (Fig. 2B). It relies on two key innovations: 1) extraction of topological features from an 107 ensemble as frequencies of path lengths between leaves, and 2) an algorithm to construct and score topologies 108 according to their consistency with observed path length frequencies.

109
Transforming topology sets into path length distributions 110 ASPEN's foundation is the equivalent representation of a topology (an acyclic, bifurcating graph) as a matrix of path 111 lengths between leaves in terms of the number of internal nodes encountered along a path. First we demonstrate 112 equivalence of graph and matrix representations by presenting a simple algorithm for interconverting between the two 113 ( Figure 4). Then we discuss how ensembles of topologies are transformed into path length frequency distributions.  Traversal is repeated as necessary until construction is completed. Because all path lengths are derived from a single 129 topology, they are guaranteed to be consistent, making the construction unambiguous.

130
Generating path length frequency distributions 131 We take advantage of the alternate matrix representation to capture the individual variation of each leaf-to-leaf path 132 length across an ensemble of topologies. Each topology is transformed into a path lengths matrix. Then path lengths 133 for each pair of leaves are aggregated into a path length distribution for that pair ( Figure 4C). Although ortholog sets 134 overwhelmingly group into monophyletic subtrees across ensemble topologies (their MRCAs have no descendant leaves 135 besides themselves), occasionally reconstructions do yield non-monophyletic ortholog sets. Because this violates an 136 underlying assumption of the reconstruction, as well as the true topology of each synthetic protein family, we preclude 137 paths compromised by this incorrect reconstruction from contributing to path length distributions: the length of any 138 leaf-to-leaf path that contains a compromised internal node is not included in the distribution for that leaf pair.

139
Path length frequencies guide topology reconstruction 140 A score reflecting consistency with extracted features 141 ASPEN uses a quantitative metric for measuring the consistency of a proposed topology with observations from an ensemble of topologies. The score assigned to a topology is expressed in terms of log frequencies of leaf-to-leaf path lengths, log(f L pair ) where L is the length of path between leaves in pair, incorporated into the topology: This scoring function rewards incorporation of frequently observed path lengths and penalizes rarely observed path 142 lengths.

143
A branch-and-bound topology construction algorithm 144 Using the bottom-up procedure for constructing a topology graph from its path lengths matrix representation 145 ( Figure 4B), we developed an algorithm that uses a branch-and-bound strategy to construct the requested number of 146 highest-scoring topologies according to the scoring function above. We describe the branching and bounding 147 procedures in the next two sections.

148
Branching By analogy with the single-topology procedure in Figure 4B, construction of internal nodes is triggered 149 by path length entries encountered during list traversal. However, this list contains every observed path length for 150 every leaf pair, together with its frequency ( Figure 4D). Unlike the single-topology case, list entries cannot be 151 assumed to be consistent with each other. In fact, many combinations of path lengths on the list cannot be 152 incorporated into one topology. For example, for hypothetical leaves A, B, and C, path lengths (A $ B, 1) and

153
(B $ C, 1) are mutually exclusive because in a bifurcating topology B can be one internal node removed from either 154

7/15
A or C, but not both. In single topology reconstruction, if a path length could be completed by the introduction of 155 an internal node, that node could be safely constructed because it was guaranteed to satisfy every other list entry.

156
Since that guarantee no longer holds, multiple topologies are constructed simultaneously by allowing the construction 157 path to branch ( Figure 5).

158
"Assemblies" are used to track simultaneous reconstruction of multiple topologies. Each assembly holds a copy of 159 the path length frequencies list, a partially constructed topology, and the current topology score according to the 160 scoring function (discussed below in the section on bounding). Reconstruction proceeds in iterations, starting with a 161 single empty assembly ( Figure 5, left). On the first iteration, the entire list is traversed and every possible extension 162 by introduction of a new node is created simultaneously in a copy of the original assembly ( Figure 5, middle). In each 163 new assembly, all path lengths completed by the new node and all path lengths incompatible with it are marked and 164 not re-examined on subsequent iterations. Remaining path lengths are not completed by the new node, but remain 165 compatible with it. On subsequent iterations the same procedure is repeated for all tracked assemblies. 166 In principle, branching and iteration alone yield every topology consistent with path lengths observed in the 167 ensemble. In practice, this results in a combinatorial explosion which must be carefully managed to allow 168 construction to proceed to completion. First, Figure 5 (right) demonstrates how branching to satisfy non-conflicting 169 path lengths can lead to collisions between diverged construction paths on later iterations. This occurs because many 170 topologies can be constructed by introducing internal nodes in multiple orders. Each branched path represents a 171 particular order of internal node introduction. In a practical implementation collisions must be managed in order to 172 prevent multiple reconstructions of the same topology by multiple paths -an enormous replication of effort.

173
Second, even if each distinct topology is constructed once, in most cases reconstructing every topology consistent 174 with observations from the ensemble, no matter how infrequent, is neither practical nor useful. Bounding, described 175 in the next section, guarantees reconstruction of only the requested number of top scoring topologies.

176
Bounding The score is used to rank completed topologies, where ranking is updated every time a new topology is finished. The number of top scoring topologies to reconstruct, X, is requested at the beginning of a reconstruction run (10,000 was used in ASPEN evaluation). Once the initial X topologies are constructed, the Xth topology score constitutes the bound. Partially constructed topologies are abandoned if no complete topology that can be derived from that construction state will score above the bound. We determine this by calculating the score for already-incorporated path lengths and projecting the best possible score for a complete topology by assuming the most frequent remaining path length will be incorporated for every unconnected leaf pair:  increases (lower Precision bins, Figure 6E-G), indicating less discriminatory power with respect to accuracy.

195
Nevertheless, ASPEN's top-ranked topology is, on average, also its most accurate across all Precision bins.

196
Top ASPEN topology beats all-sequence reconstructions 197 Next, we compared ASPEN's best topology to all-sequence single-alignment reconstructions ( Figure 6A). Like all 198 other methods, ASPEN's accuracy is a function of Precision, or difficulty of the reconstruction task. As discussed For families with high Precision, where one may reasonably expect to reconstruct an accurate topology, ASPEN's top, 219 or top few topologies are likely more accurate than any single-alignment reconstruction. One or a few of these 220 topologies can be confidently used for downstream applications. This result is far from trivial, given that ASPEN's 221 subsampling approach scales far better with the overall number of input sequences than traditional statistical 222 reconstruction methods. With the advent of affordable genome sequencing and the resulting explosion in the number 223 of sequenced and annotated species' genomes [42,43], all-sequence reconstruction of paralog divergence by statistical 224 methods has become infeasible for many families with large numbers of orthologs. Therefore, subsampling large 225 samples of orthologs to yield a Precision score can now be used to identify how likely the full sequence topology is to 226 be accurate, determining if one is working in a high or low Precision/Accuracy regime. Subsampling in the process of reconstruction proved to be extremely powerful -it identified two measures (Precision 238 and Score based on observed frequencies) of something unknowable (Accuracy) and guided a reconstruction method 239 that identifies much more accurate topologies than traditional approaches. That ASPEN reconstructions were more 240 accurate than single-alignment reconstructions, is evidence that the central hypothesis of this work is supported -241 relationships found consistently amongst the variance produced by subsampling are more likely to be reflective of true 242 protein divergence histories. We anticipate that, as a meta analysis approach to tree evaluation and reconstruction, 243 ASPEN is likely to continue to boost the accuracy of individual approaches. 244 We also conclude from this study that it is worth revisiting the reconstruction accuracy of real protein families, 245 particularly for those widely relied-upon reconstructions [20][21][22]. The reconstruction of proteins from a single 246 alignment of small numbers of orthologs may suffer from the same or worse accuracy issues we saw in single 247 alignment approaches of our synthetic family.  [45]. 100 phylogenies were generated with each average branch length of 0.5, 0.6, 258 0.7, 0.8, 0.9, and 1.0, 600 in all. The Ensembl Compara species tree topology [46] containing 66 metazoan species was 259 used for the divergence of each ortholog set. The topology was parametrized with branch lengths corresponding to 260 species divergence times at http://www.timetree.org [47,48]. For each of 15 leaves in each random phylogeny, a copy 261 of the parametrized species tree was randomly scaled in overall height and had each individual branch randomly

293
In order to compare reconstructed paralog divergence topologies we had to modify the RF metric to 294 accommodate cases when the MRCA of an ortholog set has as descendants one or more MRCAs of other ortholog 295 sets. Such topologies are poorly formed because they require inference of additional unobservable events -loss of 296 paralogs in some lineages -in order to be reconciled with a duplication/speciation divergence history. Because the 297 offending ortholog set cannot be pruned to a leaf MRCA, the resulting topology cannot be compared to properly 298 formed topologies (e.g. the true topology) using the standard RF metric. In effect, when ortholog leaves and 299 speciation internal nodes of the offending ortholog set are pruned, the resulting topology is missing a MRCA leaf, 300 because that MRCA maps to an internal node, making that node ambiguous in its duplication vs speciation status. 301 This is problematic for RF because it affects the denominator. Nevertheless, their internal nodes representing 302 pre-duplication common ancestors of the offending ortholog set/paralog and other paralogs can match, in terms of 303 induced partition of paralogs, equivalent nodes in other topologies.

304
In the modified RF ⇤ , N represents the number of paralogs (ortholog sets) in each compared topology, not the 305 number of leaves. In addition to x and y we define z as the number of MRCA leaves missing from A but not from B 306 and z 0 as the number of MRCA leaves missing from B but not from A. The modified metric is then calculated as: