{'fraction_non_alphanumeric': 0.028622487842034337, 'fraction_numerical': 0.010634557918790649, 'mean_word_length': 5.060901988636363, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 0, 'https://': 1, '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}}
A genealogical estimate of genetic relationships

The application of genetic relationships among individuals, characterized by a genetic relationship matrix (GRM), has far-reaching effects in human genetics. However, the current standard to calculate the GRM generally does not take advantage of linkage information and does not reflect the underlying genealogical history of the study sample. Here, we propose a coalescent-informed framework to infer the expected relatedness between pairs of individuals given an ancestral recombination graph (ARG) of the sample. Through extensive simulations we show that the eGRM is an unbiased estimate of latent pairwise genome-wide relatedness and is robust when computed using genealogies inferred from incomplete genetic data. As a result, the eGRM better captures the structure of a population than the canonical GRM, even when using the same genetic information. More importantly, our framework allows a principled approach to estimate the eGRM at different time depths of the ARG, thereby revealing the time-varying nature of population structure in a sample. When applied to genotyping data from a population sample from Northern and Eastern Finland, we find that clustering analysis using the eGRM reveals population structure driven by subpopulations that would not be apparent using the canonical GRM, and that temporally the population model is consistent with recent divergence and expansion. Taken together, our proposed eGRM provides a robust tree-centric estimate of relatedness with wide application to genetic studies.

temporally the population model is consistent with recent divergence and expansion. Taken together, our proposed eGRM provides a robust tree-centric estimate of relatedness with wide application to genetic studies.

Introduction
Genetic relationships among individuals, commonly characterized by a genetic relationship matrix (GRM), has fueled major advances in modern human genetics. Its applications include the detection of population structure 1,2 , adjusting for shared genetic backgrounds in genome-wide association testing 3-10 , and heritability estimation 11 . Historically, genetic relationships across pairs of individuals in a known pedigree were estimated using the expected proportion of co-inherited alleles, which neglects the variance in the distribution of alleles from meiosis [12][13][14] . The advent of high-throughput genomics has enabled estimating pairwise relationships directly from genotype data, without the need to rely on expectations determined from an inheritance model 15 .
The current standard to calculate the GRM is based on computing a weighted expectation across genotyped variants (i.e. identity-by-state or IBS) 11,12,14 . While straightforward to compute, this approach generally does not utilize linkage information between markers (though also see ref. 4,16 ). The canonical GRM is not designed to reflect the shared genealogies that connect everyone in a population and inadequately reflects the contribution of ungenotyped variation to relatedness [17][18][19] . Thus, genome-wide IBS-based relatedness is sensitive to ascertainment biases of genetic variation and only partially captures individuals' relationships compared with relatedness based on the underlying genealogies of the population. An identity-by-descent (IBD) based GRM could incorporate linkage information to infer finer-scale genetic relationships underlying the structure or demographic history of the study population. However, current bioinformatic approaches estimating shared IBD segments are subject to technical and methodological constraints that effectively limit the resolution of inferred relatedness to only the most recent branches nearing the tips of the underlying genealogies (i.e. over the last 50-100 generations) 12,13,[20][21][22][23][24] . Because of its methodological simplicity, the canonical, IBS-based, GRM continues to be the standard in statistical genetics despite its shortcomings 12,20 . Nevertheless, these shortcomings motivate the search of an approach that better captures the genealogical relatedness in a population sample.
In this study, we describe a model for pairwise relatedness using a coalescent-based framework relating everyone in a population sample. Given a coalescent tree at a locus, we define relatedness between individuals by tracing the tree backwards to a single common ancestor. The locus-specific tree provides generalized IBD information across the population sample, unlike conventional definitions of IBD that are restricted to recent branches of the tree in forms of detectable IBD segments of defined multi-generational pedigrees. The entire genealogy of a sample of individuals can be represented by a sequence of coalescent trees, encoded in a structure called the Ancestral Recombination Graph 25,26 (ARG). The ARG carries substantial linkage information as mutations on the same branch are by definition linked, and historical recombination events are encoded across the sequence of trees. In practice, the ARG is inferred through haplotypic linkage that exists in genetic data. As such, a genealogical measure of relatedness conditioned on the ARG can exploit linkage information that is commonly ignored in the canonical GRM.
Here, we propose a novel coalescent-based framework to estimate the expected genetic relatedness, or the eGRM, for pairs of individuals given the ARG of the population. Conceptually, the eGRM is based on the expected number of mutations occurring randomly on each branch of the ARG, rather than directly genotyped variants. Our framework provides two primary benefits compared with previous approaches. First, because the ARG encodes historical recombination events and the estimation of the ARG generally leverages patterns of haplotype sharing, the eGRM in practice is expected to be more robust to ungenotyped genetic variation and retains greater information of IBD relatedness among individuals than the canonical GRM. Second, and more importantly, our framework seamlessly provides insights to the time-varying nature of population structure by estimating relatedness at specific depths in the coalescent tree. To enable efficient calculation of the eGRM, our framework leverages recent computational advancements for scalable ARG inferences 27,28 , thus enabling investigation of populations in larger datasets.
We characterized the behavior of our ARG-based eGRM through extensive simulations starting from standard population genetic models. In simulations of a single, exponentially growing, population, we demonstrate that the eGRM better captures latent genome-wide relatedness compared with the canonical GRM. Importantly, we find the improved performance of eGRM is robust when performing inference using noisier ARGs inferred from a subset of common genotyped variants rather than true ARGs. It is believed that common variants are not sufficiently informative to detect recent population structure 29 . However, in simulations of a recently structured population with multiple demes, we find that principal component analysis (PCA) of the eGRM better reflects overall population structure and more accurately identifies each deme, compared to PCA of the canonical GRM. Finally, in analyses of 2,644 genotyped samples from Northern and Eastern Finland 21 , we observe that PCs derived from the eGRM reveal fine-scale structure previously not identified when using the canonical GRM. We estimate multiple partial eGRMs at multiple epochs across the history of the sample and show that time-specific patterns of population structure are qualitatively similar to simulated results of a recently structured population model, which is consistent with the known history of this region of Finland.

Method Overview: a conceptual shift in defining genetic relatedness
The eGRM, conditioned on the ARG, is conceptually different from the canonical GRM. We demonstrate this difference through a toy example on a single genealogical tree with 4 branches and 6 mutations ( Figure 1A). The canonical GRM is variant-centric and is the average of the six relatedness matrices based on each mutation. The eGRM, however, defines relatedness through tree branches that relate a pair of haplotypes. Assuming constant mutation rates across branches, the eGRM is the average of the four relatedness matrices based on each branch, weighted by their branch lengths. A single tree is shown in Figure 1A for simplicity, but the eGRM can be generalized to a sequence of trees along a chromosome by weighting each tree by its total branch length times the number of base pairs covered by each tree (Methods). In this toy example, haplotypes a and b are expected to be equally related to c in the eGRM, while in the canonical GRM b will be more closely related to c. Under the canonical GRM framework, the relative genetic distance to c is subject to the randomness and ascertainment of mutations. Instead of relying on ascertained mutations, branch lengths from the true ARG (or from the inferred ARG, reconstructed based on linkage information among nearby markers) provide an estimate of genetic relatedness that is more robust to ascertainment effects. In addition, while the eGRM is defined as a function of the genealogy, it maintains the mathematical properties of canonical GRMs (e.g., positive definiteness) as eGRM is the expectation of the canonical GRM. The eGRM is thus compatible with all downstream applications of the GRM.
To help distinguish between various eGRM estimators, here we define some useful notation. We denote the eGRM estimated conditioned on the true ARG as EK. When conditioned on an ARG inferred from genetic data using either RELATE 28 or TSINFER+TSDATE 27,30 we denote such eGRM as EKrelate and EKtsdate, respectively. We denote the canonical GRM computed using all latent genetic data as Kall and using only the observed genetic data as Kobs. In empirical data analysis, Kobs are constructed using all genotyped data passing quality controls; in simulations, Kobs are constructed using only 20% of the genetic data oversampled from the common variation of the frequency spectrum (Methods) to mimic a genotyping array. Importantly, EKrelate and EKtsdate are constructed only using the same set of the observed genetic data as Kobs.

eGRM accurately measures relatedness on a genealogical tree
To establish that the eGRM estimator better reflects genealogical relatedness compared with the canonical GRM approach, we first sought to quantify the performance of eGRM in capturing relatedness in a single tree, defined here as the TMRCA between pairs of individuals, when using the true ARG. We simulated a 1 Mb genetic region with 1,000 individuals under a single population growth model, computed EK and Kobs (see Methods; Figure 1B). Unsurprisingly, the eGRM based on the true genealogy, EK, is better correlated with TMRCA than Kobs in 97.5% of the simulations (P = 4e-252 by sign test; Figure 2A) and more accurately captures recent genetic relatedness between pairs of individuals ( Figure S1A). More importantly, eGRM constructed using genealogies inferred under RELATE or TSINFER+TSDATE on the same set of observed variants (EKrelate and EKtsdate) also showed better correlation with TMRCA than the canonical GRM in ~70% of the simulations (P < 1e-26 in all cases; Figure 2A, Figure S1B), suggesting that the eGRM is robust to noise in inferred ARGs. Our results thus demonstrate a consistent advantage of the eGRM over the canonical GRM in capturing local relatedness represented by TMRCA.
Because common variants are individually uninformative for recent relatedness, our results also suggest the eGRM framework based on predominantly common variants can provide insight for the recent part of the genealogical tree.

eGRM provides an unbiased estimate of genome-wide relatedness across realistic scenarios
While TMRCA provides an intuitive measure of the local genetic relatedness between a pair of haplotypes, the eGRM is formulated as the expectation of the latent GRM so that it adheres to the mathematical properties of a GRM necessary for many downstream statistical genetic applications [1][2][3][4][5][6][7][8][9][10][11] . Therefore, we next sought to evaluate how well eGRM is capturing the genome-wide relatedness measured by the latent GRM. Even though the mutation rate in humans is small, we reasoned that a GRM computed using all variants (Kall) in a sufficiently large genomic region and sample will precisely estimate the latent GRM by the law of large numbers. Therefore, we next sought to evaluate how well the eGRM measures genome-wide relatedness, as quantified . We observed similar performance gains when computing the eGRM using genealogies inferred by RELATE, with EKrelate attaining a highly correlated (r = 0.90 +/-0.004; Figure 2B) and approximately unbiased estimate of Kall (regression slope of 0.96, intercept 3.7e-5; Figure 2C).
We found EKtsdate demonstrated lower correlation and biased estimates of Kall ( Figure 2B, C).
Taken together, our results suggest the eGRM is an unbiased and accurate estimator of the idealized canonical GRM containing all variants.
Next, we quantified the performance of eGRM when computed using genealogies inferred from a varying proportion of observed genetic variants. We found that the correlation between Kall and EKrelate was consistently higher than the correlation between Kall and Kobs (Figure 2D left).
Moreover, for a fixed proportion of observed common SNPs (e.g., 20%; similar to SNP arrays), we observed the performance gap widened between EKrelate and Kobs as sample size increased ( Figure 2D right). Intuitively, this improvement reflects the increasing contribution from rare variants to kinship in a larger sample that would not be captured by the canonical GRM based on only variants assayed on an array. Our results imply that eGRM can in principle more effectively capture relatedness in large-scale studies.
In practice, the construction of the GRM often uses imputed variants and/or is restricted to relatively common variants after pruning of correlated variants by LD. However, in simulations we found that pruning SNPs by LD before computing the GRM (Kobs (pruned)) further decreased correlation with Kall ( Figure S2A). When using imputed variants to construct the canonical GRM (Kobs (imputed)), we found it more strongly correlated with Kall on average when compared with Kobs ( Figure S2A). However, we observed performance in this scenario depends on relatedness between individuals in the imputation reference panel with target individuals, with correlation between Kobs (imputed) and Kall decreasing with average panel-relatedness ( Figure S2B). The dependence on the availability of a closely related reference panel suggests that underrepresented populations would be at a disadvantage for genetic analysis using the canonical GRM 31 . Most importantly, across all of these scenarios, we observed our eGRM based on inferred genealogy (i.e., EKrelate) consistently exhibited better correlation with Kall than Kobs (pruned) or Kobs (imputed) (Figure S2A).

eGRM captures recent demographic events
Population structure, as the result of historical demographic processes that are encoded in the ARG, is conventionally visualized through principal components analysis (PCA) of the canonical GRM 4,32 . Because the eGRM is conditioned on the ARG encoding these historical events, we expected the eGRM to be more sensitive to population structure than the canonical GRM based on only common variants. To this end, we quantified the performance of the eGRM in capturing recent population structure through PCA. A. Motivated by recent work demonstrating that recent population structure (i.e. < 100 generations) is not well-captured by PCA computed from common SNP GRMs 29 , we simulated a stepping-stone population model with 25 demes spatially distributed in a 5x5 grid, with demes coalescing into a single deme 100 generations ago ( Figure 1B, Figure   3; see Methods). We then compared the ability of the eGRM or the canonical GRM to identify recent population structure through PCA as quantified by the separation index (SI), which measures the proportion of neighbors in multi-dimensional space that are of the same deme or cluster (Methods).
To establish a baseline, we first recapitulated previous results demonstrating that PCA based on the GRM constructed using common variants (defined as MAF ≥ 0.05; Kcommon), or the observed variants (defined as 20% of the variants, oversampling common variants; Kobs), cannot distinguish the spatial structure of the demes (SI = 0.07-0.08 for Kcommon and Kobs; Figure 3, Figure S3A). In comparison, the GRM constructed from rare variants (minor allele count = 2, 3, 4 or 5) alone (Krare; SI = 0.25) or all of the variants (Kall; SI = 0.20; Figure S3A) can better detect structure.
We repeated this PCA analysis using EKrelate computed from the same set of variants used to construct Kobs and found it to better identify recent population structure (SI = 0.22; Figure 3). We next applied a UMAP transformation to the top 10 PCs based on each of the evaluated relatedness matrices, and found cluster separation performance improved when using Krare, Kall, or EKrelate (SI = 0.80-0.82), but with little benefit when using Kcommon or Kobs (SI = 0.08-0.09; Figure   3 and Figure S3B). We also simulated a demographic history with older population split times of 200 or 500 generations. In these scenarios, EKrelate (SI = 0.36 and 0.64 when split time is 200 and 500, respectively) consistently outperformed Kobs (SI = 0.13 and 0.34, respectively) in capturing population structure (Figure S4), likely due to additional haplotypic information captured by the inferred ARG, even though Kobs also improved in its performance capturing population structure than a scenario with more recent population split because the common variants are more informative of the population structure in these simulations. Therefore, under a structured model the eGRM consistently extract more information of the population structure compared to the canonical GRM based on the same set of variants.

Partial eGRM reveals dynamic relatedness through history
By defining genome-wide relatedness as a function of coalescent trees, a major advantage of our eGRM framework is its natural generalization of relatedness constrained to a specific time window. We denote the eGRM computed from the ARG when considering only branches of a certain age as the partial eGRM (see Methods). To demonstrate the benefit of limiting relatedness calculation to certain generations, we re-analyzed our grid simulations but estimated  Figure S6). Notably, the first two PCs of EKrelate were explained in part by individuals with both parents born in the Surrendered Karelia (magenta color in Figure 4). Similar to our simulated results, when we applied UMAP to top PCs we observed improved resolution of finescale structure, with UMAP based on EKrelate continues to be more informative of the fine-scale structure within Finland than that based on the canonical GRM, regardless of the number of PCs included in UMAP ( Figure 4B).
Next, to shed light on historical migration and population movements for FinMetSeq data, we computed and analyzed the partial eGRM considering only branches for the past 0-100 generations ( Figure S5C). PCA of the partial eGRM suggests that recent structure in Northern and Eastern Finland is mainly driven by individuals from Lapland (colored red in Figure 4; the Northernmost part of Finland and home to the indigenous Finno-Urgic people, Sami), Surrendered Karelia (magenta), and Turku-Pori and Vaasa (light green and light blue, sharing major port borders with Sweden). Computing a partial eGRM further in the past exhibited patterns more similar to those found in the canonical GRM ( Figure S5C). Qualitatively, the pattern of the partial eGRM at varying time depth and its correlation with a fixed canonical GRM are more reminiscent of the pattern observed in the grid structure simulation (Figure S5B) than the pattern in a single homogenous population (Figure S5A). Together these findings further support previous claims that common variants are enriched for those that survived a bottleneck in Finland and that there are extensive internal structure due to recent population movement, isolation, and drift 21,33-35 .

Time and memory considerations for eGRM algorithm
We implemented the eGRM in a flexible Python framework using custom C extensions to accelerate core eGRM calculations. Our implementation of eGRM is memory-efficient. The main memory usage throughout the algorithm is a matrix of size × , which takes 8 ! bytes of memory when stored as doubles in C. However, outputting the resulting matrix into a NumPy array dominates the overall memory consumption (Figure S7). The time cost of computing eGRM is ( ! ), where is the number of genealogical trees (see Supplementary Text). In the case of the FinMetSeq data, the genome-wide eGRM takes ~30 hours on a single CPU to compute for 2644 samples and ~120,000 genealogical trees.

Discussion
In the current study we introduce the eGRM, a genealogical estimate of genetic relatedness. The eGRM is conceptually distinct from the canonical GRM, which is variant-or mutation-centric. As a result, analyses utilizing the canonical GRM need to be interpreted within the context of the marker ascertainment. Ascertainment could be biased because of availability of data, technical errors in data generation, or inconsistent analytical conventions across analysts. In contrast, the eGRM does not depend directly on the detection of variation (eGRM based on the true ARG does not depend on variation; but haplotypes based on a set of variant is used for inferring the ARG in practice), and thus is more robust when used in analyses with incomplete data.
A number of methods have been proposed to exploit the rich information stored in the ARG to make inference of population genetic parameters (e.g. for selection 36,37 or population history 38,39 ).
Similarly, recent theoretical work has demonstrated the relationship between mutational processes by site, or on branches and nodes of the ARG 40 . Given the ARG, our framework considers mutations as appearing uniformly at random on the ARG, and relatedness between pairs of individuals is based on the probability of shared mutation, which is proportional to the branch lengths relating the two individuals. Our decision to explore this framework over alternative paths such as manipulating a matrix of TMRCAs is driven by the conceptual shift in treating mutations as random. We expect that our genealogical framework to compute genetic relatedness will enable seamless incorporation with downstream statistical applications, such as its inclusion in a linear mixed model for controlling population structure in association testing.
Through extensive simulations, we demonstrated that the eGRM is highly correlated with TMRCA and importantly provides an approximately unbiased estimate of Kall. The former we used as the standard for true relatedness given a single tree, while the latter is an idealized GRM assuming all variants are perfectly observed. We illustrated the improvement and new insights that could be garnered using the eGRM with an application in the detection of population structure. First, common SNPs were thought to be uninformative about recent population structures because they tend to predate recent population divergence and are likely shared across all populations 29 .
However, since haplotypes derived from common variants are used to infer the ARG, we showed that the eGRM detects the recent split and can better separate the spatial structure among demes, despite relying on only the common SNPs. Second, our framework provides a means to flexibly probe into the population structure at arbitrary time depths through the partial eGRM, suggesting that it can be optimized to account for structure at varying time scales. Third, we demonstrated empirically the insights of population structure that can be learned from eGRM  Figure 4A). Finally, by examining the partial eGRMs, the structure in Northern and Eastern Finland appeared to be more recently established, and the pattern of variation from this dataset is more consistent with a recently structured population with enhanced drift, rather than the conventional belief that Finnish composes a single homogeneous population 21,42 .
We have found that eGRM inference is stable when using computationally reconstructed genealogies rather than the true ARG (e.g., EKrelate). However, we note that underlying assumptions required for accurate ARG inference may not be met in massive sample sizes. For example, it is often assumed that there are no recurrent mutations and multiple coalescent events per generation; both of these assumptions would be violated in extremely large samples.
Mutation rate is often assumed to be constant across all branches of the inferred tree, which may not be true empirically 43 , and ARG inferences currently may not appropriately model all recombination events 44 . The extent to which relaxation or violation of these assumptions impact the ARG inference and its downstream computation of the eGRM will need to be evaluated systematically. Furthermore, a major current impediment is the scalability of ARG inference. In our simulations, eGRM from RELATE-reconstructed ARGs performs better than that of TSINFER+TSDATE, perhaps due to the inferred branch length information in RELATE, but the computational time is also orders of magnitude longer. As a result, even though the eGRM computation usually takes less than 5% of the total runtime, we were unable to efficiently compute eGRM on RELATE-reconstructed ARGs beyond 10,000 diploid individuals.
Nevertheless, computational advances may well continue to make ARG inferences more accurate and scalable; until recently ARG inferences were restricted to only tens of individuals.
As faster or more accurate ARG inference algorithms become available, our method will be primed to achieve advanced usability and performance with little adjustment.
Even without a more scalable ARG inference methods, eGRM will have the potential to make an

Pairwise genetic relatedness with unobserved markers
We first describe the canonical and expected GRM (eGRM) in a haploid scenario; our framework can easily be generalized to diploid scenarios as describe in Supplementary Methods. We where is the all-ones vector and ( ) is the number of variants in . However, in practice only * < markers are observed, resulting in the × ′ observed haplotype matrix ′. The GRM computed using ′ is given by ( ′). As its definition suggests, ( ′) reflects the pair-wise relatedness conditioned on the observed haplotype data ′, and provides an incomplete picture of relatedness measure between individuals. Even though the complete haplotype matrix is unknown ( is also unknown), through linkage between the unobserved markers in and observed markers in ′, there may exist a more reasonable estimate of ( ) than ( ′) itself.
Specifically, we show that when the ancestral recombination graph (ARG) that connects the samples is given, the expectation of ( ) can be derived, which we denote as the expected GRM (eGRM) on .

eGRM ∶= ( ( )| )
In practice, the ancestral recombination graph can be inferred from the observed haplotypes ′ using recently emerging genealogical tree reconstruction tools 27,28 . Intuitively, we wish to define an eGRM where each entry represents the expected similarity between a pair of individuals should a mutation arise randomly in the ARG, after accounting for the expected similarity between a random pair of individuals in the same context given the structure of the ARG.

Expectation of pairwise relatedness given a genealogy
The recombination and coalescent history of the haploid samples can be completely represented by the ancestral recombination graph 25,26 (ARG) denoted as , which consists of a sequence of genealogical trees = { + |1 < < } across the whole genome. Each tree + = ( + , + ) is a directed binary tree with each node ∈ + representing a chromosomal segment of a sample or an ancestor, and each branch ∈ + representing the history of its child node until it coalesced into its parent node. We define ( ) as the haplotype vector associated with , that is We assume is fixed, and is generated randomly through mutations occurring on , implying expectation or variance over is conditional on by default. We overload set membership notation over and denote that is a branch in the ARG as ∈ , to mean ∈ + ∈ + ∈ for some .
Here, we define how mutations arise on . For each branch ∈ + we define ( ) as its length in

Genealogy and genotype simulations
We used two different demographic models in our simulation experiments ( Figure 1B)  We simulate genotypes and tree sequences by MSPRIME 50  We phased the genotypes by EAGLE using its default hg19 genetic map. We reconstructed the genealogical tree sequence using RELATE with all parameters same as in its official manual. We then applied our eGRM algorithm on the resulting tree sequence to compute EKrelate. The Canonical GRM was computed based on the same set of 208,681 SNPs after phasing with EAGLE.

Population structure analysis
We contrasted and visualized the information of population structure contained in GRM and eGRM through principal components analysis (PCA) and uniform manifold approximation and projection (UMAP). PCA was computed using the "linalg.eig" function in the python "numpy" library, UMAP was computed by the R "umap" package with all default parameters. To which is a real number between 0 and 1 indicating how well the metric is capturing the true classification . Note that SI is only dependent on the relative order of distances between pairs of points, making it a unified measure of clustering performance among PCA, UMAP and other distance-based methods.

eGRM Software
We have implemented the algorithms related to eGRM in a python package 'egrm', which is publicly available in PyPI. Documentation of this package as well as simulation commands used in this study can be found on its github page (https://github.com/Ephraim-usc/egrm).

Acknowledgement
We  Each mutation has a corresponding length-3 haplotype vector (e.g., mutation m1 has haplotype vector (1, 1, 0)). The "single-variant GRM" can be computed as the outer product of the centered and normalized haplotype vector (e.g., ⁄ 1 2 , ⁄ 1 2 , − 2 for mutation m1) with itself. The canonical GRM is then the unweighted average of the single-variant GRMs of the 6 mutations. The eGRM is based on the 4 branches, weighted by their lengths (i.e., the expected number of mutations on this branch). (B) Overview of simulation workflow to test the performance of eGRM. ARGs are simulated by MSPRIME based on a single-growth demographic model and a grid-like spatial structure model. Observed variants are oversampled from common variants to mimic real genotyping array data. We then compute the complete GRM (K all ), the observed GRM (K obs ), eGRM based on the true ARG (EK), eGRM based on RELATE-or TSINFER+TSDATE-reconstructed ARG (EK relate and EK tsdate ).

Figure 2 eGRM is highly and unbiasedly correlated with measures of relatedness in simulations. (A) Negative
Spearman correlation between TMRCA and K obs , EK (left) or EK relate (right) on a 1Mb non-recombining locus. Spearman correlation is used because GRM by definition normalizes according to allele frequency to upweight rare mutations, and thus is not expected to correlate linearly with TMRCA. (B) Heatmap summarizing the Pearson correlations between GRM and eGRM matrices on a 30Mb chromosome. Note that EK and EK relate are highly correlated with K all . (C) Scatter plots of the GRM and eGRM values for all pairs of individuals, using the same simulated data as in B. All simulations from A to C simulated 1000 individuals. (D) Pearson correlation with K all , with varying proportion of SNPs observed (sample size is fixed to 1000; left) or varying sample size (20% common SNPs observed; right) on a 30Mb chromosome.

Figure 3
PCA based on eGRM more effectively captures recently established spatial structure in simulation compared to the canonical GRM. A 30Mb region was simulated, with 20% of common variants observed. Each deme has a constant population size of 500, in which 50 individuals are sampled. The first two PCs based on PCA of three GRMs are shown (top): the canonical GRM based on observed SNPs (K obs ), the eGRM based on RELATE-reconstructed ARG using the observed SNPs (EK relate ), or the partial eGRM based on the subset of branches between 0 and 100 generations across the ARG (EK relate (0-100 gen)). Separation index (SI; see Method for a precise definition) is shown at the topright corner of each plot, which is the average proportion of same-label neighbors for each sample, indicating how well populations are separated. The first two features of UMAP transformation (bottom) applied to the top 10 PCs further accentuate the detected structure, as measured by SI.

Figure 4
Clustering analysis based on eGRM revealed novel population structure in the population of Northern and Eastern Finland. (A) PCA and PCA+UMAP based on either K obs or EK relate are shown. A map of Finland with regions colored is provided for reference. Main geographical locations referenced in the text are labeled (Surrendered Karelia colored in magenta, Lapland colored in red, Turku-Pori colored in light green, and Vassa colored in light blue). Scatterplot of the first two features of UMAP transformation was based on the first 24 and 58 components of PCA of K obs and EK relate , respectively. These numbers were chosen as they respectively are the number of components at which the separation index (SI) is maximized after applying the UMAP transformation. (B) Separation index achieved as successive PCs were included in UMAP transformation of PCA of K obs and EK relate .