Final Project
Austin Eastmure
2022-12-10
Prepare R session
Load necessary R packages
library(vcfR)##
## ***** *** vcfR *** *****
## This is vcfR 1.13.0
## browseVignettes('vcfR') # Documentation
## citation('vcfR') # Citation
## ***** ***** ***** *****library(vegan)## Loading required package: permute## Loading required package: lattice## This is vegan 2.5-7library(ggplot2)
library(ggpubr)Confirm your working directory and location of files
getwd()## [1] "/Users/austineastmure/Desktop/comp_bio/Final Project"list.files(pattern = "vcf")## [1] "15.31093190-31333190.ALL.chr15_GRCh38.genotypes.20170504.vcf.gz"
## [2] "all_loci-1.vcf"
## [3] "all_loci.vcf"
## [4] "vcf_noNA.csv"
## [5] "vcf_num_df.csv"
## [6] "vcf_num_df2.csv"
## [7] "vcf_num.csv"Set SNP data up for R
Load the vcf data
my_vcf <- "15.31093190-31333190.ALL.chr15_GRCh38.genotypes.20170504.vcf.gz"
vcf <- vcfR::read.vcfR(my_vcf,
convertNA = T)## Scanning file to determine attributes.
## File attributes:
## meta lines: 130
## header_line: 131
## variant count: 8146
## column count: 2513
##
Meta line 130 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
## Character matrix gt rows: 8146
## Character matrix gt cols: 2513
## skip: 0
## nrows: 8146
## row_num: 0
##
Processed variant 1000
Processed variant 2000
Processed variant 3000
Processed variant 4000
Processed variant 5000
Processed variant 6000
Processed variant 7000
Processed variant 8000
Processed variant: 8146
## All variants processedConvert raw VCF file to genotype scores
vcf_num <- vcfR::extract.gt(vcf,
element = "GT",
IDtoRowNames = F,
as.numeric = T,
convertNA = T)
write.csv(vcf_num, file = "vcf_num.csv", row.names = F)
#Confirm presence of files
list.files()## [1] "07-mean_imputation.html"
## [2] "07-mean_imputation.Rmd"
## [3] "08-PCA_worked.html"
## [4] "08-PCA_worked.Rmd"
## [5] "09-PCA_worked_example-SNPs-part1.html"
## [6] "09-PCA_worked_example-SNPs-part1.Rmd"
## [7] "10-PCA_worked_example-SNPs-part2.html"
## [8] "10-PCA_worked_example-SNPs-part2.Rmd"
## [9] "1000_genomes_snps"
## [10] "1000genomes_people_info2-1.csv"
## [11] "15.31093190-31333190.ALL.chr15_GRCh38.genotypes.20170504.vcf.gz"
## [12] "all_loci-1.vcf"
## [13] "all_loci.vcf"
## [14] "bird_snps_remove_NAs.html"
## [15] "Final Project.Rproj"
## [16] "Final_Project_AE.html"
## [17] "Final_Project_AE.Rmd"
## [18] "Final_Project_Script.R"
## [19] "final_report_AE.html"
## [20] "final_report_AE.Rmd"
## [21] "final_report_template.Rmd"
## [22] "removing_fixed_alleles.Rmd"
## [23] "rsconnect"
## [24] "Script.R"
## [25] "SNPs_cleaned.csv"
## [26] "Software_Checkpoint_Loading_VCF_file_into_R.R"
## [27] "transpose_VCF_data.html"
## [28] "transpose_VCF_data.Rmd"
## [29] "vcf_noNA.csv"
## [30] "vcf_num_df.csv"
## [31] "vcf_num_df2.csv"
## [32] "vcf_num.csv"
## [33] "walsh2017morphology.csv"Transpose original VCF orientation to R dataframe orientation
vcf_num_t <- t(vcf_num)Make into dataframe
vcf_num_df <- data.frame(vcf_num_t)Get sample names
sample <- row.names(vcf_num_df)Add sample into dataframe
vcf_num_df <- data.frame (sample, vcf_num_df)
#Check working directory
getwd()## [1] "/Users/austineastmure/Desktop/comp_bio/Final Project"#Save the cvs
write.csv(vcf_num_df,
file = "vcf_num_df.csv",
row.names = F)
#Confirm presence of file
list.files()## [1] "07-mean_imputation.html"
## [2] "07-mean_imputation.Rmd"
## [3] "08-PCA_worked.html"
## [4] "08-PCA_worked.Rmd"
## [5] "09-PCA_worked_example-SNPs-part1.html"
## [6] "09-PCA_worked_example-SNPs-part1.Rmd"
## [7] "10-PCA_worked_example-SNPs-part2.html"
## [8] "10-PCA_worked_example-SNPs-part2.Rmd"
## [9] "1000_genomes_snps"
## [10] "1000genomes_people_info2-1.csv"
## [11] "15.31093190-31333190.ALL.chr15_GRCh38.genotypes.20170504.vcf.gz"
## [12] "all_loci-1.vcf"
## [13] "all_loci.vcf"
## [14] "bird_snps_remove_NAs.html"
## [15] "Final Project.Rproj"
## [16] "Final_Project_AE.html"
## [17] "Final_Project_AE.Rmd"
## [18] "Final_Project_Script.R"
## [19] "final_report_AE.html"
## [20] "final_report_AE.Rmd"
## [21] "final_report_template.Rmd"
## [22] "removing_fixed_alleles.Rmd"
## [23] "rsconnect"
## [24] "Script.R"
## [25] "SNPs_cleaned.csv"
## [26] "Software_Checkpoint_Loading_VCF_file_into_R.R"
## [27] "transpose_VCF_data.html"
## [28] "transpose_VCF_data.Rmd"
## [29] "vcf_noNA.csv"
## [30] "vcf_num_df.csv"
## [31] "vcf_num_df2.csv"
## [32] "vcf_num.csv"
## [33] "walsh2017morphology.csv"##Clean the data
Merge data with population meta data
#Load population metadata
pop_meta <- read.csv(file = "1000genomes_people_info2-1.csv")
#Merge meta data with SNP data
##Make sure column "sample" appears in meta data and SNP data
names(pop_meta)## [1] "pop" "super_pop" "sample" "sex" "lat" "lng"names(vcf_num_df)[1:10]## [1] "sample" "X1" "X2" "X3" "X4" "X5" "X6" "X7"
## [9] "X8" "X9"##Merge the two sets of data
vcf_num_df2 <- merge(pop_meta,vcf_num_df, by = "sample")Quality Assurance: Check dimensions of data before and after the merge
nrow(vcf_num_df) == nrow(vcf_num_df2)## [1] TRUECheck the names of the new dataframe
names(vcf_num_df2)[1:15]## [1] "sample" "pop" "super_pop" "sex" "lat" "lng"
## [7] "X1" "X2" "X3" "X4" "X5" "X6"
## [13] "X7" "X8" "X9"Check working directory
getwd()## [1] "/Users/austineastmure/Desktop/comp_bio/Final Project"Save the csv
write.csv(vcf_num_df2, file = "vcf_num_df2.csv", row.names = F)
#Check for presence of file
list.files()## [1] "07-mean_imputation.html"
## [2] "07-mean_imputation.Rmd"
## [3] "08-PCA_worked.html"
## [4] "08-PCA_worked.Rmd"
## [5] "09-PCA_worked_example-SNPs-part1.html"
## [6] "09-PCA_worked_example-SNPs-part1.Rmd"
## [7] "10-PCA_worked_example-SNPs-part2.html"
## [8] "10-PCA_worked_example-SNPs-part2.Rmd"
## [9] "1000_genomes_snps"
## [10] "1000genomes_people_info2-1.csv"
## [11] "15.31093190-31333190.ALL.chr15_GRCh38.genotypes.20170504.vcf.gz"
## [12] "all_loci-1.vcf"
## [13] "all_loci.vcf"
## [14] "bird_snps_remove_NAs.html"
## [15] "Final Project.Rproj"
## [16] "Final_Project_AE.html"
## [17] "Final_Project_AE.Rmd"
## [18] "Final_Project_Script.R"
## [19] "final_report_AE.html"
## [20] "final_report_AE.Rmd"
## [21] "final_report_template.Rmd"
## [22] "removing_fixed_alleles.Rmd"
## [23] "rsconnect"
## [24] "Script.R"
## [25] "SNPs_cleaned.csv"
## [26] "Software_Checkpoint_Loading_VCF_file_into_R.R"
## [27] "transpose_VCF_data.html"
## [28] "transpose_VCF_data.Rmd"
## [29] "vcf_noNA.csv"
## [30] "vcf_num_df.csv"
## [31] "vcf_num_df2.csv"
## [32] "vcf_num.csv"
## [33] "walsh2017morphology.csv"Omit invariant features
#Load invar_omit() function
invar_omit <- function(x){
cat("Dataframe of dim",dim(x), "processed...\n")
sds <- apply(x, 2, sd, na.rm = TRUE)
i_var0 <- which(sds == 0)
cat(length(i_var0),"columns removed\n")
if(length(i_var0) > 0){
x <- x[, -i_var0]
}
return(x)
}Check which columns have invariant data
names(vcf_num_df2)[1:10]## [1] "sample" "pop" "super_pop" "sex" "lat" "lng"
## [7] "X1" "X2" "X3" "X4"#We don't want to put these columns through invar_omit(). We'll use negative indexing to skip
##New data frame to store output, no metadata
vcf_noinvar <- vcf_num_df2
##Run invar_omit() on numeric data
vcf_noinvar[, -c(1:6)] <- invar_omit(vcf_noinvar[, -c(1:6)])## Dataframe of dim 2504 8146 processed...
## 2055 columns removed#add back metadata
vcf_noinvar <- data.frame(vcf_num_df2[, c(1:6)], vcf_noinvar)Create an object to store the number of invariant columns removed
my_meta_N_invar_cols <- 2055Remove low-quality data
#Load find_NAs()
find_NAs <- function(x){
NAs_TF <- is.na(x)
i_NA <- which(NAs_TF == TRUE)
N_NA <- length(i_NA)
return(i_NA)
}
#for() loop to search for NAs
#N_rows
#number of rows (individuals)
N_rows <- nrow(vcf_noinvar)
#N_NA
#vector to hold output (number of NAs)
N_NA <- rep(x = 0, times = N_rows)
#N_SNPs
#total number of columns (SNPs)
N_SNPs <-ncol(vcf_noinvar)
cat("This may take a minute...")## This may take a minute...#For the for() loop
for(i in N_rows){
#For each row, find the location of
#NAs with find_NAs()
i_NA <-find_NAs(vcf_noinvar[i,])
#then determine how many NAs with
#length()
N_NA_i <- length(i_NA)
#then save the output to
#storage vector
N_NA[i] <- N_NA_i
}
#Check if any row has >50% NAs
#It will probably be 0 for 1000 genomes data
cutoff50 <- N_SNPs*0.5
percent_NA <- N_NA/N_SNPs*100
any(percent_NA > 50)## [1] FALSEWhat is the average number of NAs per row?
mean(percent_NA)## [1] 0Save the mean percent of NAs per row
my_meta_N_invar_rows <- mean(percent_NA)Imputation of NAs
Load imputation function
mean_imputation <- function(df){
cat("This may take some time...")
n_cols <- ncol(df)
for (i in n_cols) {
column_i <- df[, i]
#Get the mean of the current column
mean_i <- mean(column_i, na.rm = TRUE)
#Get the NAs in the current column
NAs_i <- which(is.na(column_i))
#report the number of NAs
N_NAs <- length(NAs_i)
#replace the NAs in the current column
column_i[NAs_i] <- mean_i
#replace the original column with the
#update columns
df[, i] <- column_i
}
return(df)
}
#run the functionWe’ll only run this on numeric columns
names(vcf_noinvar)[1:10]## [1] "sample" "pop" "super_pop" "sex" "lat"
## [6] "lng" "sample.1" "pop.1" "super_pop.1" "sex.1"#make a new copy of the data
vcf_noNA <- vcf_noinvar
vcf_noNA[, -c(1:6)] <- mean_imputation(vcf_noinvar[, -c(1:6)])## This may take some time...#this may take some timePrepare for PCA
Scale the data
#new copy of data
vcf_scaled <- vcf_noNA
#
vcf_scaled[, -c(1:12)] <- scale(vcf_noNA[,-c(1:12)])Run the PCA
vcf_pca <- prcomp(vcf_scaled[,-c(1:12)])PCA diagnostics
Examine the default screeplot
screeplot(vcf_pca)
Calculate explained variation
#Load PCA variation function
PCA_variation <- function(pca_summary, PCs = 2){
var_explained <- pca_summary$importance[2,1:PCs]*100
var_explained <- round(var_explained,3)
return(var_explained)
}Extract PCA variation data and calculate percentage variation
#Get summary information
vcf_pca_summary <- summary(vcf_pca)
#Extract raw variation data. You will probably need to consider >100 PCs
var_out <- PCA_variation(vcf_pca_summary, PCs = 2000)Calculate the cut off for the rule of thumb
#number of dimensions in the data
N_columns <- ncol(vcf_scaled)
#The value of the cutoff
cut_off <- 1/N_columns*100
#Calculate the number of PCs which exceed the cut off
# which values below the cutoff
i_cut_off <- which(var_out < cut_off)
#what is first value below cutoff?
i_cut_off <- min(i_cut_off)Save the first value below the cutoff
my_meta_N_meanNA_rowsPCs <- i_cut_off
#Extract the amount of variation explained by the first 3 PCs
my_meta_var_PC123 <- var_out[c(1,2,3)]Plot Percentage variation
#make barplot
barplot(var_out,
main = "Percent variation (%) Scree plot",
ylab = "Percent variation (%) Explained",
names.arg = 1:length(var_out))
abline(h = cut_off, col = 2, lwd = 2)
abline(v = i_cut_off)
legend("topright",
col = c(2,1),
lty = c(1,1),
legend = c("Vertical line: cutoff",
"Horizontal line: 1st value below cut off"))
Plot cumulative percentage variatio
Another way to look at this is to calculate the cumulative amount of variation explained Point 1 on the plot= PC1s variation explained. Point 2 = PC1 + PC2 Point 3 = PC1 + PC2 + PC3 and so on
cumulative_variation <- cumsum(var_out)
plot(cumulative_variation, type = "l")
Plot PCA Results
Note that if the arrows aren’t plotted, it is not a biplot
Calculate scores
Get the scores:
#call vegan::scores()
vcf_pca_scores <- vegan::scores(vcf_pca)Combine the scores with the species information into a dataframe
# call data.frame()
vcf_pca_scores2 <- data.frame(super_pop = vcf_noNA$super_pop, vcf_pca_scores)Look for information on the variation explained by the PCs
my_meta_var_PC123[1]## PC1
## 2.087my_meta_var_PC123[2]## PC2
## 1.655my_meta_var_PC123[3]## PC3
## 1.075Plot the results
Plot PC1 versus PC2
#Plot the scores, with super population color-coded
#make color and shape = "super_pop"
ggpubr::ggscatter(data = vcf_pca_scores2,
y = "PC2",
x = "PC1",
color = "super_pop",
shape = "super_pop",
main = "PCA_Scatterplot",
xlab = "PC1 (2% of variation)",
ylab = "PC2 (1.7% of variation)")
Note how in the plot the amount of variation is explained by each PC is
shown in the axis labels
Plot PC2 versus PC3
#Plot the scores, with super population color-coded
#make color and shape = "super_pop"
ggpubr::ggscatter(data = vcf_pca_scores2,
y = "PC3",
x = "PC2",
color = "super_pop",
shape = "super_pop",
main = "PCA_Scatterplot",
xlab = "PC2 (1.7% of variation)",
ylab = "PC3 (1% of variation)")
Note how in the plot the amount of variation is explained by each PC is
shown in the axis labels
Plot PC1 versus PC3
#Plot the scores, with super population color-coded
#make color and shape = "super_pop"
ggpubr::ggscatter(data = vcf_pca_scores2,
y = "PC3",
x = "PC1",
color = "super_pop",
shape = "super_pop",
main = "PCA_Scatterplot",
xlab = "PC1 (2% of variation)",
ylab = "PC3 (1% of variation)")
Note how in the plot the amount of variation is explained by each PC is
shown in the axis labels