Final Report: Analysis of 1000 Genomes Data with PCA
Madelyn Hoffman
2022-12-16
Introduction
This report summarizes the analysis workflow and results of an analysis of SNPs from the 1000 Genomes Project.
Data preparation
Obtaining and loading data
Single Nucleotide Polymorphism (SNPs) data in VCF format were obtained from the 1000 Genomes Project.
SNPs were downloaded using the Ensembl Data Slicer from chromosome 4 between genomic coordinates 13111147 and 13351147. This represents 0.126% of the chromosome. A total of 7573 variants genotyped in 2504 individuals were downloaded.
The VCF file was loaded into R using the vcfR package
(function read.vcfR) and converted to counts of the minor
allele using the function vcfR::extract.gt().
Meta-data and sample information
Population meta-data was mostly obtained from ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/integrated_call_samples_v3.20130502.ALL.panel. Population codes are defined here: ftp.1000genomes.ebi.ac.uk/vol1/ftp/README_populations.md. This data was merged with SNP data so that rows now contain sample ID, population, super population, sex, latitude, longitude, and genotypes for each person.
Data cleaning
These SNPs were then screened for any SNPs that were
invariant (fixed), resulting in removal of
1912 SNPs (features). This was done using the
invar_omit() function by Nathan Brouwer.
NOTE: The original workflow code for removing invariant SNPs
contained and error that resulted in no columns actually being removed
(Brouwer, personal communication). The code was updated and a
reduction in the size of the dataframe after omitting invariant columns
confirmed by checking the dimensions of the dataframes before and after
this process using dim().
The data were then screened for rows (people) with >50% NAs. There were no NAs in the data, so no rows were removed due to the presence of excessive NAs. Similarly, because no NAs were present no imputation was required.
The data were then centered and scaled using R’s
scale() function. (Alternatively a SNP-specific centering
technique common in other studies could have been applied).
The data were then saved in .csv format using
write.csv() for PCA analysis.
After final processing the data contained 5661 SNPs and 2504 samples (people).
Data Analysis
The code below carries out a PCA on the data and presents the results. The key steps are:
- Load the data with
read.csv(). - Process the data with
prcomp(). - Extract PCA scores.
- Carry out PCA diagnostics, including construction of a scree plot.
- Plot PCs 1 through 3 as scatterplots as pairwise scatterplot.
- Plots PCS 1 through 3 as a 3D scatterplot.
Packages
The following packages were used in this analysis:
# plotting:
library(ggplot2)
library(ggpubr)
# scores() function
library(vegan)## Loading required package: permute## Loading required package: lattice## This is vegan 2.6-4# 3D scatter plot
library(scatterplot3d)Loading data
Load the fully processed data:
NOTE: With 5661 SNPs, this CSV is 257 megabytes. There are more specialized packages for doing PCA with datasets this big. I do not recommend working with more than 10,000 SNPs with basic R functions as we have done in class.
vcf_scaled <- read.csv(file = "prepared-data.csv")Check the dimensions of the data to confirm this is the correct data:
dim(vcf_scaled)## [1] 2504 5667Principal Components Analysis
The data are scaled and ready for analysis. The first six columns contains character data and need to be omitted.
head(vcf_scaled[,c(1:10)])## sample pop super_pop sex lat lng X1 X2
## 1 HG00096 GBR EUR male 52.48624 -1.890401 -0.01998402 -0.01998402
## 2 HG00097 GBR EUR female 52.48624 -1.890401 -0.01998402 -0.01998402
## 3 HG00099 GBR EUR female 52.48624 -1.890401 -0.01998402 -0.01998402
## 4 HG00100 GBR EUR female 52.48624 -1.890401 -0.01998402 -0.01998402
## 5 HG00101 GBR EUR male 52.48624 -1.890401 -0.01998402 -0.01998402
## 6 HG00102 GBR EUR female 52.48624 -1.890401 -0.01998402 -0.01998402
## X3 X4
## 1 -0.06330894 -0.04472137
## 2 -0.06330894 -0.04472137
## 3 -0.06330894 -0.04472137
## 4 -0.06330894 -0.04472137
## 5 -0.06330894 -0.04472137
## 6 -0.06330894 -0.04472137PCA
Principal Components Analysis was run using
prcomp().
vcf_pca <- prcomp(vcf_scaled[,-c(1:6)])Get the PCA scores, which will be plotted.:
vcf_pca_scores <- vegan::scores(vcf_pca) Combine the scores with the sample information into a data frame.
# call data.frame()
vcf_pca_scores2 <- data.frame(population = vcf_scaled$super_pop,
vcf_pca_scores)
# set as a factor
vcf_pca_scores2$population <- factor(vcf_pca_scores2$population)PCA diagnostics
The following steps help us understand the PCA output and determine how many PCs should be plotted and/or used in further analyses such as scans for natural selection, cluster analysis, and GWAS.
Default scree plot
A default R scree plot was created with screeplot().
This plot does not provide extra information for assessing the
importance of the PCs.
screeplot(vcf_pca,
xlab = "Principal Components")
Advanced scree plot
The original workflow and function for making a more advanced scree plot lacked flexibility (Brouwer, personal communication).The following function and workflow simplifies things.
- Run PC (done above)
- Call function
PCA_variation()(below) on PCA output. - (Do NOT call summary() on PCA output)
- Call function
screeplot_snps()on the output ofPCA_variation()to make an advanced scree plot - Call function
PCA_cumulative_var_plot()to show the cumulative variation explained as more PCs are considered
NEW Functions
PCA_variation() function
This function extracts information needed to make a more advanced, annotated scree plot.
# This is a NEW function
PCA_variation <- function(pca){
# get summary information from PCA
pca_summary <- summary(pca)
# extract information from summary
## raw variance for each PC
variance <- pca_summary$importance[1,]
## % variance explained by each PC
var_explained <- pca_summary$importance[2,]*100
var_explained <- round(var_explained,3)
## cumulative % variance
var_cumulative <- pca_summary$importance[3,]*100
var_cumulative <- round(var_cumulative,3)
# prepare output
N.PCs <- length(var_explained)
var_df <- data.frame(PC = 1:N.PCs,
var_raw = variance,
var_percent = var_explained,
cumulative_percent = var_cumulative)
# return output
return(var_df)
}screeplot_snps() function
This functions makes a more advanced scree plot better suited for PCS on for SNPs.
# This is a NEW function
screeplot_snps <- function(var_df){
total_var <- sum(var_df$var_raw)
N <- length(var_df$var_raw)
var_cutoff <- total_var/N
var_cut_percent <- var_cutoff/total_var*100
var_cut_percent_rnd <- round(var_cut_percent,2)
i_above_cut <- which(var_df$var_percent > var_cut_percent)
i_cut <- max(i_above_cut)
ti <- paste0("Cutoff = ",
var_cut_percent_rnd,
"%\n","Useful PCs = ",i_cut)
plot(var_df$var_percent,
main =ti, type = "l",
xlab = "PC",
ylab = "Percent variation",
col = 0)
segments(x0 = var_df$PC,
x1 = var_df$PC,
y0 = 0,
y1 = var_df$var_percent,
col = 1)
segments(x0 = 0,
x1 = N,
y0 = var_cut_percent,
y1 = var_cut_percent,
col = 2)
}PCA_cumulative_var_plot() function
This makes a plot complementary to a scree plot. A scree plot plots the amount of variation explained by each PC. This plot plots a curve of cumulative amount of variation explained by the PCs.
# This is a NEW function
PCA_cumulative_var_plot <- function(var_df){
plot(cumulative_percent ~ PC,
data = var_out,
main = "Cumulative percent variation\n explained by PCs",
xlab = "PC",
ylab = "Cumulative %",
type = "l")
total_var <- sum(var_df$var_raw)
N <- length(var_df$var_raw)
var_cutoff <- total_var/N
var_cut_percent <- var_cutoff/total_var*100
var_cut_percent_rnd <- round(var_cut_percent,2)
i_above_cut <- which(var_df$var_percent > var_cut_percent)
i_cut <- max(i_above_cut)
percent_cut_i <- which(var_out$PC == i_cut )
percent_cut <- var_out$cumulative_percent[percent_cut_i]
segments(x0 = i_cut,
x1 = i_cut,
y0 = 0,
y1 = 100,
col = 2)
segments(x0 = -10,
x1 = N,
y0 = percent_cut,
y1 = percent_cut,
col = 2)
}Advanced screeplot analysis
Extract information
Extract information on the variance explained by each PC.
var_out <- PCA_variation(vcf_pca)Look at the output of PCA_variation()
head(var_out)## PC var_raw var_percent cumulative_percent
## PC1 1 13.665919 3.299 3.299
## PC2 2 10.894489 2.097 5.396
## PC3 3 10.010252 1.770 7.166
## PC4 4 8.898001 1.399 8.564
## PC5 5 8.626755 1.315 9.879
## PC6 6 8.581747 1.301 11.180Advanced screeplot
This advanced scree plot shows the amount of variation explained by all PCs. It marks with a horizontal line what the cutoff is for the amount of Percent variation explained that is useful. The title indicates the percentage value of the cutoff and which PC is the last PC below that value. Though only the first few PCs can be plotted, PCs below the cut off value (“useful PCs) should probably used for further machine learning algorithms.
Make the scree plot with screeplot_snps()
screeplot_snps(var_out)
Cumulative variation plot
The cumulative variation plot shows how much variation in the data explained in total as more and more PCs are considered.The vertical red line shows the cutoff value from the scree plot (above). The horizontal line indicates what the total percentage of variation explained by these useful PCs is.
Make cumulative variation plot with
PCA_cumulative_var_plot()
PCA_cumulative_var_plot(var_out)
PCA Scatterplots
The object created above var_out indicates how much
variation is explained by each of the Principal components. This
information is often added to the axes of scatterplots of PCA
output.
head(var_out)## PC var_raw var_percent cumulative_percent
## PC1 1 13.665919 3.299 3.299
## PC2 2 10.894489 2.097 5.396
## PC3 3 10.010252 1.770 7.166
## PC4 4 8.898001 1.399 8.564
## PC5 5 8.626755 1.315 9.879
## PC6 6 8.581747 1.301 11.180PC1 explains 3.299% percent of the variation, PC2 explains 2.097%, and PC3 explains 1.770%. In total, the first 3 PCs explain ~7.166% of the variability in the data. The scree plot indicate that the first ~647 PCs are useful explain ~76.700% of the variation in the data. In further analysis such as GWAS the first 647 PCs should therefore be used.
Plot PC1 versus PC2
Plot the scores, with super-population color-coded
ggpubr::ggscatter(data = vcf_pca_scores2,
y = "PC2",
x = "PC1",
color = "population",
shape = "population",
main = "PCA Scatterplot",
ylab = "PC2 (2.097% of variation)",
xlab = "PC1 (3.299% of variation)")
Note how in the plot the amount of variation explained by each PC is shown in the axis labels.
Plot PC2 versus PC3
Plot the scores, with super population color-coded
ggpubr::ggscatter(data = vcf_pca_scores2,
y = "PC3",
x = "PC2",
color = "population",
shape = "population",
main = "PCA Scatterplot",
ylab = "PC3 (1.770% of variation)",
xlab = "PC2 (2.097% of variation)")
Note how in the plot the amount of variation explained by each PC is shown in the axis labels.
Plot PC1 versus PC3
Plot the scores, with super population color-coded
ggpubr::ggscatter(data = vcf_pca_scores2,
y = "PC3",
x = "PC1",
ellipse = T,
color = "population",
shape = "population",
main = "PCA Scatterplot",
ylab = "PC3 (1.770% of variation)",
xlab = "PC1 (3.299% of variation)")
Note how in the plot the amount of variation explained by each PC is shown in the axis labels.
3D scatterplot
The first 3 principal components can be presented as a 3D scatterplot.
colors_use <- as.numeric(vcf_pca_scores2$population)
scatterplot3d(x = vcf_pca_scores2$PC1,
y = vcf_pca_scores2$PC2,
z = vcf_pca_scores2$PC3,
color = colors_use,
xlab = "PC1 (3.299%)",
ylab = "PC2 (2.097%)",
zlab = "PC3 (1.770%)")
Based on the PC2 vs. PC1 scatterplot, which shows the most variation in the data, there is little to no apparent grouping of super populations. There is significantly more spread in the African super population. There are a lot of samples all grouped in the same vicinity, which include several from the African super population. The PC3 vs. PC2 and PC3 v. PC1 scatterplots yield similar results. While the African super population has more variety, the SNPs studied from genomic coordinates 13111147 and 13351147 on chromosome 4 do not show inherent clustering of any of the five super populations that distinguishes them from the others.