In general, the standard practice for correcting for population stratification in genetic studies is to use principal components analysis (PCA) to categorize samples along different

*ethnic axes*. Price et al. published on this in 2006, and since then PCA plots are a common component of many published GWAS studies. One key advantage to using PCA for ethnicity is that each sample is given coordinates in a multidimensional space corresponding to the varying components of their ethnic ancestry. Using either full GWAS data or a set of ancestral informative markers (AIMs), PCA can be easily conducted using available software packages like EIGENSOFT or GCTA. HapMap samples are sometimes included in the PCA analysis to provide a frame of reference for the ethnic groups.
Once computed, each sample will have values that correspond to a position in the new coordinate system that effectively clusters samples together by ethnic similarity. The results of this analysis are usually plotted/visualized to identify ethnic outliers or to simply examine the structure of the data. A common problem however is that it may take more than the first two principal components to identify groups.

To illustrate, I will plot some PCs generated based on 125 AIMs markers for a recent study of ours. I generated these using GCTA software and loaded the top 5 PCs into R using the read.table() function. I loaded the top 5, but for continental ancestry, I've found that the top 3 are usually enough to separate groups. The values look something like this:

new_ruid pc1 pc2 pc3 pc4 pc5

1 11596 4.10996e-03 -0.002883830 0.003100840 -0.00638232 0.00709780

2 5415 3.22958e-03 -0.000299851 -0.005358910 0.00660643 0.00430520

3 11597 -4.35116e-03 0.013282400 0.006398130 0.01721600 -0.02275470

4 5416 4.01592e-03 0.001408180 0.005077310 0.00159497 0.00394816

5 3111 3.04779e-03 -0.002079510 -0.000127967 -0.00420436 0.01257460

6 11598 6.15318e-06 -0.000279919 0.001060880 0.00606267 0.00954331

I loaded this into a dataframe called pca, so I can plot the first two PCs using this command:

plot(pca$pc1, pca$pc2)

We might also want to look at the next two PCs:

plot(pca$pc2, pca$pc3)

Its probably best to look at all of them together:

pairs(pca[2:4])

So this is where my mind plays tricks on me. I can't make much sense out of these plots -- there should be four ethnic groups represented, but its hard to see who goes where. To look at all of these dimensions simultaneously, we need a 3D plot. Now 3D plots (especially 3D

*scatterplots)*aren't highly regarded -- in fact I hear that some poor soul at the University of Washington gets laughed at for showing his 3D plots -- but in this case I found them quite useful.

Using a library called rgl, I generated a 3D scatterplot like so:

plot3d(pca[2:4])

Now, using the mouse I could rotate and play with the cloud of data points, and it became more clear how the ethnic groups sorted out. Just to double check my intuition, I ran a model-based clustering algorithm (mclust) on the data. Different parameters obviously produce different cluster patterns, but I found that using an "ellipsoidal model with equal variances" and a cluster size of 4 identified the groups I thought should be there based on the overlay with the HapMap samples.

fit <- Mclust(pca[2:4], G=4, modelNames = "EEV")

plot3d(pca[2:4], col = fit$classification)

plot(pca[2:3], cex = fit$uncertainty*10)I had to scale the uncertainty factor by 10 to make the questionable points more visible in this plot, shown as the hollow circles. We will likely drop these samples from any stratified analyses. We can export the cluster assignment by accessing the fit$classification column, and we have our samples assigned to an ethnic group.

Great! This just gave me new ideas for the visualization of my data. But how did you get these 3d graphs moving?

ReplyDeleteWill wrote this post - I'm guessing he wrote a script to generate PNGs over many different angles, then used ImageMagick/montage to combine them all together into an animated gif. That's how I'd do it at least. Will?

ReplyDeleteThe rgl package is rather extensive. I used a function called movie3d (with another function called rotate3d) to generate a series of images. If have the imagemagick tool Stephen mentioned installed, it will automagickly generate an animation for you. Otherwise you can use the "convert" command from imagemagick to create the animation from a stack of images.

ReplyDeleteCould you post the code used to generate the animation?

ReplyDeleteThanks

Why don't you colour your data points in the scatterplot matrices? Instead of using pairs, give the splom function a go (it is part of the lattice package); than you should also be able to colour your labels. I don't really see that the 3D makes too much of a difference until you coloured your points. Also, keep in mind that for papers you can't have a nice rotating image so you're still stuck trying to show something that is 3 dimensional in 2D when you publish (so you'll always occlude something); the splom doesn't occlude anything since it is already in 2D.

ReplyDeleteHi,

ReplyDeleteI was wondering if somebody could help me. I'm just new to R and bioinformatics and I'm still finding it hard to work everything out. I'm currently doing a GWAS study and having done an MDS-plot I noticed there appears to be a large proportion of cryptic relatedness/population stratification. I wanted to investigate this further and so I have tried to follow the excellent example provided in this blog. I have conducted a PCA analysis using the GCTA software and my first question is, how do I load my data into the dataframe called pca, and read in just the the top 5 PCs. Secondly I would really appreciate if somebody could please provide a step by step guide on how the 3D plot was generated and how I would go about coloring my case and control populations. I apologize if these seem rather mundane questions but being just new to this field it can be very confusing.. I would greatly appreciate any help that can be offered.

Thanks,

Erik

Erik,

ReplyDeleteI would start by breezing through Peter Dalgaard's book Introductory Statistics with R. It will give you some pretty good help with the basics of R.

excellent, post, thanks a lot!

ReplyDeleteI deal myself with PCA and PLS (chemo-, rather than bio-informatics) and will try to contribute something back to the community soon...