Final Project Workflow
Ethan Frank
2022-12-14
Preliminaries
Load the vcfR package and other packages with
library().
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.6-4library(ggplot2)
library(ggpubr)Make sure that your working directory is set to the location of the file.
setwd("/Users/ethanfrank/Desktop/R/Final Project")
getwd()## [1] "/Users/ethanfrank/Desktop/R/Final Project"list.files()## [1] "09-PCA_worked_example-SNPs-part1.Rmd"
## [2] "10-PCA_worked_example-SNPs-part2.Rmd"
## [3] "1000genomes_people_info2-1.csv"
## [4] "1540_final_project_Final_Report_template.pdf"
## [5] "1540_final_report_flowchart (1).pdf"
## [6] "1540_week14_PCA_SNP_workflow.pdf"
## [7] "17.20809577-21049577.ALL.chr17_GRCh38.genotypes.20170504.vcf.gz"
## [8] "Final Project Report.Rmd"
## [9] "Final Project Workflow.Rmd"
## [10] "final_report_template.Rmd"
## [11] "Final-Project-Workflow.docx"
## [12] "Final-Project-Workflow.Rmd"
## [13] "vcf_num_df.csv"
## [14] "vcf_num_df2.csv"
## [15] "vcf_num.csv"
## [16] "vcf_scaled.csv"list.files(pattern = "vcf")## [1] "17.20809577-21049577.ALL.chr17_GRCh38.genotypes.20170504.vcf.gz"
## [2] "vcf_num_df.csv"
## [3] "vcf_num_df2.csv"
## [4] "vcf_num.csv"
## [5] "vcf_scaled.csv"Extract Data from .vcf
Reads metadata and data from .vcf file and argument to convert “.” in data to NA. Creates R object out of .vcf data.
vcf <- vcfR::read.vcfR("17.20809577-21049577.ALL.chr17_GRCh38.genotypes.20170504.vcf.gz", convertNA = TRUE)## Scanning file to determine attributes.
## File attributes:
## meta lines: 130
## header_line: 131
## variant count: 8483
## column count: 2513
##
Meta line 130 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
## Character matrix gt rows: 8483
## Character matrix gt cols: 2513
## skip: 0
## nrows: 8483
## 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: 8483
## All variants processedExtract Genotype Scores
Extracts genotype data from R object and converts to numeric values. Stores numeric values into another R object.
vcf_num <- vcfR::extract.gt(vcf,
element = "GT",
IDtoRowNames = F,
as.numeric = T,
convertNA = T)Save New Data
Save the new data to .csv
write.csv(vcf_num, file="vcf_num.csv", row.names = F)
list.files(pattern = "vcf")## [1] "17.20809577-21049577.ALL.chr17_GRCh38.genotypes.20170504.vcf.gz"
## [2] "vcf_num_df.csv"
## [3] "vcf_num_df2.csv"
## [4] "vcf_num.csv"
## [5] "vcf_scaled.csv"Transpose Data
Flips rows and columns of data and stores into R object.
vcf_num_t <- t(vcf_num) Converts R object data into data frame.
vcf_num_df <- data.frame(vcf_num_t) Store row names into object.
sample <- row.names(vcf_num_df)Rewrite data frame with row names.
vcf_num_df <- data.frame(sample, vcf_num_df)Save vcf_num_df into .csv file.
write.csv(vcf_num_df, file = "vcf_num_df.csv", row.names = F)
list.files(pattern = "csv")## [1] "1000genomes_people_info2-1.csv" "vcf_num_df.csv"
## [3] "vcf_num_df2.csv" "vcf_num.csv"
## [5] "vcf_scaled.csv"Clean Data
Read .csv file from 1000genomes.
pop_meta <- read.csv(file = "1000genomes_people_info2-1.csv")List names of pop_meta.
names(pop_meta)## [1] "pop" "super_pop" "sample" "sex" "lat" "lng"List names of columns 1-10 of vcf_num_df.
names(vcf_num_df[1:10])## [1] "sample" "X1" "X2" "X3" "X4" "X5" "X6" "X7"
## [9] "X8" "X9"Merge pop_meta and vcf_num_df by “sample” and store in vcf_num_df2.
vcf_num_df2 <- merge(pop_meta, vcf_num_df, by = "sample")Check if row lengths are equal.
nrow(vcf_num_df) == nrow(vcf_num_df2)## [1] TRUEList names of columns 1-15 of vcf_num_df2.
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"Save the new data to .csv
write.csv(vcf_num_df2, file="vcf_num_df2.csv", row.names = F)
list.files(pattern = "csv")## [1] "1000genomes_people_info2-1.csv" "vcf_num_df.csv"
## [3] "vcf_num_df2.csv" "vcf_num.csv"
## [5] "vcf_scaled.csv"Omit Invariant Features
Create invar_omit function which takes a data frame object argument. Calculates standard deviation of each column. Stores each column which sd is 0 into i_var0 object. Removes columns which sd is 0. Returns manipulated data frame.
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)
} Run invar_omit function on vcf_num_df2 and store into vcf_noinvar.
vcf_noinvar <- vcf_num_df2
vcf_noinvar <- vcf_noinvar[, -c(1:6)]
vcf_noinvar <- invar_omit(vcf_noinvar)## Dataframe of dim 2504 8483 processed...
## 2041 columns removedRewrite data frame object with first 6 columns from vcf_num_df2.
vcf_noinvar <- data.frame(vcf_num_df2[, c(1:6)], vcf_noinvar)Store number of invariant columns in object.
N_of_invar_cols <- 2041Remove NAs
Creates a function with a data frame object as an argument that returns the number of NAs present in the data frame object.
find_NAs <- function(x) {
NAs_TF <- is.na(x)
i_NA <- which(NAs_TF == TRUE)
N_NA <- length(i_NA)
cat("Results:", N_NA, "NAs present\n.")
return(i_NA)
}Creates multiple objects. N_rows is the number of rows in vcf_noinvar. N_NA is a vector to hold the number of NAs. N_SNPs is the number of columns in vcf_noinvar.
N_rows <- nrow(vcf_noinvar)
N_NA <- rep(x = 0, times = N_rows)
N_SNPs <- ncol(vcf_noinvar)For loop loops through each row and finds location of NAs and stores them into i_NA. Finds the number of NAs and stores them in N_NA_i. Save output to our original vector N_NA.
# for (i in 1:N_rows) {
# i_NA <- find_NAs(vcf_noinvar[i,])
# N_NA_i <- length(i_NA)
# N_NA[i] <- N_NA_i
# }0 NAs were found so code was commented out for performance.
Cut the number of columns in half and store them in cutoff50. Stores percentage of NAs in percent_NA. Finds which NAs are above 50% threshhold.
# cutoff50 <- N_SNPs*0.5
# percent_NA <- N_NA/N_SNPs*100
# any(percent_NA > 50)
# mean(percent_NA)
# n_meanNA_rows <- mean(percent_NA)0 NAs were found so code was commented out for performance.
Mean Imputation
Create mean imputation function which takes a data frame object argument and replaces NAs with mean value from column.
mean_imputation <- function(x) {
cat("This may take some time...")
n_cols <- ncol(x)
for (i in 1:n_cols) {
column_i <- x[, i]
mean_i <- mean(column_i, na.rm = TRUE)
NAs_i <- which(is.na(column_i))
N_NAs <- length(NAs_i)
column_i[NAs_i] <- mean_i
x[, i] <- column_i
}
return(x)
}Run mean_imputation() on vcf_noinvar.
vcf_noNA <- vcf_noinvar
vcf_noNA[, -c(1:6)] <- mean_imputation(vcf_noinvar[, -c(1:6)])## This may take some time...Scale Data for PCA
Scale data to prepare for PCA
vcf_noNA <- vcf_noinvar
vcf_scaled <- vcf_noNA
vcf_scaled[, -c(1:6)] <- scale(vcf_noNA[, -c(1:6)])Save scaled data to .csv.
write.csv(vcf_scaled, file = "vcf_scaled.csv", row.names = F)
list.files(pattern = "csv")## [1] "1000genomes_people_info2-1.csv" "vcf_num_df.csv"
## [3] "vcf_num_df2.csv" "vcf_num.csv"
## [5] "vcf_scaled.csv"Run PCA
Run prcomp() on scaled data and store in vcf_pca.
vcf_pca <- prcomp(vcf_scaled[, -c(1:6)])Run screeplot() on vcf_pca.
screeplot(vcf_pca)
PCA Variation
Create PCA_Variation function to return % variation.
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)
}vcf_pca_summary <- summary(vcf_pca)
var_out <- PCA_variation(vcf_pca_summary, PCs=500)
N_columns <- ncol(vcf_scaled)
cut_off <- 1/N_columns*100
i_cut_off <- which(var_out < cut_off)
i_cut_off <- min(i_cut_off)## Warning in min(i_cut_off): no non-missing arguments to min; returning InfPlot PCA percent variation on screeplot.
N_meanNA_rowsPCs <- i_cut_off
var_PC123 <- var_out[c(1, 2, 3)]
barplot(var_out,
main = "Percent Variation",
ylab = "Percent Variation explained",
names.arg = 1:length(var_out))
abline(h = cut_off, col = 4, lwd = 2)
abline(v = i_cut_off)
legend("topright", col = c(1, 4), lty = c(1, 1),
legend = c("Vertical line: cutoff", "Horizontal line: 1st value below cutoff"))
Plot cumulative variation of var_out.
cumulative_variation <- cumsum(var_out)
plot(cumulative_variation, type = "l")
PCA Results
Extract scores from vcf_pca and store into vcf_pca_scores. Create data frame from vcf_pca_scores.
vcf_pca_scores <- vegan::scores(vcf_pca)
vcf_pca_scores2 <- data.frame(super_pop = vcf_noNA$super_pop, vcf_pca_scores)PC1 Variation
var_PC123[1]## PC1
## 3.417PC2 Variation
var_PC123[2]## PC2
## 2.55PC3 Variation
var_PC123[3]## PC3
## 1.928Plot PCA Results on Scatterplot
PC1 v PC2
ggpubr::ggscatter(data = vcf_pca_scores2,
x = "PC1",
y = "PC2",
color = "super_pop",
shape = "super_pop",
main = "PC1 v PC2 Scatterplot",
xlab = "PC1 3.417% variation",
ylab = "PC2 2.55% variation")
PC1 v PC3
ggpubr::ggscatter(data = vcf_pca_scores2,
x = "PC1",
y = "PC3",
color = "super_pop",
shape = "super_pop",
main = "PC1 v PC3 Scatterplot",
xlab = "PC1 3.417% variation",
ylab = "PC3 1.928% variation")
PC2 v PC3
ggpubr::ggscatter(data = vcf_pca_scores2,
x = "PC2",
y = "PC3",
color = "super_pop",
shape = "super_pop",
main = "PC2 v PC3 Scatterplot",
xlab = "PC2 2.55% variation",
ylab = "PC3 1.928% variation")