{'fraction_non_alphanumeric': 0.02595199563242147, 'fraction_numerical': 0.009768557334217248, 'mean_word_length': 5.001404165691552, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Using Genotype Abundance to Improve Phylogenetic Inference

Abstract Modern biological techniques enable very dense genetic sampling of unfolding evolutionary histories, and thus frequently sample some genotypes multiple times. This motivates strategies to incorporate genotype abundance information in phylogenetic inference. In this article, we synthesize a stochastic process model with standard sequence-based phylogenetic optimality, and show that tree estimation is substantially improved by doing so. Our method is validated with extensive simulations and an experimental single-cell lineage tracing study of germinal center B cell receptor affinity maturation.


Introduction
Although phylogenetic inference methods were originally designed to elucidate the relationships between groups of organisms separated by eons of diversification, the last several decades have seen new phylogenetic methods for populations that are evolving on the timescale of experimental sampling (Drummond et al. 2003). This development is being spurred by new experimental techniques that enable deep sequencing at single-cell resolution, some of which enable quantification of original abundance. For bulk sequencing, random barcodes can be used to quantify PCR template abundance (Jabara et al. 2011;Kivioja et al. 2012;Brodin et al. 2015). More recently, cell isolation (Shapiro et al. 2013) or combinatorial techniques (Cusanovich et al. 2015;Howie et al. 2015;DeWitt et al. 2016) have provided sequence data at single-cell resolution. With such data, a given unique genotype-among many in the data-is represented in a measured number of cells. The abundance of a genotype can be read out as the number of cells bearing that genotype. Here, we demonstrate that incorporating genotype abundance improves phylogenetic inference for densely sampled evolutionary processes in which it is common to sample genotypes more than once.
We are motivated by the setting of B cell development in germinal centers. B cells are the cells that make antibodies, or more generally immunoglobulins. Immunoglobulins are encoded by genes that undergo a stage of rapid Darwinian mutation and selection called affinity maturation . During affinity maturation, immunoglobulin is in its membrane-bound form, known as the B cell receptor (BCR). The biological function of this process is to develop BCRs with high-affinity for a pathogen-associated antigen molecule, and later excrete large quantities of the associated antibody.
This affinity maturation process occurs in specialized sites called germinal centers in lymph nodes, which have specific cellular organization to enable B cells to compete among each other to bind a specific antigen (proliferating more readily if they do) while mutating their BCRs via a mechanism called somatic hypermutation (SHM). Using microdissection, researchers can extract germinal centers from model animals and sequence the genes encoding their BCR directly (Kuraoka et al. 2016;Tas et al. 2016). Lymph node samples are also available through autopsy (Stern et al. 2014) or fine needle aspirates from living subjects (Havenar-Daughton et al. 2016). Such samples provide a remarkable perspective on an ongoing evolutionary process.
Indeed, these samples contain a population of cells with BCRs that differentiated via SHM at various times and have various cellular abundances. Because the natural selection process in germinal centers appears permissive to a variety of BCR-antigen affinities (Kuraoka et al. 2016;Tas et al. 2016), earlier-appearing BCRs are present at the same time as laterappearing BCRs. The collection of descendants from a single founder cell in this process naturally form a phylogenetic tree. However, it is a tree in which each genotype is associated with a given abundance, and such that older ancestral genotypes are present along with more recent appearances. Reconstruction of phylogenetic trees from BCR data may benefit from methods designed to account for these distinctive features.
Standard sequence-based methods for inferring phylogenies fall into several classes according to their optimality criteria. Maximum likelihood methods posit a probabilistic substitution model on a phylogeny and find the tree that maximizes the probability of the observed data under this model (Felsenstein 1973(Felsenstein , 1981(Felsenstein , 2003. Bayesian methods augment likelihood with a prior distribution over trees, branch lengths, and substitution model parameters, and approximate the posterior distribution of all the above variables by Markov chain Monte Carlo (MCMC) (Huelsenbeck et al. 2001;Drummond and Bouckaert 2015). Maximum parsimony methods use combinatorial optimization to find the tree that minimizes the number of evolutionary events (Eck and Dayhoff 1966;Kluge and Farris 1969;Fitch 1971). Parsimony methods often result in degenerate inference, in which multiple trees achieve the same minimal number of events (i.e., maximum parsimony) (Maddison 1991). Additional approaches include distance matrix methods, which summarize the data by the distances between sequence pairs, and phylogenetic invariants, which select topologies based on the value of polynomials calculated on character state pattern frequencies. None of the above methods incorporate genotype abundance information, and it is standard for data with duplicated genotypes to be reduced to a list of deduplicated unique genotypes before a phylogeny is inferred.
In this article, we show that genotype abundance is a rich source of information that can be productively integrated into phylogenetic inference, and we provide an opensource implementation to do so. We incorporate abundance via a stochastic branching process with infinitely many types for which likelihoods are tractable, and show that it can be used to resolve degeneracy in parsimony-based optimality. We first validate the procedure against simulations of germinal center BCR diversification. We also empirically validate our method using an experimental lineage tracing approach combining multiphoton microscopy and single-cell BCR sequencing, allowing us to study individual germinal center B cell lineages from brainbow mice. Beyond the setting of BCR development, we foresee direct application to tumor phylogenetics in single-cell studies of cancer evolution (reviewed by Schwartz and Schaffer 2017), and single-cell implementations of lineage tracing based on genome editing technology (McKenna et al. 2016).

Genotype-Collapsed Trees
Given sequence data obtained from a diversifying cellular lineage tree ( fig. 1a), our goal is to infer the genotype-collapsed tree (GCtree) defining the lineage of distinct genotypes and their observed abundances ( fig. 1b). The GCtree is constructed from the lineage tree by collapsing subtrees composed of cells with identical genotype to a single node annotated with its final cellular abundance. Our data consist of the genotypes sampled at least once in the GCtree, along with their associated abundances. Under the infinite types assumption that every mutant daughter generates a novel genotype, each genotype can be identified with one subtree in the original lineage tree. We are not claiming any originality in the GCtree definition, but it is useful to have a word for this object.
We note that, unlike standard phylogenetic trees where only leaf nodes represent observed genotypes, GCtree internal nodes represent observed genotypes if they are annotated with nonzero abundance. Although not leaves per se in the GCtree, a nonzero abundance represents a clonal sublineage that resulted in a nonzero number of leaves of that genotype in the lineage tree. A node in the GCtree, along with its descending edges, summarizes the lineage outcome for a given genotype as its mutant offspring clades and the number The corresponding genotype-collapsed tree (GCtree) describes the descent of distinct genotypes, and is our inferential goal. (c) Genotype abundance informs topology inference. Two hypothetical GCtrees, equally optimal with respect to the sequence data, propose two possible parents of the green genotype-the gray and yellow genotypes (the yellow genotype was not sampled and thus has a small circle with no number inside). Intuitively, the abundance information indicates that the tree on the left is preferable because the more abundant parent is more likely to have generated mutant descendants. of its clonal leaves. Because this summary does not completely specify the genotype's clonal lineage structure ( fig. 2c), several branching structures may be consistent with a given node, and we have no information with which to distinguish between the various lineage trees consistent with a GCtree. Hence, our goal is to infer the GCtree topology.
Parsimony with a Prior BCR sequence data from a germinal center sample have the following characteristics from the perspective of phylogenetics: genotypes have abundances, there is a limited amount of mutation between genotypes, and ancestral genotypes are present along with later ones. The latter two features suggest maximum parsimony as a useful tool because of the limited amount of mutation and because ancestral genotypes can be assigned to internal nodes of the tree (although recent Bayesian methods can do such assignment as well; Gavryushkina et al. 2014Gavryushkina et al. , 2017. For these reasons, parsimony has been used extensively in B cell sequence analysis (Barak et al. 2008;Stern et al. 2014). Because having many duplicate sequences inhibits efficient tree space traversal, these studies have inferred trees using the unique genotypes (BCR sequences). This ignores the varying cellular abundances of the observed genotypes.
Here, we wish to use a branching process model to rank trees that are equally optimal according to sequence-level optimality criteria. Indeed, maximum parsimony often results in degenerate inference: there are many trees that are maximally optimal (Maddison 1991). We refer to these trees as a parsimony forest. In later sections, we show, using in silico and empirical data, that parsimony degeneracy is common and often severe for BCR sequencing data, and that parsimony forests exhibit substantial variation in phylogenetic accuracy.
It is common practice to arbitrarily select one tree in the parsimony forest at random, without regard for this variability in inference accuracy. Instead, we will rank trees in the parsimony forest with an auxiliary likelihood that incorporates abundance information, thereby resolving this degeneracy.
Genotype abundance is an additional source of information for phylogenetics, using the simple intuition that more abundant genotypes are more likely to have more mutant descendant genotypes. This intuition makes sense because relative sample abundance is a reasonable estimator of relative total historical abundance, and total historical abundance is related to the number of mutant offspring-that is, genotypes with larger abundance are likely to have more mutant descendant genotypes simply because there are more individuals available to mutate. The number of mutant offspring genotypes is in turn related to the number of surviving mutant offspring sampled. Thus, given two equally parsimonious trees, this intuition would prefer the tree that has more mutant descendants of a frequently observed node ( fig. 1c). We formalize this intuition using a stochastic process model for the phylogenetic development of germinal centers, and integrate this model with sequence-based tree optimality via empirical Bayes.
In this stochastic process model, a GCtree node i has a random number T i 2 N of mutant children (i.e., descending edges) and a random abundance A i 2 N. We will index nodes in a "level order" as follows, which is well defined given an embedding of the tree into the plane. Index 1 refers to the root node, and 2 through 1 þ T 1 refer to the children of the root node. The level-order continues in order through all tree nodes of the same level before nodes at the next level. Adopting this level-ordering convention, a GCtree containing N nodes is specified by integer-valued random vectors giving FIG. 2. Modeling sequences equipped with abundances. (a) Both genotype sequence data G and genotype abundance data A inform tree topology T. As illustrated in this probabilistic graphical model, we assume independence between G and A conditioned on T rather than a fully joint model of G, A, and T. This facilitates using standard sequence-based phylogenetic optimality for G, augmented with a branching process (with parameters h) for A. (b) For the binary infinite-type Galton-Watson process, h ¼ ðp; qÞ. Four possible branching events characterize the offspring distribution common to all nodes. A node may bifurcate (with probability p) or terminate, and upon bifurcating its descendants each may be a mutant (with probability q). (c) A GCtree node specifies a genotype's clonal leaf count and number of descendant genotypes, but not lineage details. The likelihood of a GCtree node marginalizes over consistent lineage branching outcomes. (d) GCtree likelihood factorizes into the product of likelihoods for each genotype. the (planar) topology T ¼ ðT 1 ; . . . ; T N Þ, and abundances A ¼ ðA 1 ; . . . ; A N Þ. We also have the observed genotype sequences associated with each node G ¼ ðG 1 ; . . . ; G N Þ.
A complete diversification model would give a joint distribution on T, G, and A. As an approximation to such a model, facilitating use of existing sequence-based optimality methods, we propose a generative model containing conditional independences as follows ( fig. 2a). First, we model abundances A and tree topology T as being drawn from a branching process likelihood, conditioned on parameters h (characterizing birth, death, and mutation rates in the underlying lineage tree): PðA; TjhÞ. This stochastic process likelihood will capture the intuition (described earlier) that more abundant genotypes are likely to have more mutant descendant genotypes. Next, we assume that genotype sequences G are generated by a mutation model conditioned on the fixed tree T, independent of A. This sequence-based optimality is captured by a distribution over G dependent only on T: PðGjTÞ. The lack of direct dependence of G on A constitutes an approximation to a more realistic sequence-valued branching process. However, this formulation has the advantage that it allows us to leverage standard sequence-based phylogenetic optimality in the specification of PðGjTÞ. In a later section (In silico validation), we validate this approximation with simulations that do not assume this conditional independence.
In an empirical Bayes treatment (see Materials and Methods for details), a maximum likelihood estimate for the branching process parameters,ĥ, can be obtained by marginalizing T, and this in turn can be used to approximate a posterior over T conditioned on the data G and A (as well aŝ h). Using parsimony as our sequence-based optimality, one can rank trees in the parsimony forest (denoted T G ) according to the GCtree likelihood. We encode the parsimony criteria in PðGjTÞ by assigning uniform weight to the trees in T G , and zero to the other trees. This gives the following approximate maximum a posteriori tree: (1) where the point estimateĥ is given bŷ Next we turn to explicitly defining the GCtree likelihood PðA; TjhÞ.

A Stochastic Process Model of Abundance
To compute likelihoods PðA; TjhÞ for GCtrees ( fig. 1b), we model the lineage tree ( fig. 1a) as a subcritical infinite-type binary Galton-Watson (branching) process (Harris 2002) in which extinct leaf nodes correspond to observed cells. All mutations in an infinite-type process result in a novel genotype, embodying the assumption that each genotype can be identified with one subtree. Subcriticality ensures that the branching process terminates in finite time, so an explicit sampling time is not needed. The process is initiated with a single cell (a naive germinal center B cell before affinity maturation ensues), and runs to eventual extinction. This model is highly idealized and unable to capture many biological realisms of B cell affinity maturation and the sampling process. However, as we show in our validations, it is useful as a minimal model for leveraging genotype abundance information in a tractable likelihood. The offspring distribution for our process, governing reproduction and mutation for all lineage tree nodes at all time steps, is specified by two parameters: the binary branching probability p, and the mutation probability q. Because the offspring distribution is independent of type, subcriticality simply requires that the expected number of offspring of any node is <1, in this case equivalent to p < 0.5. In this case a "mutation" is an event that causes the evolving lineage to change to a novel genotype (under the infinite-types assumption). Thus the corresponding offspring distribution supports four distinct branching events ( fig. 2b). Letting C and M denote the (random) number of clonal and mutant offspring of any given node in the lineage tree, respectively, the offspring distribution is We can compute the likelihood of a hypothetical binary lineage tree simply by evaluating equation (3) at each node in the tree and multiplying the results. The likelihood for a GCtree is then given by summing over all possible binary lineage trees that are consistent with that GCtree (i.e., that give the same GCtree when collapsing by genotype), thus marginalizing out the details of intragenotype branching events that give rise to the same abundance. Here, we show how to calculate the GCtree likelihood directly for the simple offspring distribution (eq. 3). Other work (Bertoin 2009) has described how to calculate statistics of the infinite-type branching process with a general subcritical offspring distribution.
First consider the likelihood for an individual node in the GCtree, say the root node, in the context of the branching process described earlier. A GCtree node i is specified by its abundance A i and the number of edges descending from it T i (both random variables). There are, in general, multiple distinct branching process realizations for genotype i that result in A i ¼ a clonal leaves and T i ¼ s mutations off the genotype i lineage subtree ( fig. 2c). Determining the likelihood of node i in the GCtree under this process, which we denote by f as ðp; qÞ ¼ PðA i ¼ a; T i ¼ sjh ¼ ðp; qÞÞ, requires marginalizing over all such genotype lineage subtrees. In Materials and Methods, we derive a recurrence for f as ðp; qÞ by marginalizing over the outcome of the branching event at the root of the lineage subtree for genotype i, and show that the GCtree node likelihood f as ðp; qÞ can be computed by dynamic programming.
A complete GCtree containing N nodes is specified by level-ordering the nodes as described earlier T ¼ ðT 1 ; . . . ; T N Þ; A ¼ ðA 1 ; . . . ; A N Þ. Because the same offspring distribution generates the lineage branching of each genotype subtree, the same recurrence can be applied to all GCtree nodes. Specifically, we show in Materials and Methods that the joint distribution over all nodes in a GCtree factorizes by genotype ( fig. 2d).
Using dynamic programming and factorization by genotype, the computational complexity of the GCtree likelihood is OðmaxðAÞmaxðTÞ þ NÞ. Ranking parsimony trees with GCtree requires a polynomial increase in runtime compared with finding the parsimony forest, which is itself NP-hard (Foulds and Graham 1982). Supplementary figure S1, Supplementary Material online, depicts runtime from simulations of various size, and shows that, in practice, this increased runtime is negligible.
A computational implementation of the inference method above is available at http://github.com/matsengrp/gctree; last accessed February 23, 2018. The GCtree inference subprogram accepts sequence data in FASTA or PHYLIP format, determines a parsimony forest from the unique sequences using the dnapars program from the PHYLIP package (Felsenstein 2005), determines the genotype-collapsed form of these trees and outputs tree visualizations using the ETE package (Huerta-Cepas et al. 2016), and ranks them according to their GCtree likelihood using the sequence abundances. Bootstrap analysis is also implemented, providing confidence values of each split in the maximum likelihood GCtree. The GCtree maximizing the branching process likelihood (with optional bootstrap support) is the inference result. Next we show that resolving parsimony degeneracy using GCtree substantially increases both accuracy and precision of phylogenetic inference.

In Silico Validation
To explore the accuracy and robustness of GCtree inference, we developed a simulation subprogram to generate random lineages starting with a naive BCR sequence. For simulated lineages, true trees can be compared against those inferred with the GCtree inference subprogram. The stochastic process model used in GCtree inference is intended as a minimal model (in terms of biological realism) that captures the intuition that genotype abundance is relevant to phylogenetic reconstruction. Experimental data need not obey our simplifying assumptions, thus, we set out to test GCtree's robustness to deviations of the data generating process from the inferential model.
A simulation process was implemented that includes biological realisms of B cells undergoing SHM (and violates inferential assumptions). These realisms of simulationdetailed in Materials and Methods-include: branching process multifurcations (controlled by a parameter k, the expected number of children of a node in the cell lineage tree), sequence context sensitive mutations (Dunn-Walters et al. 1998;Spencer and Dunn-Walters 2005) (with a baselineline mutation rate k 0 , and a context-specific mutational model with 5mer mutabilities taken from Yaari et al. 2013), explicit sampling time (t, or N representing the number of cells desired in the sampled generation), incomplete sampling (the number of cells to sample n N), and repeated genotypes allowed (deviation from the infinite-type assumption). This constitutes a more challenging validation than simply simulating under the same assumptions that had been invoked for tractability of the inferential framework.
Our in silico validation workflow is demonstrated in figure 3a for a small simulation that resulted in a parsimony forest with just two equally parsimonious trees. The output of the simulation software consists of FASTA data (sequences and their abundances), visualizations of the lineage tree and its GCtree equivalent, and a file containing the true GCtree structure. The GCtree inference subprogram can then be run on the FASTA data, and the resulting inferred GCtree compared with the true GCtree (in this case they were identical). To calibrate simulation parameters, we defined summary statistics on sequence data with abundance information, and tuned parameters to produce data similar to experimental BCR sequencing data under these statistics (see Materials and Methods).
Our validation shows that using abundance information via a branching process likelihood can substantially improve inference results ( fig. 3b). For each simulation, we ranked otherwise degenerately optimal parsimony trees using GCtree. For each parsimony forest, we compared the GCtrees in the forest to the true GCtree for that simulation using the Robinson-Foulds (RF) distance (Robinson and Foulds 1981) as a measure of tree reconstruction accuracy. The maximum likelihood GCtree tends to be closer to the true tree than other equally parsimonious trees, which vary widely in accuracy, showing that GCtree is able to leverage abundance data to resolve parsimony degeneracy and improve the accuracy of tree reconstruction in this simulation regime.

Empirical Validation
We next performed a biological validation by investigating if GCtree improves inference according to biological criteria using real germinal center BCR sequence data. The BCR is a heterodimer encoded by the immunoglobulin heavy chain (IgH) and immunoglobulin light chain (IgL) loci. Both loci undergo V(D)J recombination, and then evolve in tandem during affinity maturation. By obtaining matched sequences from both loci using single-cell isolation, we have two independent data sets to inform the same phylogeny of distinct cells (each of which is associated with a single IgH sequence and single IgL sequence). Performing separate and independent IgH and IgL tree inference, we can then validate GCtree by comparing the inferred IgH tree to the inferred IgL tree. If the GCtree likelihood (eq. 4) meaningfully ranks equally parsimonious trees, then the Using Genotype Abundance to Improve Phylogenetic Inference . doi:10.1093/molbev/msy020 two MLE trees (IgH and IgL) would be expected to be more correct reconstructions than the other parsimony trees. Thus, we are to expect that the two MLE trees are more similar to each other (in terms of the lineage of distinct cells) than other pairs of IgH and IgL parsimony trees (which, if they are more distorted phylogenies, should show less concordance in the partitioning of the distinct cells). Conversely, if the GCtree likelihood is not meaningfully ranking trees, we expect that the MLE IgH and IgL trees will not be significantly closer to each other than other pairs of IgH and IgL parsimony trees.
We used data from a previously reported experiment in which multiphoton microscopy and BCR sequencing were combined to resolve individual germinal center B cell lineages from mouse lymph nodes 20 days after subcutaneous immunization with alum-adsorbed chicken gamma globulin (Tas et al. 2016) (see Materials and Methods). Brainbow mice were used for multicolor cell fate mapping, enabling B cells and their progeny to be permanently tagged with different fluorescent proteins. In situ photo-activation followed by fluorescence-activated cell sorting yielded B cells from a color-dominant germinal center ( fig. 4a, left). BCR sequences were obtained for 48 cells in this lineage by single-cell mRNA sequencing of the IgH and IgL loci, resulting in 32 distinct IgH and 26 distinct IgL genotypes due to SHM mutations acquired through affinity maturation. The unmutated naive IgH and IgL V(D)J rearranged sequences (not observed) were inferred with partis using each set of 48 sequences (IgH and IgL) as a clonal family using germline genetic information Matsen 2016a, 2016b). These naive sequences were used as outgroups for rooting parsimony trees.
GCtree results are depicted in figure 4b. Parsimony analysis resulted in degeneracy for both loci, with 13 equally FIG. 3. In silico validation of GCtree inference. (a) Demonstrating the simulation-inference-validation workflow, a small simulation resulted in two equally maximally parsimonious trees, and the one inferred using GCtree was correct. The initial sequence was a naive BCR V gene from the experimental data described in Materials and Methods. Branch lengths in the cell lineage tree (left) correspond to simulation time steps, while those in collapsed trees correspond to sequence edit distance. (b) About 100 simulations were performed with parameters calibrated using the BCR sequencing data and summary statistics described in Materials and Methods. Of 100 simulations, 66 resulted in parsimony degeneracy, with an average degeneracy of 12 and a maximum degeneracy of 124. For each of these 66, we show the distribution of Robinson-Foulds (RF) distance of trees in the parsimony forest to the true tree. "RF" denotes a modified Robinson-Foulds distance: since nonzero abundance internal nodes in GCtrees represent observed taxa, RF distance was computed as if all such nodes had an additional descendant leaf representing that taxon. GCtree MLEs (red) tend to be better reconstructions of the true tree than other parsimony trees (gray boxes). Four simulations resulted in a tie for the GCtree MLE, and the two tied trees in these cases are both displayed in red. Aggregated data across all simulations are depicted on the right, clearly indicating superior reconstructions from GCtree. parsimonious trees for IgH, and 9 for IgL. Empirical Bayes point estimation according to equation (2) yieldedp ¼ 0:495; q ¼ 0:388 (IgH) andp ¼ 0:495;q ¼ 0:304 (IgL). GCtree likelihoods (eq. 4) were computed to rank the equally parsimonious trees, and the MLE trees are shown with support values among 100 bootstrap samples (see Materials and Methods). Because the binary Galton-Watson process assigns probability zero to a GCtree node with frequency zero and one mutant descendant, the unobserved naive root node (which had one descendant after rerooting and collapsing identical genotypes in all parsimony trees) was given a unit pseudocount. These sequences were analyzed with partis Matsen 2016a, 2016b) to infer naive (preaffinity-maturation) ancestor sequences using germline genetic information, and trees were inferred with GCtree. (b) GCtree inference was performed separately for IgH and IgL loci, resulting in parsimony degeneracies of 13 and 9, respectively. Maximum likelihood GCtrees for each locus are indicated in red and the GCtrees with annotated abundance are shown. Roots are labeled with the gene annotations of the naive state inferred using partis. Small unnumbered nodes indicate inferred unobserved ancestral genotypes. Numbered edges indicate support in 100 bootstrap samples. (c) All possible pairings of IgH and IgL parsimony trees were compared in terms of the Robinson-Foulds distance between the IgH and IgL trees, labeled by cell identity. IgH and IgL parsimony trees are ordered by GCtree likelihood rank in columns and rows, respectively. Grid values show RF distance between each IgH/IgL pair. MLE trees result in more consistent cell lineage reconstructions between IgH and IgL (smaller RF values). (d) For each locus, distributions of bootstrap support values are shown for the tree inferred by GCtree and for a majority rule consensus tree of all trees in the parsimony forest. The latter contain more partitions with very low support. (e) Using additional data from a second germinal center from the same lymph node that had the same naive BCR sequence, GCtree correctly resolves the two germinal centers as distinct clades (as did other lower ranked parsimony trees).
We then compared the concordance between pairs of heavy and light trees. Since both IgH and IgL loci have been recorded from the same set of 48 cells, the units of cell abundance in an IgH GCtree map to the units of cell abundance from an IgL GCtree (i.e., each cell identity among the 48 is associated with an IgH genotype and an IgL genotype). We can then consider the consistency of a given IgH tree and a given IgL tree in terms of the lineage of the 48 cell identities. For each possible pairing of an IgH parsimony tree with a IgL parsimony tree, we computed the RF distance (Robinson and Foulds 1981) between the two trees using the cell identities (rather than the genotype sequences) to define splits. We observed that the GCtree MLE based on IgH sequences and GCtree MLE based on IgL sequences form the most concordant pair among all pairs of parsimony trees ( fig. 4c). Moreover, pairs of parsimony trees that contained at least one GCtree MLE tree ranked consistently higher in terms of their similarity.
We assessed confidence in GCtree partitions by comparing bootstrap support values in GCtree trees to those from the majority-rule consensus parsimony trees made using the consense program from the PHYLIP package (Felsenstein 2005). We observed the latter contained an excess of very low confidence partitions ( fig. 4d and supplementary fig. S4, Supplementary Material online). These results demonstrate that parsimony reconstructions for real BCR data sets suffer from degeneracy, and that GCtree likelihood can correctly resolve this degeneracy by incorporating abundance information ignored by previously published methods.
Finally, using data collected from a second germinal center from the same lymph node, we tested GCtree's ability to correctly group cells from each germinal center into separate clades when run on combined data from both germinal centers. The two germinal center sequence data sets appeared to have the same naive BCR sequence (IgH and IgL), indicating they were both seeded from the same B cell lineage. Concatenating the IgH and IgL sequences for each cell in each germinal center, we used GCtree to infer a single tree for all cells from both germinal centers ( fig. 4e and supplementary fig. S5, Supplementary Material online). GCtree correctly resolved the two germinal centers as distinct clades (we note that all the parsimony trees had this feature, regardless of likelihood rank). This demonstrates the phylogenetic resolvability of germinal centers with the same naive BCR diversifying under selection for the same antigen specificity.

Discussion
We have shown that genotype abundance information can be productively incorporated in phylogenetic inference. By augmenting standard sequence-based phylogenetic optimality with a stochastic process likelihood, we were able to implement abundance-aware inference as a processing step downstream of results from an existing and widely used parsimony tree inference tool. We have shown that our method-implemented in the publicly available GCtree package-is useful for inferring B cell receptor affinity maturation lineages. Although branching processes have been used previously to infer parameters of BCR evolution (Kleinstein et al. 2003;Magori-Cohen et al. 2006) and construct SHM lineage trees from error-prone bulk sequencing reads (Sok et al. 2013), to our knowledge, we are the first to use branching processes to sharpen phylogenetic inference for BCRs sequenced at single-cell resolution from germinal centers.
We believe GCtree will find use in other settings where sequence data from dense quantitative sampling of diversifying loci are available. Studies of cancer evolution are increasingly performed with single-cell resolved sequencing, however most tumor phylogenetics approaches use standard phylogenetic methods (reviewed by Schwartz and Schaffer 2017) that do not model genotype abundance. Exceptions include OncoNEM (Ross and Markowetz 2016) and SCITE (Jahn et al. 2016), both of which leverage single-cell data for tumor phylogenetic inference that is robust to genotyping errors and missing data, but do not aim to capture the intuition that genotype abundance and the number of direct mutant descendants are related. Single-cell implementations of lineage tracing based on genome editing technology (McKenna et al. 2016) may also benefit from reconstruction methods that model the abundance of observed editing target states, since cell types may vary widely in rates of proliferation.
Using parsimony as our sequence-based optimality resulted in particularly simple results, as the tree space necessary to explore is exactly the degenerate parsimony forest. However, our empirical Bayes formulation is agnostic to the particular choice of sequence-based optimality, so in the future, we envision augmenting likelihood-based sequence optimality. This will require more computationally expensive tree space search and sampling schemes.
In contrast to GCtree, a fully Bayesian approach to incorporate genotype abundance could use the full set of sequences (without deduplication) in a Bayesian phylogenetics package-such as BEAST (Drummond and Bouckaert 2015)-with a birth-death process prior. This would not enforce the infinite-type assumption, so a set of identical sequences could be placed in disjoint subtrees. However, such an approach will not scale well with many identical sequences: trees that only differ by exchange of identical sequences will create islands of constant posterior in tree space. Methods do not currently exist for tree space traversal that avoids moves within such islands. Even if such methods existed, they would need to be combined with algorithms to infer trees with sampled ancestors (Gavryushkina et al. 2014(Gavryushkina et al. , 2017 as well as multifurcations (Lewis et al. 2005(Lewis et al. , 2015; even just this combination is not currently available. Although our methods can be applied to other sequencebased optimality functions besides parsimony, it is important to recognize that GCtree (and indeed any tree inference procedure that deduplicates repeated sequences) contains an inherent weak parsimony assumption: that each unique genotype arose from mutation just once in the lineage and therefore corresponds to a single subtree in the lineage tree, and thus a single node in the GCtree. Thus it is important to continue to assess the impact of this weak parsimony assumption with simulation that does not make this assumption, as done here.
The GCtree framework can also be extended to nonneutral models. For example, one could consider a model in which each genotype obtains a random fitness encoded by branching process parameters h that are fixed within a given genotype but randomly drawn by the genotype founder cell upon mutation from its parent. This will likely necessitate modeling genotype birth time explicitly, rather than restricting to extinct subcritical processes, since a genotype with small abundance may be a result of low fitness or just young age. One might also consider extending the offspring distribution to separately model synonymous and nonsynonymous mutations. Synonymous mutations do not change fitness, while nonsynonymous mutations change fitness as described earlier.
Another direction of extension is to incorporate mutation models specialized to the case of BCR evolution, such as the S5F model (Yaari et al. 2013) used in our simulation study.

An Empirical Bayes Framework for Incorporating Genotype Abundance in Phylogenetic Optimality
Here, we more fully develop the empirical Bayes perspective on our estimator for the model depicted in figure 2a. This graphical model implies the factorization PðG; A; T; hÞ ¼ PðGjTÞPðA; TjhÞPðhÞ: A hierarchical Bayes treatment would assign a prior PðhÞ (such as uniform over the unit square for the model h ¼ ðp; qÞ) and compute the posterior over trees conditioned on the data, marginalizing over h: Rather than attempting this integral over PðA; TjhÞ, each evaluation of which requires dynamic programming, we first seek a maximum likelihood estimate for h marginalizing T: Using this point estimate, an approximate posterior over trees is PðTjG; A;ĥÞ / PðGjTÞPðA; TjĥÞ: This formulation embodies an optimality over trees conditioned on both genotype sequence data G and genotype abundance data A. Evaluation ofĥ with equation (6) in general requires summation over the space of all trees consistent with the data. A simple application of this formalism is to augment parsimony-based tree optimality with abundance data. Let T G denote the degenerate set of maximally parsimonious trees given G (each of which has the same total genotype sequence distance over its edges). Encode parsimony optimality as a PðGjTÞ assigning uniform weight to each tree in T G , and zero elsewhere. In this case, equation (2) becomeŝ and equation (7) becomes PðTjG; A;ĥÞ / PðA; TjĥÞ; t 2 T g 0; t 6 2 T g : ( With equation (9), we have a framework using abundance information to distinguish among the otherwise equally optimal trees presented by a parsimony analysis. In our application, we use a subcritical infinite-type binary Galton-Watson branching process model for the lineage tree, and describe a recursion for computing GCtree likelihoods PðA; TjĥÞ by dynamic programming to marginalize over compatible lineage trees.

Dynamic Programming to Marginalize Lineage Tree Structure
We derive a recurrence for f as ðp; qÞ ¼ PðA i ¼ a; T i ¼ sjh ¼ ðp; qÞÞ by marginalizing over the outcome {C, M} of the branching event at the root of the lineage subtree for genotype i (the first cell of type i). We will use that a and s are the sum over two iid processes for the left and right clonal branches. We temporarily suppress the parameters h ¼ ðp; qÞ, writing f as for notational compactness. In the case fC ¼ 2; M ¼ 0g, f a 0 s 0 f aÀa 0 ;sÀs 0 : As this is the convolution of f as with itself, we denote it as f Ã2 as . Marginalizing over all outcomes {C, M}, we have f Ã2 a0 pð1 À qÞ 2 a > 1; s ¼ 0; f a;sÀ1 2pqð1 À qÞ þ f Ã2 as pð1 À qÞ 2 otherwise; where d ÁÁ denotes the Kronecker delta function. In light of the first case, the convolutional square may be written as showing that there are no terms containing f as on the RHS of equation (11). The GCtree node likelihood f as is thus amenable to computation by straightforward dynamic programming.

The GCtree Likelihood Factorizes by Genotype
We argue that the joint distribution over all nodes in a GCtree factorizes by genotype ( fig. 2d): Since s 1 is the number of children of node 1 (the root node), the children of the root node are indexed in level order by 2; . . . ; 1 þ s 1 . Let K i denote the set of indices of the nodes of the subtree rooted at node i, so K 2 ; . . . ; K 1þs 1 refer to sister subtrees rooted on each of the s 1 children of the root. Using the definition of conditional probability, and since sister subtrees are independent, we have Pðfða j ; s j Þ : j 2 K i gÞ; where random variable notation has been dropped for notational compactness. Now, within each subtree factor, we may reindex in level order (that is, level order in that subtree) starting from 1. We then pull out factors f a 2 s 2 ; . . . ; f a 1þs 1 s 1þs 1 corresponding to the root nodes of the sister subtrees (children of the original root). We obtain equation (12) by applying this logic recursively. Restoring the offspring distribution parameters, we recognize this as the distribution needed in equations (1) and (2) to rank trees in a parsimony forest: where f a i s i ðp; qÞ is computed by dynamic programming using equation (11). Numerical validation of the GCtree likelihood is summarized in supplementary figure S3, Supplementary Material online, using 10,000 Galton-Watson process simulations at each of several parameter values. The likelihood accurately recapitulates tree frequencies, and simulation parameters are recoverable by numerical maximum likelihood estimation.

Simulation Details
To provide for a more challenging in silico validation study, several biological realisms were built into our simulation that defied simplifying assumptions in the GCtree inference methodology.

Arbitrary Offspring Distribution
The recursion (eq. 11) used to compute GCtree likelihood components specifies a binary branching process, and such an approach would in general require an offspring distribution with bounded support on the natural numbers. Our simulation implements an arbitrary offspring distribution with no explicit bounding. In the results that follow, we used a Poisson distribution with parameter k for the expected number of offspring of each node in the lineage tree.

Context Sensitive Mutation
To generate mutant offspring, all offspring sequences (drawn from a Poisson as described earlier) were subjected to a sequence-dependent mutation process. The SHM process is known to introduce mutations in a sequence contextdependent manner, with certain hot-spot and cold-spot motifs (Dunn-Walters et al. 1998;Spencer and Dunn-Walters 2005). We used a previously published 5-mer context model S5F (Yaari et al. 2013) to compute the mutabilities l 1 , . . . ; l ' of each position 1, . . . ; ' within a sequence of length ' based on its local 5-mer context. This model also provided substitution preferences among alternative bases given the 5mer context. To compute mutabilities for beginning and ending positions without a complete 5-mer context, we averaged over missing sequence context.
Although existing code can simulate a mutational process parameterized by S5F on branches of a fixed tree with a prespecified number of mutations on each branch (Gupta et al. 2015), in our simulations, we wanted the number of mutations on the branches to be determined by the sequence mutability as it changes via mutation across the tree. For example, as an initial mutation hotspot motif acquires mutations down the tree, its mutability typically degrades as it diverges from the original motif. We defined the mutability of the sequence as a whole by the average over its positions i¼1 l i . We defined a baseline mutation expectation parameter k 0 as a simulation parameter, and the number of mutations m any given offspring sequence received was drawn from a Poisson distribution. The Poisson parameter was modulated by the sequence's mutability m $ Poisðl 0 k 0 Þ, so that more mutable sequences tended to receive more mutations. Given m > 0, the positions in the sequence to apply mutations were chosen sequentially as follows. A site j to apply the first mutation was drawn from a categorical distribution using the site-wise mutabilities to define relative probability of choosing each site j $ Catðl 1 ; . . . ; l ' Þ. We mutated the site using a categorical distribution over the three alternative bases parameterized by the substitution preferences defined by the site's context. We then updated mutabilities l 0 and l 1 , . . . ; l ' as necessary to account for contexts that had been altered by the mutation. This process was repeated m times.
Since the mutability of each node in the lineage tree will depend on the mutation outcome of its parent, the GCtree likelihood components will not factorize by genotype. Because the probability of mutation is sequence-dependent, the topology of the GCtree will be sequence-dependent. Therefore, the generative assumptions of the empirical Bayes inference do not hold in this simulation scheme, nor does the offspring distribution equivalence across lineage tree nodes specified by equation (3).

Sampling Time
Our inference model specifies a subcritical branching process run until extinction, and sampling of all terminated nodes (leaves). Our simulation more realistically assigns a discrete time of sampling parameter t (number of time steps from root), and thus does not need to constrain the offspring distribution to achieve subcriticality. At the specified time, extant nodes can be sampled, so all genotypes that terminated or mutated at a prior times are not observed. Alternatively, a parameter N specifying the desired number of simulated observed sequences may be passed, in which case the simulation runs until a time such that at least N sequences exist (unless terminated). Genotypes born at different times will be sampled under a process with different effective sampling times since birth. Thus this sampling time parameter also increases dependence between genotypes, further distancing the simulation model from the inferential model.

Incomplete Sampling
We introduce imperfect sampling efficiency with a parameter n for the number of simulated sequences that end up in the simulated sample data (FASTA), requiring n N. This violates the inferential assumption of complete sampling, and renders the true genotype abundances latent variables (which a more complete likelihood approach might aim to marginalize out).

Repeated Genotypes
Our simulation is seeded with an initial naive BCR sequence, from which randomly mutated offspring are created. Because there is no built-in restriction that the same sequence cannot arise along different branches (or mutations could be reversed), the model assumption of infinite types-such that identical sequences can be associated with a single genotype subtree-does not necessarily hold. When this assumption is violated the tree must necessarily be incorrect.

Calibrating Simulation Parameters Using Summary Statistics
We defined several summary statistics on sequences equipped with abundances which were used to calibrate simulation parameters representative of a regime similar to experimental data. We chose these statistics to reflect information relevant to tree inference, but not actually require tree inference, so as to avoid circularity. Denote g 0 2 G as the naive BCR (root genotype) and d H ðÁ; ÁÞ as the Hamming distance function between two sequences. Given simulation or experimental data G and A, we characterize the degree of mutation (from naive BCR) in the lineage by the set of Hamming distances of the observed genotypes from the naive genotype: fd H ðg; g o Þ; g 2 Gg. For a given genotype g i 2 G, we can compute its number of Hamming neighbors in the data g i ¼ jfg j 2 G : d H ðg i ; g j Þ ¼ 1gj.
A simulation is specified by parameters ðk; k 0 ; Nðor tÞ; nÞ, a mutability model (here S5F; Yaari et al. 2013), and an initial sequence. We found parameters ðk ¼ 1:5; k 0 ¼ 0:25; N ¼ 100; n ¼ 65Þ produced simulations that were comparable to experimental data under these statistics. The experimental data used for comparison, consisting of 65 total BCR V gene sequences from a single germinal center lineage, is described in the following section. Supplementary figure S2, Supplementary Material online, depicts these summary statistics for 100 simulations, compared with experimental BCR data.

Germinal Center BCR Sequencing
Germinal center B cell lineage tracing and B cell receptor sequencing was performed as previously described (Tas et al. 2016). Full length IgH and IgL sequences from lymph node 2 germinal centers 1 and 2 from this reference were used for empirical validation results, whereas V gene sequences only (which are not dependent on partis-inferred naive sequences) were used for calibrating simulation parameters.

Bootstrap Support
We computed bootstrap support values for edges on a GCtree extending the standard approach (Felsenstein 1985): we resampled columns from the alignment G 100 times with replacement, generating an inferred GCtree (maximum GCtree likelihood among equally parsimonious trees) for each. Each edge is equivalent to a bipartition of observed genotypes obtained by cutting the edge; such a bipartition is typically referred to as a split. We compute the number of bootstrapped trees that contain the same split, and annotate the edge with this number. Because resampling the alignment G can produce repeated genotypes, there can exist ambiguity about how to perform genotype collapse of a parsimony tree. We simply group genotypes in the bootstrap analysis that collapse to identical genotypes. For example, if two observed sister genotypes with resampled sequences are both identical in sequence to their mutual parent, both have a claim on collapsing into the parent. When collapsing this tree, both genotypes will be associated with this collapsed node, rather then just a single one.

Data Availability
Germinal center BCR sequence data can be found in supplementary Database S1 of Tas et al. (2016), lymph node 2 and germinal center 1.

Software Availability
The GCtree source code is available at github.com/matsengrp/gctree and accepts sequence alignments in FASTA or PHYLIP format as input. It is open-source software under the GPL v3.

Supplementary Material
Supplementary data are available at Molecular Biology and Evolution online.