Heterozygote

In the eigenanalysis of the Shriver data, we examine no more than two markers as independent regression variables for each marker we analyze, insisting that any marker that enters the regression be within 100,000 bases of the marker being analyzed. This slightly sharpens the results. Varying these parameters made little difference.
For all STRUCTURE runs, we ran with a burn-in of 10,000 iterations with 20,000 follow-on iterations, and no admixture model was used. Computations were carried out on a cluster of Intel Xeon compute nodes, each node having a 3.06-GHz clock.
For our coalescent simulations, we assumed a phylogenetic tree on the populations, and at each simulated marker, ran the coalescent back in time to the root of the tree. At this point we have a set of ancestors A of the sampled chromosomes. We now assume that the marker is biallelic and that the population frequency f of the variant allele in the ancestral population is distributed uniformly on the unit interval. Sample the frequency f and then choose an allele for each ancestor of A, picking the allele for each ancestor with probability f. Now retain the marker if it is polymorphic in our samples. This process is mathematically equivalent to having a very large outgroup population diverging from the sampled populations at the phylogenetic root, with the population panmictic before any population divergence, and ascertaining by finding heterozygotes in the outgroup. If our simulated samples have n individuals, our procedure yields a sample frequency that is approximately uniform on (1,2,…,2n − 1).
For the admixture analysis that created the plot of Figure 8 we had a population C that was admixed with founder populations A and B. For each individual of C, we generated a mixing value x that is Beta-distributed B(3.5,1.5). Now for each marker independently, the individual was assigned to population A with probability x or B with probability 1 − x.

Table 1

Drop-in and drop-out measurements for 12 samples.

Sample	Input(ng)	Mean coverage	Drop-in	Drop-out
1	20	657	0.0195	0
2	20	40	0.0227	0
3	20	20	0.0237	0.0040
4	20	8	0.0243	0.0537
5	20	4	0.0220	0.2597
6	20	2	0.0167	0.6160
7	20	1	0.0067	0.9223
8	1	105	0.0179	0
9	0.25	59	0.0143	0.0196
10	0.125	35	0.0129	0.0822
11	0.031	12	0.0092	0.5558
12	0.015	9	0.0091	0.7461

Drop-in is measured as the proportion of reads indicating an allele that is not present in the genotype the read reports for. Drop-out is measured as the proportion of heterozygous loci where there are reads for only one allele type

The assessment of our observation model in Sect. 3.1 uses the Coriell sample NA12878, a genomic reference material (Coriell Institute) sequenced in Tillmar et al. [13 (link)] on a Illumina MiSeq instrument. Results for samples with varying amounts of DNA and varying allelic depths are shown in Table 1.
The simlation uses 3929 SNPs from Tillmar et al. [13 (link)], describing a SNP panel with autosomal SNPs evenly spread across the chromosomes. Genetic positions are downloaded from Ruther’s repository [22 (link)]. From Tillmar et al, we further use genotype data for the Coriell sample NA12878 to obtain coverage statistics. Allele frequencies are extracted for individuals with European ancestry (CEU) from the 1000 Genomes project [23 ]. We generate founder alleles through the population model in Sect. 2.3 with $\theta =0.01$ and $\gamma =0.001$ . We continue to drop alleles through the pedigree using the inheritance model in Sect. 2.4, with crossover probabilities derived from the genetic positions alluded to above. Next, to mimic low coverage data (lcNGS) based on reduced-quality samples, we use the model in Sect. 2.2 with $m=10$ and $e=0.02$ to generate sequence read data. The allelic depths are drawn independently for each locus using a discretized Gamma distribution, first with expectation 10 and standard deviation 2 for Figs. 2 and 3 and then with expectation 3 and standard deviation 1 for Figs. 4 and 5.Fig. 1

Our five example pedigrees. For each we indicate the numbering of the persons, the numbering of the parent–child relationships and which persons are tested (filled symbols)

In the simulation study in Sect. 3.3 we focus on whether two persons are second cousins (see Fig. 1) or unrelated. For each relationship 1000 cases are simulated and a Likelihood Ratio in favour of relatedness is computed using three different methods: Our proposed method, an amended version where linkage is ignored, and an amended version where genotypes are called. To make it optimally competitive, the calling algorithm uses the same likelihood as in our model, combining it with prior probabilities for genotypes based on allele frequencies and selecting the genotype that maximizes the resulting posterior. In other words, the called genotype is the one that maximizes the product of the population frequency of the genotype and the likelihood of the data given the genotype, where the likelihood is computed as in Sect. 2.2.
For each simulated case we also estimate Jacquard coefficients using NgsRelate [18 (link)]. We use VCF-files as input, with PL-fields derived from the same data likelihoods we use in our proposed method. The Euclidean distances from the estimated point $k=(k_0,k_1,k_2)$ of non-inbred coefficients to corresponding points representing the second cousin relationship or unrelatedness are computed. Comparing the difference in distances to a cutoff value yields a classification of cases into related or unrelated. Varying the cutoff value yields receiver operating characteristic (ROC) curves seen in Figs. 3 and 5. For comparison, the figures also show results for other methods, converted to ROC curves using the LR as cutoff.

SNP genotypes were denoted as 0/0 for homozygous reference alleles, 0/1 for heterozygous alleles, and 1/1 for homozygous alternate alleles (0: reference allele; 1: alternate allele). Association analysis of logistic regression was performed using the Python package statsmodels (Seabold et al., 2010 ). An additive model was used for the association between the SNPs and AMD. For additive logistic regression analysis, homozygous reference alleles, heterozygous alleles, and homozygous alternate alleles were respectively defined as the values 0, 1, and 2. The clinical data mining and management of the SQL server in TCVGH was conducted using Microsoft Azure Data Studio. Patient comorbidities included hypertension (ICD-9-CM codes 401.xx—405.xx), coronary artery disease (410.xx—414.xx), cardiac dysrhythmias (427.xx, 785.0, and 785.1), cerebrovascular diseases (433.xx—438.xx), chronic respiratory diseases (490—496), and hyperlipidemia (272.x). Individuals with any comorbidity were identified through diagnoses performed during at least two ambulatory visits to TCVGH. Statistical significance was defined as a p-value < 0.05.
Survival analysis was assessed by the Kaplan–Meier estimate using the R package survival (Therneau and Grambsch, 2000 ). Observation time was defined as the period of duration from the first outpatient visit for a comorbidity to the first time receiving a diagnosis for AMD. The survival curve was plotted by the R package survminer (https://CRAN.R-project.org/package=survminer). Log-rank tests for significant differences in survival time between the two groups were performed using the survdiff function in the survival package. A Cox proportional hazard (PH) model was used to estimate the hazard ratio (HR) using the coxph function in the survival package. For Cox PH model, homozygous reference alleles (0/0), heterozygous alleles (0/1), and homozygous alternate alleles (1/1) were respectively defined as the values 0, 1, and 2.