A Worked Workflow for PCA and Kinship Estimation

From Genotypes to a Joint Analysis

Principal component analysis (PCA) and genomic relationship estimation both summarize patterns across individuals, but they answer different questions. PCA gives directions of greatest variation in the genotype data and scores for each individual along those directions. A genomic relationship matrix (GRM) summarizes pairwise genotype similarity. A practical workflow can use the same carefully selected markers for both, provided that the marker set and scaling choices are made explicit.

The calculations below use diploid, biallelic autosomal markers coded as allele dosages $g_{ij}\in\{0,1,2\}$, where $i=1,\ldots,n$ indexes individuals and $j=1,\ldots,m$ indexes markers. The allele frequency $p_j$ is estimated in the analysis cohort after applying marker and sample quality checks. The examples are hypothetical and deliberately small; real analyses require many more markers and more extensive diagnostics.

Definition: For a non-missing dosage $g_{ij}$, the frequency-standardized genotype is $$ z_{ij}=\frac{g_{ij}-2p_j}{\sqrt{2p_j(1-p_j)}}. $$ The standardized genotype matrix $Z$ has individuals in rows and markers in columns. A commonly used standardized GRM is $K=ZZ^\mathsf{T}/m$. Its entry $K_{ik}$ is the average standardized genotype cross-product for individuals $i$ and $k$.

This scaling subtracts the dosage mean $2p_j$ and divides by the Hardy–Weinberg dosage standard deviation $\sqrt{2p_j(1-p_j)}$. It is a modeling scale, not a guarantee that a marker’s observed sample variance is exactly one. Monomorphic markers have a zero denominator and must be excluded; very rare alleles can produce unstable scaling and are commonly filtered according to the analysis objective. The usual interpretation also assumes dosage coding is comparable across individuals and that allele frequencies are estimated from an appropriate set of samples.

Prepare the Marker Set Before Computing Components

Quality control and marker selection precede PCA and GRM construction. Check sample and marker missingness, allele frequencies, genotype coding and the consistency of the sample identifiers. The earlier material on missing genotypes explains why different missingness patterns can create spurious structure. For a PCA workflow, one common choice is to replace a missing dosage by its marker mean $2p_j$, which becomes zero after standardization. This prevents missingness itself from contributing a nonzero standardized value, but it does not correct systematic differences in missingness between groups; those still require investigation.

Linkage disequilibrium (LD) pruning reduces redundant information from correlated markers. For two markers, calculate the squared sample correlation $r^2$ between their dosage columns. A pruning procedure typically considers markers within specified genomic windows and removes one of a pair when $r^2$ exceeds a chosen threshold. The window, threshold and rule for choosing which marker to retain should be recorded. Pruning reduces the influence of dense correlated regions; it does not make all retained markers independent.

Worked Example: Standardizing Dosages and Checking LD

Consider four hypothetical individuals, A through D, genotyped at three markers. Their dosages are $M_1=(0,0,2,2)$, $M_2=(0,1,1,2)$, and $M_3=(0,0,2,2)$. At each marker, the sample allele frequency is $p_j=0.5$, so the standardization denominator is $\sqrt{2(0.5)(0.5)}=\sqrt{0.5}\approx0.7071$.

For $M_1$, subtracting $2p_1=1$ gives $(-1,-1,1,1)$. Dividing by $\sqrt{0.5}$ gives $z_1=(-\sqrt{2},-\sqrt{2},\sqrt{2},\sqrt{2})$, approximately $(-1.414,-1.414,1.414,1.414)$. The standardized values for $M_2$ are $(-\sqrt{2},0,0,\sqrt{2})$. Marker $M_3$ duplicates $M_1$.

To check LD between $M_1$ and $M_2$, use the centered dosage columns $(-1,-1,1,1)$ and $(-1,0,0,1)$. Their cross-product is $2$, and their sums of squares are $4$ and $2$. Thus $r=2/\sqrt{4\cdot2}=1/\sqrt{2}$, and $r^2=0.5$. By contrast, $M_1$ and $M_3$ are identical, so their correlation is 1 and $r^2=1$. With an illustrative pruning threshold of $0.8$, retain $M_1$ and $M_2$, and remove one of the duplicate pair $M_1,M_3$. The first pair does not exceed the threshold.

Compute PCA Scores From the Retained Genotypes

After pruning, center and standardize each retained marker as specified, then calculate the eigenvectors and eigenvalues of $Z^\mathsf{T}Z$, or equivalently use a singular value decomposition of $Z$. The eigenvectors describe marker directions; projecting each individual’s row of $Z$ onto an eigenvector gives that individual’s component score. Eigenvector signs are arbitrary: changing the sign reverses all scores on that component without changing the analysis.

Result: If $v_\ell$ is a unit-length eigenvector of $Z^\mathsf{T}Z$ with eigenvalue $\lambda_\ell$, the component scores are $Zv_\ell$. The fraction of the total standardized sum of squares accounted for by component $\ell$ is $\lambda_\ell/\sum_q\lambda_q$, where the sum is over all eigenvalues.

This fraction follows directly from the trace: the sum of eigenvalues of $Z^\mathsf{T}Z$ is its trace, which is the sum of the squared entries of $Z$. Because the eigenvectors are an orthonormal basis, the total sum of squares is partitioned among the component directions. This is a descriptive variance summary, not a measure of how much ancestry or relatedness a component represents.

Worked Example: Finding the First Principal Component

Use the two retained standardized marker columns above. Their squared lengths are $8$ for $M_1$ and $4$ for $M_2$. Their cross-product is $4$, since the centered dosage cross-product $2$ is divided by $0.5$ when both columns are standardized. Therefore

Z^\mathsf{T}Z= \begin{pmatrix} 8&4\\ 4&4 \end{pmatrix}.

The characteristic equation is $(8-\lambda)(4-\lambda)-16=0$, or $\lambda^2-12\lambda+16=0$. Its roots are $\lambda_1=6+2\sqrt{5}\approx10.472$ and $\lambda_2=6-2\sqrt{5}\approx1.528$. Their sum is 12, matching the matrix trace. The first component accounts for $10.472/12\approx0.8727$, or 87.27% of the total standardized sum of squares.

For the leading eigenvalue, the eigenvector equation gives $4v_1+4v_2=\lambda_1v_2$, so $v_2/v_1=(\lambda_1-8)/4\approx0.6180$. A unit-length leading eigenvector is approximately $v_1=(0.85065,0.52573)^\mathsf{T}$. Projecting the four rows of $Z$ onto it gives scores approximately $(-1.9465,-1.2030,1.2030,1.9465)$ for A through D. These scores reflect the computed leading eigenvector, rather than an arbitrary direction in marker space.

Estimate and Inspect the Genomic Relationship Matrix

With the same $Z$ and $m$ retained markers, calculate $K=ZZ^\mathsf{T}/m$. This averages standardized cross-products across markers. In many genomic analyses, off-diagonal entries are treated as estimates on an additive relationship scale, while pedigree kinship coefficients use a different convention: for example, the expected additive relationship of a parent and offspring is twice their kinship coefficient. The matrix $K$ is a genomic estimate, not a pedigree coefficient, and its finite-sample values can be affected by allele-frequency choices, population structure, marker selection and related individuals in the cohort.

Worked Example: Calculating the GRM Entries

For the two retained markers, the standardized rows for A through D are $(-\sqrt{2},-\sqrt{2})$, $(-\sqrt{2},0)$, $(\sqrt{2},0)$, and $(\sqrt{2},\sqrt{2})$. There are $m=2$ markers. For example, the A–B entry is the row cross-product divided by two: $[(-\sqrt{2})(-\sqrt{2})+(-\sqrt{2})(0)]/2=1$. The A–C entry is $-2/2=-1$.

Calculating all entries gives

K=\frac{ZZ^\mathsf{T}}{2}= \begin{pmatrix} 2&1&-1&-2\\ 1&1&-1&-1\\ -1&-1&1&1\\ -2&-1&1&2 \end{pmatrix}.

This tiny, deliberately structured example produces large diagonal values and strongly negative entries. They are the arithmetic result for these two markers, not reliable estimates of pedigree relationships. The cohort is too small, and the markers separate the same four individuals into contrasting genotype patterns. In a real analysis, inspect the distribution of diagonal and off-diagonal entries and whether broad patterns track known sample groups or technical variables before assigning a relationship interpretation.

Visual Diagnostics, Alternatives and Limits

A useful workflow uses several complementary displays rather than relying on a single plot. A scree plot of ordered eigenvalues helps show whether variation is concentrated in a few components. A score plot of PC1 against PC2 can reveal clusters, gradients or isolated samples for follow-up; the plot alone does not establish their cause. A heat map of $K$, with samples ordered by group or by clustering, can make blocks of unusually high or low similarity visible. Also inspect marker and sample missingness, allele-frequency distributions, component loadings and whether the scores align with batch, platform or other recorded variables.

If a component changes markedly when a few markers are removed, check whether those markers are in a high-LD region or have unusual quality. If scores align with missingness or batch, revisit quality control rather than interpreting the pattern as ancestry. If the GRM diagonal is far from the scale expected in the analysis, check allele-frequency estimates, scaling, retained marker count and sample identities. In small cohorts, noisy estimates and sample composition can have a large effect.

PCA on unstandardized centered dosages gives greater weight to markers with larger observed variance; frequency standardization instead gives markers a variance scale based on $2p_j(1-p_j)$. LD pruning is one way to limit redundant contributions, while retaining all markers can be reasonable when a method explicitly accounts for LD or when the aim differs. A GRM based on standardized dosages is also not the only relatedness estimator: other estimators can use different allele-frequency choices or be designed to reduce sensitivity to population structure. These choices are not interchangeable; state the target and assumptions before comparing results.

Validate the inputs.
Confirm sample identifiers, dosage coding, marker quality and missingness patterns; remove markers that cannot be standardized.

Choose markers and scale.
Estimate allele frequencies, standardize dosages, and apply a documented LD pruning rule when appropriate.

Compute and inspect both summaries.
Obtain component scores and eigenvalues, calculate $K=ZZ^\mathsf{T}/m$, and compare plots and matrix patterns with sample and technical information.

A frequent mistake is to treat a PCA score plot as a direct ancestry label or to read every large GRM entry as a close pedigree relationship. The correction is to interpret each as a summary conditional on the chosen samples, markers, allele frequencies and scaling. PCA describes dominant genotype variation; the GRM describes average standardized similarity. Neither alone identifies why a pattern arose, and neither replaces checks for technical artifacts or an explicitly chosen relationship model.

Check Your Understanding

Answer each question, then compare with the solutions below.

For a marker with $p_j=0.25$, what are the standardization denominator and the standardized value for dosage $g_{ij}=2$?
Two dosage columns have centered cross-product 3 and sums of squares 9 and 4. Calculate $r$ and $r^2$. Would the pair exceed an $r^2$ pruning threshold of 0.8?
What is the first-component fraction if the leading eigenvalue is 15 and the eigenvalues sum to 25?
Why can mean-imputing a missing dosage to $2p_j$ fail to resolve a batch-related PCA pattern?

Solutions. (1) The denominator is $\sqrt{2(0.25)(0.75)}=\sqrt{0.375}\approx0.6124$. The centered dosage is $2-0.5=1.5$, so the standardized value is $1.5/0.6124\approx2.449$. (2) $r=3/\sqrt{9\cdot4}=0.5$, so $r^2=0.25$; it does not exceed 0.8. (3) $15/25=0.60$, or 60%. (4) Mean imputation sets a missing standardized value to zero, but does not remove differences in how often or why samples are missing. If missingness differs by batch, that pattern can still affect the analysis through observed genotypes or sample selection and must be investigated directly.