GRAPE requires various reference data to perform preprocessing, quality control, phasing, and imputation. Reference data is also required for the simulation workflow. In order to facilitate collection of reference data, we created a separate workflow that automates this step. It can be run by specifying reference command to the GRAPE pipeline launcher. The workflow downloads data, unpacks it, and performs the required post-processing procedures. If phasing and genotype imputation are required, one should also add the additional flags--phase and --impute to the command. It affects the amount of downloaded data. If these flags are specified, the workflow downloads additional reference dataset to make phasing and imputation possible.

There is another option to download all required reference data as a single file. This file is prepared by us and preloaded on our side in the cloud. It can be done by specifying the additional flag --use-bundle to the workflow. This way is faster, since all the post-processing procedures have already been performed. Reference files consist of three main groups. •

Files for the preprocessing. These files include genetic recombination maps for mapping SNPs coordinates from base pairs (bp) to centimorgans (cM); files with the SNPs information from the 1000 Genomes Project for the SNPs quality control; reference genome of hg37 build; liftOver chain file.

GRAPE have a versatile and configurable preprocessing workflow. One part of the preprocessing is required and must be performed before the relationship inference workflow. It is launched by the preprocess command of the GRAPE. Along with some necessary technical procedures, preprocessing includes the following steps. •

[Required] SNPs quality control by minor allele frequency (MAF) and the missingness rate. We discovered that blocks of rare SNPs with low MAF value in genotype arrays may produce false positive IBD segments. To address this problem, we filter SNPs by minor allele frequency. We remove SNPs with a MAF value less than 0.02. Additionally, we remove multiallelic SNPs, insertions/deletions, and SNPs with the high missingness rate, because such SNPs are inconsistent with IBD detection tools.

There are three relationship inference workflows implemented in GRAPE. These workflows are activated by the find command of the launcher. Workflow selection is made by the --flow parameter. 1.

IBIS + ERSA, --flow ibis. During this workflow IBD segments detection is performed by IBIS, ⁸ and estimation of relationship degree is carried out by means of ERSA algorithm. ¹² This is the fastest workflow.

We added a simulation workflow into GRAPE to perform a precision/recall analysis of the pipeline. It’s accessible by simulate command of the pipeline launcher and incorporates the following steps: (1) pedigree simulation with unrelated founders; here we use the Ped-sim simulation package ²⁴ ; (2) relatedness degrees estimation; (3) comparison between true and estimated degrees. The source dataset for the simulation is taken from CEU (Northern Europeans from Utah) population data of 1000 Genomes Project. ²⁵ As CEU data consists of trios, we picked no more than one member of each trio as a founder. We also ran GRAPE on selected individuals to remove all cryptic relationships up to the 6th degree. Then, we randomly assigned sex to each individual and used sex-specific genetic maps to take into account the differences in recombination rates between men and women. ²⁶ Results of our precision/recall analysis for the simulated datasets are presented in the corresponding Results section.

Distribution of IBD segments among non-related (background) individuals within a population may be quite heterogeneous. ¹³ There may exist genome regions with extremely high rates of overall matching, which are not inherited from the recent common ancestors. ²⁷ Instead, these regions more likely reflect other demographic factors of the population. The implication is that IBD segments detected in such regions are expected to be less useful for estimating recent relationships. Moreover, such regions are potentially prone to false-positive IBD segments.

GRAPE can use two different approaches to address this issue. The first one is based on genome regions exclusion mask, wherein some genome regions are completely excluded from the consideration. This approach was proposed by authors of the ERSA algorithm, see Ref. 13. The mask was computed based on whole-genome sequencing data for European individuals. The computed mask is built-in into ERSA 1.0 algorithm and is used by GRAPE by default.

The second approach is based on the so-called IBD segments weighing. This idea reminds one proposed in the Ancestry DNA Matching White Paper ²⁸ (see description of their Timber algorithm there). The key idea is to down-weight IBD segment, i.e. reduce the IBD segment length, if the segment crosses regions with a high rate of matching. The approach can be briefly described as follows. At first, one should evaluate IBD segments B j for some background population. After that, we break chromosomes into windows W i of fixed length and compute total overlap length c i between found IBD segments B j and each window W i , c i = ∑ j | W ij | , where W ij is an overlap between the segment B j and the window W i . Obtained overlap lengths are transformed into weights w i which are assigned to the windows, w i = f c i ∈ 0 1 . Here f is a weighing function that reflects the following heuristic: if the overlap length c k for a window W k is a relative outlier among all c i , then the value of f c k is close to 0 , otherwise it’s close to 1 . Given that weights for all windows are computed, for each IBD segment G one can compute its weighted length G w using the formula G w = ∑ i w i | G i | , (1) where | G i | denotes the overlap length between IBD segment G and a window W i . Weighted lengths of IBD segments then are used by the ERSA algorithm (while the ERSA mask is disabled).

GRAPE provides an ability to compute the weight mask from the VCF file with presumably unrelated individuals. It breaks each chromosome into 1cM windows to compute overlaps. After that, GRAPE detects outliers among overlap lengths by means of the minimum covariance determinant (MCD) algorithm, ²⁹ and then determines the outliers upper bound h . This upper bound is used to compute weights, w i = h / c i , if c i > h ; and w i = 1 otherwise. Computed mask can be further used as a parameter for relationship inference workflow (see the --weight-mask parameter). As an example, Figure 2 depicts the weight mask computed for the individuals of East Asian Ancestry taken from 1000 Genomes Project. ²⁵ To detect IBD segments, IBIS workflow was used with parameters --ibis-seg-len 5, --ibis-min-snp 400. High matching regions ( w i ≈ 0 ) are marked with the blue color. For comparison, ERSA 1.0 masked regions are depicted with the pick hatching.

GRAPE can be run inside a Docker container. This way is recommended. Another option is to run the pipeline from scratch with all the dependencies pre-installed. We successfully ran the pipeline on Ubuntu 18.04 Linux distribution ³⁰ with eight CPUs and 32 GB of RAM to evaluate the performance on huge datasets with up to 100k samples and millions of SNPs.

GRAPE can utilize multiple cores. For that, one should specify the cores number via the --cores parameter. The default number of cores is equal to the total number of available CPUs minus 1.

The pipeline can be run using Funnel, ³¹ a lightweight task scheduler that implements Task Execution Schema ²⁰ developed by GA4GH. ²² The scheduler can work in various environments, from a regular virtual machines to Kubernetes cluster with the support of resource quotas. We provide several examples of the task specifications for Funnel. Each sample represents a JSON file with the task description. These files are available in the GRAPE GitHub repository within the corresponding funnel subfolder.

To estimate the performance of the GRAPE pipeline we used a machine with eight CPU cores, Intel(R) Xeon(R) CPU E5-2686 v4 @ 2.30 GHz, and 32 GB of RAM. The fastest IBIS + ERSA (--flow ibis) relatedness inference workflow takes about 22 hours to process 100k individuals dataset, see Figure 3. IBIS tool has a quadratic time complexity with respect to the total number of individuals. The addition of KING (--flow ibis-king) increases the total running time by roughly 50%. Performance analysis confirmed that IBIS is a simple and efficient tool. It allows the pipeline to process hundred of thousands of individuals in a reasonable amount of time.

As the first step, reference data must be downloaded. We suppose that reference data to be stored in/media/ref directory. Here and below we specify the --real-run flag. Without this flag GRAPE performs a dry-run.

Listing 1. Reference downloading for the IBIS + ERSA workflow.

As the second step, preprocessing is performed. We suppose that the input file is located at/media/input.vcf.gz, and it’s in hg38 build. Input file location specified by the flag --vcf-file. GRAPE working directory is/media/data. It’s specified by the --directory flag.

Listing 2. Preprocessing for the IBIS + ERSA workflow.

The third step is the relationship inference. It is launched with the the find command. We use IBIS + ERSA workflow that corresponds to the ibis value of the --flow parameter (default).

Listing 3. Relationship inference with the IBIS + ERSA workflow.

GRAPE has an ability to specify additional parameters to the ERSA and IBIS algorithms to control the sensitivity and the false positive rate. These parameters are described below. Default GRAPE parameters are quite conservative. They provide a low false positive rate and low sensitivity in high (9+) degrees. •

[IBIS] --ibis-seg-len. Minimum length of the IBD segment to be found by IBIS. Higher values reduce false positive rate and give less distant matches (default = 7 cM).

The first and second steps are the same as for the previous use case. The third step is launched with specifying of the ibis-king flow parameter. For this case GRAPE performs an estimation of the first three degrees of relationship with KING. Other degrees are estimated with ERSA (see Figure 1). The resulted report of individual relationships for this case contains additional king_degree and kinship columns. KING algorithm has no additional parameters.

Listing 4. Relationship inference with the IBIS + ERSA & KING workflow.

At first, reference data must be downloaded. The GERMLINE tool works with phased data only. So, if input data is unphased, one should download additional reference dataset to perform phasing. For that purpose we use a reference panel from 1000 Genomes Project. This panel takes a considerable amount of disk space (∼25 GB), and requires significant time to download. To download this panel one should specify the --phase flag while using the reference command.

Listing 5. Reference downloading.

The second step is the data preprocessing. We use the --phase flag to apply the phasing procedure during this stage.

Listing 6. Preprocessing for the GERMLINE + ERSA & KING workflow.

The third step is the relationship inference. One should specify --flow germline to use the GERMLINE tool for the IBD segments detection. ERSA parameters are the same as for the IBIS + ERSA workflow. Parameters of the GERMLINE are not configurable. Currently, we use the following sets of GERMLINE parameters: -min_m 2.5, -err_hom 2, -err_het 1. See GERMLINE documentation for the parameters description. ⁷

Listing 7. Relationship inference with GERMLINE + ERSA & KING workflow.

To perform simulation, at the first step one should download full reference dataset by using the reference command with the --phase and --impute flags enabled. The second step is the simulation workflow. For that, one should use the simulate command of the GRAPE launcher.

Listing 8. Evaluation of the IBIS + ERSA workflow on a simulated dataset.

Along with the parameters for preprocessing, ERSA and IBIS tools, simulation workflow has two additional parameters listed below. •

--sim-params-file. File with parameters of simulation for the Ped-sim tool. For more information see. ²⁴ We have prepared several simulation parameters files, and stored them in the GRAPE repository on GitHub.

The output files include a list of kinship matches found in the simulated dataset, precision/recall plots, and a confusion matrix to compare the detected degrees of relationship with the true degrees. Detailed information on the computed metrics is presented in the Results section.

The first step is to download the reference dataset (see the previous section). The second step is to run the simulate command, specifying germline-king flow and the --phase flag.

Listing 9. Evaluation of the GERMLINE + ERSA & KING workflow on a simulated dataset.

GRAPE has a dedicated command to compute the weight mask, compute-weight-mask. It performs IBD segments detection for the input VCF file and then analyse IBD segments distribution to compute the weight mask. The resulting files consist of a weight mask file in JSON format and a visualization of the mask (see Figure 2).

Listing 10. Computation of the IBD segments weighing mask.

For the example above the resulting files are store in the/media/background/weight-mask/directory.

To apply weight mask during the relatedness detection (find command), one should specify the mask file with the --weight-mask parameter. When used, the ERSA 1.0 exclusion mask is disabled.

Listing 11. Usage of the IBD segments weighing mask.

Finally, GRAPE was evaluated on simulated datasets. We performed the simulation using unrelated founders from 1KGP Affymetrix genotype chip data. This chip contains approximately 900k SNPs. The simulation was carried out with the Ped-sim package. Using this tool, we produced several pedigree structures with eight generations, and a maximum degree of relationships equal to 14. Then we joined all of the pedigrees into one dataset, and performed the precision/recall analysis of the GRAPE pipeline for the different flows. For each degree i of relationships we computed precision and recall metrics: Precision i = TP i TP i + FP i ; Recall i = TP i TP i + FN i . (2)

Here TP i , FP i , FN i are the numbers of true positive, false positive, and false negative relationship matches predicted for the degree i . In our analysis we used non-exact (fuzzy) interval metrics. For the 1st degree, we require an exact match. For the 2nd, 3rd, and 4th degrees, we allow a degree interval of ± 1 . For example, for the 2nd true degree we consider a predicted 3rd degree as a true positive match. For the 5th+ degrees, we use the ERSA confidence intervals which are typically 3-4 degrees wide. For 10th+ degrees, these intervals are 6-7 degrees wide. We also plot a confusion matrix for the predicted vs true degrees.

Results of the simulation for the IBIS + ERSA & KING workflow is presented in Figure 4. The following set of parameters was used: --ibis-seg-len 7, --ibis-min-snp 500, --zero-seg-count 0.5, --zero-seg-len 5, --alpha 0.01. Confusion matrix is presented on the right panel in Figure 4. There -1 stands for no-relationship. The pipeline shows recall above 90%+ for degrees from 1 to 5. It detects all relatives with 1-4 degrees. GRAPE found no false positive matches, i.e. it does not find any relationships among unrelated individuals, Recall − 1 = 1 . Precision is above 90% among all detected degrees. We set these parameters as default for this GRAPE workflow. These parameters are quite conservative, i.e. they provide high precision but low sensitivity.

One can relax GRAPE parameters to get more relatedness matches for high degrees. On the other hand, the number of false positive matches increases as well. Figure 5 shows the simulation results for the IBIS + ERSA & KING workflow for slightly relaxed parameters: --ibis-seg-len 5, --ibis-min-snp 400, --zero-seg-count 0.1, --zero-seg-len 5, --alpha 0.01. The number of detected relationships increases significantly for higher ( ≥ 8) degrees. But false positive matches arise as well, i.e. Recall − 1 ≠ 1 .

Simulation results for the GERMLINE + ERSA & KING workflow is presented in Figure 6. We used the same ERSA parameters as for Figure 4: --zero-seg-count 0.5, --zero-seg-len 5, --alpha 0.01. In comparison to Figure 4, one can see that GERMLINE slightly decreases recall for 6th and 7th degrees, but improves it for higher ( ≥ 8) degrees. Precision is above 95% among all detected degrees. No false positive matches were found, i.e. Recall − 1 = 1 . Our experiments showed that, in comparison to IBIS, GERMLINE is better suited for a careful analysis of relatively small cohorts with phased data. The IBIS algorithm is less sensitive, but for segments with length above 7 cM, it produces the same results as GERMLINE, while working with unphased data.

Our simulation experiments showed that both weighing and exclusion approaches reduce the false-positive rate and significantly improves the overall performance of the pipeline. By the way, weighing mask can be better adapted to specific ancestries, and after additional parameters tuning may slightly outperform ERSA 1.0 approach. In Figure 7 comparison between between the ERSA 1.0 exclusion mask and the GRAPE weight mask is presented for a simulated dataset with the founders of East Asian Ancestry from 1KGP. IBIS workflow was used. Weight mask was taken from Figure 2. Parameter --zero-seg-count was varied while using weight mask to achieve the same level of precision. One can see, that with approximately the same precision, the weighing approach gives several percentages higher recall.

In the end, we compared GRAPE with TRIBES. ²³ TRIBES is an earlier open-source pipeline for relatedness detection. The pipeline combines the GERMLINE algorithm for IBD segments detection and the calculation of the genome proportion with zero alleles inferred IBD (IBD0) for each pair to detect the relatedness. If the data is not phased, TRIBES provides an ability to phase data with the EAGLE tool. This part of TRIBES is similar to one of the corresponding GRAPE workflows of IBD segments detection. In contrast to GRAPE, TRIBES estimates degrees of relationship according to expected IBD0 segments proportion ranges. GRAPE uses the ERSA algorithm, which, to our knowledge, is a more advanced approach.

We have run the TRIBES pipeline on the same simulated datasets. Since the simulated datasets contain unphased data, we also applied built-in phasing procedure from TRIBES. The result of the analysis is presented in Figure 8. TRIBES has demonstrated high detection power for distant relationships of up to the 12th degree. Given that many 13+ degree relatives do not share any IBD segments, this is near a theoretical limit. However, TRIBES produces a huge number of false positive matches, see confusion matrix in the right panel of the Figure 8. Since TRIBES lacks options that allow users to control false positives rates by varying pipeline parameters, it becomes a crucial drawback. This obstacle does not allow the TRIBES pipeline to be adapted for applications, where desired precision/recall rate may vary depending on different business or research objectives.

Publicly available datasets were used to test GRAPE. These datasets are available from: •

1000 Genomes Project ²⁵ ;

GRAPE pipeline can be accessed from the following pubic available resources: •

GitHub: https://github.com/genxnetwork/grape.

Archived source code at time of publication: https://doi.org/10.5281/zenodo.6482561. ³⁷

License: GPLv3.