# Principal components as covariates in a linear model

I'm working with some genetics data, performing linear regressions, and have been advised to control for population structure by performing principal components analysis. My model at the moment is of the structure:

Minor Allele Frequency ~ Age + Age^2


I've looked on many different sites and read many papers but I'm still not understanding how to carry this out, could someone explain to me? I understand I must add the top P.C.'s to the model as covariates but I'm not sure how exactly.

I am looking at the genetic effects of aging and have data on 20,000 individuals from 40-105 years old. My main goal is to see if there are any SNPs that significantly increase/decrease in MAF over age. From some explanatory 'MAF vs. Age' plots there are 'U' shape curves, which is why I added age^2 to my model. My QQ plot of p-values suggested there is still population structure (I already QC'd the data) and so I should add PCs into the models to correct this problem.

Assuming "Frequency" in your post is the phenotype of interest, you would just change your model to:

$$Phenotype \sim Age + Age^2 + PC_1 + PC_2 + ... PC_k + genotype$$

with $k$ top principal components estimated from the genotype matrix. Standard empirical choices could be between 3-10 top PCs.

If $age$ itself is the phenotype, one can apply a similar model in a reverse regression style: $$genotype \sim Age + Age^2 + PC_1 + PC_2 + ... PC_k$$

but more appropriate would be explicit modeling of age as a time-to-event outcome, i.e.: (using R-style notation) $$Surv(Age, death) \sim PC_1 + PC_2 + ... PC_k + genotype$$ Survival analysis on genome scale is available in ProbABEL, but the model choice there is limited. You might want to ask again on Bioinformatics SE if you're looking for some other tools.

The background of using PCA here is that human populations have different genotype distributions, but also differ in phenotype because of non-genetic reasons (e.g. specific environmental exposures). Thus, genotype-phenotype associations are confounded by population structure, and conditioning on PCs as a proxy for the structure should control for this. Granted, it is a bit crude, and state of the art now is to fit a mixed linear model with a covariance matrix estimated from the genotypes - although the PC control seems to work fine at the level of precision needed in most such studies.

Some references:
Price AL, Patterson NJ, Plenge RM, Weinblatt ME, Shadick NA, et al. (2006) Principal components analysis corrects for stratification in genome-wide association studies. Nat Genet 38: 904–909. https://www.nature.com/articles/ng1847 - the original paper, unfortunately behind a paywall

Zhang F, Wang Y, Deng HW. (2008) Comparison of Population-Based Association Study Methods Correcting for Population Stratification. PLoS One. 2008;3(10):e3392. http://journals.plos.org/plosone/article?id=10.1371/journal.pone.0003392

Price AL, Zaitlen NA, Reich D, Patterson N. (2010) New approaches to population stratification in genome-wide association studies. Nat Rev Genet 11:459–463 https://www.nature.com/articles/nrg2813

• It is not clear to me how PCA would control for dependency. Do you have a reference for that? – Frans Rodenburg Aug 4 '18 at 2:22
• @FransRodenburg intuitively, each PC is assumed to reflect a latent variable along which the populations differ. In an experimentally homogeneous sample, biggest determinant of human genetic variation is ancestry, so top PCs of a genotype sample usually align with geographic axes - e.g. Nordic vs. Southern European. Assuming that the same axis is causal for phenotype variation (let's say through sun exposure, in this example), this looks like a valid approach for me. Edited post to add references. – juod Aug 4 '18 at 4:02
• I see, I'm not familiar enough with GWAS to assess the validity of such an assumption, but thank you for clarifying (+1). – Frans Rodenburg Aug 4 '18 at 4:09
• Thanks juod, that makes things a lot clearer for me. Sorry I should have clarified the 'frequency' is Minor Allele Frequency (MAF) for SNP i (i=1 to 100,000) - would this affect the model you stated? I'm also thinking of creating age bins but am afraid there will be too many covariates in the model? – Erika_Hammerl246 Aug 7 '18 at 18:20
• @Erika_Hammerl246 then I think you should update the question with more info. I really don't understand what's your phenotype / outcome / overall goal - why would you be regressing SNPs against age? Or is it some study of overall mortality? – juod Aug 7 '18 at 20:10