A likelihood ratio approach for identifying three-quarter siblings in genetic databases

The detection of family relationships in genetic databases is of interest in various scientific disciplines such as genetic epidemiology, population and conservation genetics, forensic science, and genealogical research. Nowadays, screening genetic databases for related individuals forms an important aspect of standard quality control procedures. Relatedness research is usually based on an allele sharing analysis of identity by state (IBS) or identity by descent (IBD) alleles. Existing IBS/IBD methods mainly aim to identify first-degree relationships (parent–offspring or full siblings) and second degree (half-siblings, avuncular, or grandparent–grandchild) pairs. Little attention has been paid to the detection of in-between first and second-degree relationships such as three-quarter siblings (3/4S) who share fewer alleles than first-degree relationships but more alleles than second-degree relationships. With the progressively increasing sample sizes used in genetic research, it becomes more likely that such relationships are present in the database under study. In this paper, we extend existing likelihood ratio (LR) methodology to accurately infer the existence of 3/4S, distinguishing them from full siblings and second-degree relatives. We use bootstrap confidence intervals to express uncertainty in the LRs. Our proposal accounts for linkage disequilibrium (LD) by using marker pruning, and we validate our methodology with a pedigree-based simulation study accounting for both LD and recombination. An empirical genome-wide array data set from the GCAT Genomes for Life cohort project is used to illustrate the method.


Introduction
The detection of related individuals in genetic databases is of great interest in various areas of genetic research. Most obviously, it is of interest in forensic studies aiming at identifying relationships between individuals such as paternity tests (Evett and Weir, 1998) (Mo et al., 2016, Wang, 2004. Good high-resolution techniques for detecting related individuals are also of interest for genealogical research on family reconstruction (Staples et al., 2014). In conservation genetics, careful selection of unrelated individuals for breeding programs is needed (Oliehoek et al., 2006), requiring the estimation of pairwise genetic relationships. In genome-wide association studies (GWAS) that have become popular during the past two decades (Visscher et al., 2017), standard quality control filters are applied prior to genetic association analysis. The presence of cryptic relatedness violates the assumption of independent individuals in association modeling. For this reason, removing related individuals in the genetic database prior to the GWAS analysis is a common quality control step (Anderson et al., 2010). Many methods for relatedness research are described in the literature. Most of them are based on the principle of allele sharing. Two individuals can share 0, 1, or 2 alleles for a diploid genetic marker. These alleles can be identical by state (IBS) or identical by descent (IBD). A scatterplot of the mean (x IBS ) and standard deviation (s IBS ) of the number of IBS alleles over variants can be used to identify related pairs (Abecasis et al., 2001). Alternatively, a scatterplot of the proportions of sharing 0, 1, or 2 IBS alleles (p 0 , p 1 , p 2 ) is also often used to detect related pairs (Rosenberg, 2006). In genetic studies, the probabilities of sharing 0, 1, and 2 IBD alleles (k 0 , k 1 , k 2 ) can be estimated and used for relationship inference, since their theoretically expected values are known for the standard relationships (see Table 1). For example, parent-offspring pairs have (k 0 , k 1 , k 2 ) = (0, 1, 0) and full siblings have (k 0 , k 1 , k 2 ) = (0.25, 0.50, 0.25). For a given pair of individuals, these probabilities can be estimated by maximum likelihood (Milligan, 2003, Thompson, 1975, 1991, by the method of moments (Purcell et al., 2007) or with robust estimators (Manichaikul et al., 2010). From these probabilities, the kinship coefficient, defined as ϕ = k 1 /4 + k 2 /2, can be obtained. The kinship coefficient can be used to remove individuals with first degree (parent-offspring (PO) or full siblings (FS)) and second-degree relationships (half-siblings, avuncular or grandparent-grandchild) by retaining only pairs with ϕ < 1/16. In addition, third-degree relationships (first cousins (FC)) can be eliminated by retaining only pairs with ϕ < 1/32 (Anderson et al., 2010). All these methods have in common that the inference of the family relationships is based on the judgment of the analyst of the point estimates ( k 0 ;k 1 ;k 2 ;φ) or of a graphical representation ((x IBS ,s IBS ), (p 0 , p 1 , p 2 ) or (k 0 ;k 1 ;k 2 )) ( Galvan-Femenia et al., 2017).
The sample size used in genetic studies, GWAS in particular, is progressively increasing owing to large human sequencing projects that involve genetic data from hundreds of thousands of individuals such as UK Biobank (Bycroft et al., 2018), gnomAD (Karczewski et al., 2020, TOPMed (Taliun et al., 2019), and DiscovEHR (Staples et al., 2018) among others. With such large databases, it becomes increasingly likely that in-between 1st and 2nd degree, and in-between 2nd and 3rd-degree relationships are found. Such in-between relationships are mostly ignored in a relatedness analysis, which typically mostly focus on 1st, 2nd, and 3rd-degree relationships. In this paper, we therefore develop a likelihood ratio (LR) approach that will allow us to identify three-quarter siblings (3/4S), a family relationship whose individuals share fewer alleles than 1st-degree relationships but more alleles than 2nd-degree relatives (Table 1). A 3/4S pair has one common parent, whereas their unshared parents have a first-degree relationship (FS or PO; see Graffelman et al. 2019 Fig. S2). The IBD probabilities for 3/4S are (k 0 , k 1 , k 2 ) = (3/8, 1/2, 1/8) and their kinship coefficient is ϕ = 3/ 16. A detailed derivation of these probabilities is shown in Appendix A. A 3/4S relationship is not very common, but is more likely to be present in GWAS studies with everincreasing sample sizes. The 3/4S relationship has received very little attention in the literature, and the aim of this paper is to develop tools that distinguish 3/4S from full siblings and second-degree relatives.
The remainder of this paper is structured as follows. Section "Methods and materials" develops a LR approach for identifying three-quarter siblings. Section "Simulations" evaluates the LR approach in a simulation study. Section "Case study" illustrates our approach with genome-wide SNP array data from the GCAT Genomes for Life project cohort. Finally, we end the article with a discussion of the proposed methodology.

Methods and materials
Overview of the likelihood of a relationship A detailed derivation of the likelihood of having a given relationship is given by Wagner et al. (2006). In brief, let n be the number of individuals in a non-inbred homogeneous population and assuming absence of population structure. We consider biallelic genetic variants with alleles A and B having allele frequencies p and q, respectively. Let G 1 /G 2 be the genotypes for a pair of individuals, let k m with m = 0, 1, 2 be their IBD probabilities (shown in Table 1) and let R be their family relationship. Then, the probability of observing G 1 /G 2 , given R is: The terms P(G 1 /G 2 |m = 0), P(G 1 /G 2 |m = 1) and P(G 1 / G 2 |m = 2) are the probabilities of observing each pair of genotypes given the number of IBD alleles ( Table 2).

The LR approach for identifying three-quarter siblings
The LR approach has been widely used for relatedness research during the last decades (Boehnke and Cox, 1997, Heinrich et al., 2016, Katki et al., 2010, Kling and Tillmar, 2019, Thompson, 1986, Weir et al., 2006. In brief, the LR approach is based on the contrast of two hypotheses, one in the numerator, H i ; and the other one in the denominator, H j . The larger the LR, the more plausible is H i ; whereas the smaller the LR, the more plausible is H j . For relatedness research, we consider the ratio of the probabilities from Eq. 1 according to the hypothesis of the R i and R j relationships. That is: Here we consider the FS, 3/4S, 2nd, and unrelated (UN) relationships and calculate three LR having FS, 3/4S, or 2nd in the numerator and having the UN relationship in the denominator. The common denominator makes the LR values comparable in order to distinguish 3/4S from FS and 2nd degree. The inference of relatedness for each pair of individuals is based on the largest LR value in the FS~UN, 3/4S~UN, and 2nd~UN ratios. The LRs are shown in Table 3, depending on the observed genotypes of a pair of individuals. Most of these ratios are derived in Heinrich et al. (2016), and the new results for 3/4S are shown in Appendix B. The e parameter from the PO~UN ratio in Table 3 is a small number (i.e., 0.001) used to account for genotype errors and de novo mutations if the genotype combination does not occur. In this way, the LR cannot be zero. For S biallelic SNPs, the LR can be obtained by multiplying the LR s across independent markers and by dividing by the number of SNPs. It is convenient to work in a logarithmic scale such that: which corresponds to the logarithm of the geometric mean of the LRs. Obtained LRs are subject to uncertainty. To assess this uncertainty, we propose to apply bootstrap resampling (Efron and Tibshirani, 1994). This allows the construction of 95% bootstrap confidence intervals for the LRs, which are helpful to assess which relationship is the most likely one for a given pair.

Materials
We test our method for detecting 3/4S with data from the GCAT Genomes for Life cohort project . In brief, the GCAT project is a prospective study that includes~20K participants recruited from the general population of Catalonia, a Western Mediterranean region in the Northeast of Spain. A subset of 5459 participants was genotyped using the Infinium Expanded Multi-Ethnic Genotyping Array (MEGAEx) (ILLUMINA, San Table 2 Possible pairs of biallelic genotypes and the probability of each pair given the number of alleles identical by descent (m).
We assume that the order of the genotypes is irrelevant, i.e., the probabilities for G 1 /G 2 and G 2 /G 1 are the same.  The considered LR are PO, FS, 3/4S, 2nd, or FC relationships in the numerator and the UN relationship in the denominator. The LR values depend on the observed genotypes of a pair of individuals and the allele frequencies p and q of the population under study. The e parameter is used to account for genotype errors and de novo mutations if the genotype combination does not occur (Heinrich et al., 2016). We assume that the order of the genotypes is irrelevant, i.e., the LR for G 1 /G 2 and G 2 /G 1 is the same.

Simulations
In this section, we evaluate the likelihood ratio approach to distinguish 3/4S from FS and 2nd relationships by using simulated data. Pedigrees were simulated from the genetic data of the individuals of the GCAT project, using the pedsim method of Caballero et al. (2019). We apply this method in order to account for recombination by using sexspecific genetic maps (Bherer et al., 2017) and also a crossover interference model (Campbell et al., 2015). The simulations were carried out as follows. First, we identified 4147 potentially unrelated individuals with kinship coefficient <0.025. From these individuals, we retained 537,488 autosomal SNPs with minor allele frequency (MAF) > 0.01, Hardy-Weinberg exact mid p value > 0.05 (Graffelman and Moreno, 2013) and missing call rate zero. Genotypes of the unrelated individuals were phased with SHAPEIT4 (Delaneau et al., 2019) and were used as input for the ped-sim method. Then, we simulated 500 pedigrees containing one FS pair and 500 pedigrees containing one 3/4S pair. In total, we used 3000 random GCAT individuals as founders to generate 3000 artificial individuals. The number of simulated related pairs were 4,000 PO, 500 FS, 500 3/4S and 3,500 2nd degree from a total of 17,997,000 of pairs. To estimate the IBD probabilities and the kinship coefficient for these simulated pairs we used 27,087 SNPs obtained by retaining variants with MAF > 0.40 and by LD pruning, requiring markers to have low pairwise correlation (r 2 < 0.20). Figure 1 shows the ðk 0 ;k 1 Þ-plot for these simulated pairs of individuals. The IBD probabilities were estimated with the PLINK software (Purcell et al., 2007). As expected, the estimated IBD probabilities are close to the expected theoretical values from Table 1 for most pairs of individuals. In Fig. 1, the 3/4S relationships show good separation from 2nd-degree relationships but mix to some extent with FS pairs. Estimated IBD probabilities appear to be centered on their expected values for FS, 3/4S, and 2nd-degree pairs, and have larger variance then PO and UN pairs. The discriminative power of our method crucially depends on the variance of these estimated probabilities (Hill and Weir, 2011).
Boxplots of the kinship estimator recently proposed by Goudet & Weir (Goudet et al. (2018), Weir and Goudet (2017)) shown in Fig. 2 clearly show a difference in median for 3/4S and 1st-and 2nd-degree relationships, though the distribution of the kinship coefficient of the 3/4S overlaps with those of 1st and 2nd-degree pairs. Also, kinship coefficients can be identical for different relationships, as is the case for PO and FS. Therefore, according to Eq. (3), we calculate the FS~UN, 3/4S~UN, and 2nd~UN likelihood ratios for 500 2nd, 500 3/4S, and 500 FS simulated pairs. Figure 3 shows that FS pairs mostly have the largest LR values in the FS~UN ratio, 3/4S pairs mostly have the largest LR values in the 3/4S~UN ratio and 2nd-degree  pairs mostly have largest LR in the 2nd~UN. Note the plotted data profile shaped in a "greater-than" sign (">") pattern suggesting the inference of 3/4S for most 3/4S pairs. In fact, the correct classification rate of the LR approach for the 2nd, 3/4S and FS simulated pairs is 500/500 = 1, 479/ 500 = 0.958 and 475/500 = 0.95, respectively. When comparing the correct classification rate of the LR approach with the LR-kinbiplot approach (Graffelman et al., 2019) based on 500 FS, 500 3/4S, 3,500 2nd, and 5,000 UN simulated pairs (Fig. S1), we observe slightly lower classification rates for 3/4S (478/500 = 0.956) and FS (468/500 = 0.936) using linear discriminant analysis and slightly better classification rates for 3/4S (481/500 = 0.962) and FS (483/500 = 0.966) when using quadratic discriminant analysis as a predictive model. These simulations show the proposed LR approach to be useful for distinguishing 3/4S relationships from FS and 2nd-degree relationships, and to have similar performance to the previously proposed LR-kinbiplot approach.

Case study
In this section, we apply the likelihood ratio approach to genome-wide SNP array data from the aforementioned GCAT project. Graffelman et al. (2019 , Table 5 and Fig. 7) suggested this database to contain eight 3/4S pairs using a log-ratio biplot approach combined with discriminant analysis (LR-kinbiplot). Figures 4 and 5 show the ðk 0 ;k 1 Þ-plot and boxplots of kinship estimates of the GCAT data. The IBD probabilities were estimated with the PLINK software, whereas the kinship coefficient was estimated by the estimator proposed by Weir and Goudet (2017). The colors for the FS, 3/4S, and 2nd-degree pairs in both Figures were assigned according to the relationship inferred by the LR approach. Figure 4 shows, like the simulations, a larger variance for FS pairs, and smaller variances for PO and UN pairs.    Figure 6 shows the LR ratio values for the three relationships (FS~UN, 3/4S~UN and 2nd~UN ratios) on the horizontal axis, for the presumably FS, 3/4S and 2nd pairs from the GCAT project. The LR analysis reveals eight 3/4S pairs (black color) that have the 'greater-than' sign (">") shaped pattern, because the largest LR values are obtained for the 3/4S~UN ratio. All inferred FS pairs (blue color) have a monotonously increasing shaped pattern ("/") since the largest LR values are obtained for the FS~UN ratio; and all 2nddegree pairs have a monotonously decreasing pattern ("\") since the largest LR values are obtained for the 2nd~UN ratio. Table 4 shows the LR values for each pair which confirm that there are eight 3/4S pairs in concordance with the LR-kinbiplot approach. We used bootrapping to assess the amount of uncertainty in the LRs. The bootstrap distribution of the LR for the eight hypothesized 3/4S pairs is shown in Fig. 7. This plot shows seven pairs having the entire bootstrap distributions for the two relationships completely separated, and these pairs therefore clearly do not correspond to FS pairs. For one pair (20) the 3/4S relationship is most likely, for having on average the largest LR; however, given the overlap of the two distributions, the evidence for a 3/4S relationship is less compelling for this pair.

Discussion
In this paper, we show that the likelihood ratio approach is useful for distinguishing three-quarter siblings from FS and 2nd-degree relationships. Figure 4 shows that in a standard ðk 0 ;k 1 Þ-plot, 3/4S can easily go unnoticed as FS pairs. The LR approach can be of great help to detect such cases. The LR approach developed in this paper confirmed eight 3/4S pairs previously uncovered by a log-ratio biplot (LR-kinbiplot) approach (Graffelman et al., 2019) for genome-wide SNP array data from the GCAT cohort. The assessment of the precise relationship of a pair based on the numerical values of the LRs, or on a plot of the LRs, ignores the uncertainty in these statistics. We found bootstrap procedures to be extremely useful for quantifying this uncertainty, and consider it to be an invaluable tool for the sensible interpretation of the pairwise LR statistics.
The estimated relationships for the GCAT cohort were to some extent confirmed by an analysis of the surnames of the participants, respecting their privacy. In Spain, people have a double surname, usually the first from the father and the second from the mother. This implies that FS and 3/4S pairs share two surnames, whereas 2nd-degree relationships share only one. All identified 3/4S pairs were confirmed to share two surnames, supporting that these pairs are not 2nd degree.
The proposed LR approach multiplies the likelihoods over loci, under the assumption of independence. The existence of LD between variants violates this assumption. In order to approximately meet the requirement of independence, LD pruning of neighboring variants in a window is therefore recommended (Kling and Tillmar, 2019). This pruning can be done in PLINK (Purcell et al., 2007) or with other software (Calus and Vandenplas, 2018). A future improvement of the LR approach could use Markov chain algorithms Wigginton, 2005, Kling et al., 2015) that allow efficient likelihood computations on blocks of tightly linked markers.
The LR approach developed in this paper assumes known allele frequencies and non-inbred individuals. The first assumption seems reasonable given the large sample size used in this study. Inbreeding could be accounted for by the use of nine condensed Jacquard coefficients (Hanghoj et al., 2019, Jacquard, 1974 in the development of the likelihood ratio. Inbreeding could yield other levels of relationship in-between FS, 3/4S, and 2nd degree. The ðk 0 ;k 1 Þ-plot of the GCAT data in Fig. 4 reveal closeness of the 3/4S and FS pairs, and suggests intermediate relationships like seven-eighths siblings (7/8S) might also exist in the data. Indeed, the full range of 2ND, 5/8S, 3/4S, 7/8S, and FS relationships could be present in the data. It is easily shown that 5/8S and 7/8S have a kinship coefficient of 5/32 and 7/32, respectively. Figure 4 also shows evidence of some pairs in-between a 2nd a 3rd-degree relationship. In future work, the likelihood ratio approach presented in this paper could be further refined to identify all these relationships more precisely. In-between relationships, like the 3/4S relationship studied in this paper, essentially stress that relatedness is a continuous rather than a discrete concept.

Appendix B
Here we show the LR of 3/4S~UN for a biallelic SNP whose alleles are A and B. Let p and q be the allele frequencies for A and B of the population under study. For a pair of individuals, we show the LR computation for four genotype pairs: AA/AA, AA/AB, AA/BB and AB/AB. The LR for the remaining genotype pairs (AB/AA, AB/BB, BB/AA, BB/AB, and BB/BB) are equivalent or can be obtained analogously. The IBD probabilities for 3/4S are (k 0 , k 1 , k 2 ) = (3/8, 1/2, 1/8) and for UN pairs are (k 0 , k 1 , k 2 ) = (1, 0, 0). Then, according to Tables 1 and 2 and Eqs (1) and (2), the LR for 3/4S~UN is derived as follows: