Workflow Final Project
Elizabeth Taylor
2022-12-09
Preliminaries
Load the vcfR 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
21.31658131-31898131.ALL.chr21_GRCh38.genotypes.20170504.vcf.gz.
getwd()## [1] "/Users/elizabethtaylor/1540/Final Project"list.files(pattern = "vcf")## [1] "21.31658131-31898131.ALL.chr21_GRCh38.genotypes.20170504.vcf.gz"
## [2] "vcf_num_df.csv"
## [3] "vcf_num_df2.csv"
## [4] "vcf_num.csv"
## [5] "vcf_scaled.csv"Data preparation
Loading the Data
Load vcf file from “21.31658131-31898131.ALL.chr21_GRCh38.genotypes.20170504.vcf.gz” to an object called vcf so that it can be used in R. The working directory must be the same as where the .vcf file is saved.
my_vcf <- "21.31658131-31898131.ALL.chr21_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: 7102
## column count: 2513
##
Meta line 130 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
## Character matrix gt rows: 7102
## Character matrix gt cols: 2513
## skip: 0
## nrows: 7102
## 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: 7102
## All variants processedExtract numeric genotype scores (allele counts)
The vcfR::extract.gt() function is acting on the data stored in vcf to retrieve the numeric genotype scores and save them to the object vcf_num
vcf_num <- vcfR::extract.gt(vcf,
element = "GT",
IDtoRowNames = F,
as.numeric = T,
convertNA = T)
#return.alleles = F)Save the CSV
write.csv(vcf_num, file = "vcf_num.csv", row.names = F)Check the directory for the new csv file.
list.files()## [1] "1000genomes_people_info2.csv"
## [2] "21.31658131-31898131.ALL.chr21_GRCh38.genotypes.20170504.vcf.gz"
## [3] "Final Project.Rproj"
## [4] "final_report_template.Rmd"
## [5] "Final_Report.Rmd"
## [6] "rsconnect"
## [7] "vcf_num_df.csv"
## [8] "vcf_num_df2.csv"
## [9] "vcf_num.csv"
## [10] "vcf_scaled.csv"
## [11] "Workflow_Final_Project.html"
## [12] "Workflow_Final_Project.Rmd"Transpose Data
Transpose the data so that the genotype scores no longer have SNPs in columns and samples in rows
vcf_num_t <- t(vcf_num) Turn the vcf_num_t matrix into a dataframe and save it to the object vcf_num_df
vcf_num_df <- data.frame(vcf_num_t) Get the names of the samples (each person)
sample <- row.names(vcf_num_df)Add sample into the dataframe
vcf_num_df <- data.frame(sample, vcf_num_df)###Save new csv file Check directory and save the new csv for the dataframe
getwd()## [1] "/Users/elizabethtaylor/1540/Final Project"write.csv(vcf_num_df, file = "vcf_num_df.csv", row.names = F)Verify csv exists
list.files()## [1] "1000genomes_people_info2.csv"
## [2] "21.31658131-31898131.ALL.chr21_GRCh38.genotypes.20170504.vcf.gz"
## [3] "Final Project.Rproj"
## [4] "final_report_template.Rmd"
## [5] "Final_Report.Rmd"
## [6] "rsconnect"
## [7] "vcf_num_df.csv"
## [8] "vcf_num_df2.csv"
## [9] "vcf_num.csv"
## [10] "vcf_scaled.csv"
## [11] "Workflow_Final_Project.html"
## [12] "Workflow_Final_Project.Rmd"Clean the data
Load the data and merge the data with the population meta data in another csv file
pop_meta <- read.csv(file = "1000genomes_people_info2.csv")check for “sample” 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 data sets and check the dimesions from before and after
vcf_num_df2 <- merge(pop_meta, vcf_num_df, by = "sample")
nrow(vcf_num_df) == nrow(vcf_num_df2)## [1] TRUE###Save the csv file Check the names of the new df and save the new dataframe as a csv file
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"getwd()## [1] "/Users/elizabethtaylor/1540/Final Project"write.csv(vcf_num_df2, file = "vcf_num_df2.csv", row.names = F)Check for csv file
list.files()## [1] "1000genomes_people_info2.csv"
## [2] "21.31658131-31898131.ALL.chr21_GRCh38.genotypes.20170504.vcf.gz"
## [3] "Final Project.Rproj"
## [4] "final_report_template.Rmd"
## [5] "Final_Report.Rmd"
## [6] "rsconnect"
## [7] "vcf_num_df.csv"
## [8] "vcf_num_df2.csv"
## [9] "vcf_num.csv"
## [10] "vcf_scaled.csv"
## [11] "Workflow_Final_Project.html"
## [12] "Workflow_Final_Project.Rmd"Handling invariant features
Omit the invariant features by creating an 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 for the columns that have character data. Those columns no not get run with invar_omit(), so creat a new dataframe without those columns.
names(vcf_num_df2)[1:10]## [1] "sample" "pop" "super_pop" "sex" "lat" "lng"
## [7] "X1" "X2" "X3" "X4"vcf_noinvar <- vcf_num_df2Run invar_omit() on the numeric data
#vcf_noinvar[, -c(1:6)] <- invar_omit(vcf_noinvar[, -c(1:6)])
dim(vcf_num_df2)## [1] 2504 7108vcf_noinvar <- vcf_noinvar[, -c(1:6)]
vcf_noinvar <- invar_omit(vcf_noinvar)## Dataframe of dim 2504 7102 processed...
## 1846 columns removedvcf_noinvar <- data.frame(vcf_num_df2[, c(1:6)], vcf_noinvar)
dim(vcf_noinvar)## [1] 2504 5262Store the number of invariant columns removed in an object
my_meta_N_invar_cols <- 1846Remove low-quality data
Search for NAs
find_NAs <- function(x){
NAs_TF <- is.na(x)
i_NA <- which(NAs_TF == TRUE)
N_NA <- length(i_NA)
return(i_NA)
}Use find_NAs inside of a for loop
#number of rows (individuals)
N_rows <-nrow(vcf_noinvar)
#number of NAs
N_NA <- rep(x=0, times = N_rows)
#number of columns (snps)
N_SNPs <- ncol(vcf_noinvar)
for(i in 1:N_rows){
i_NA <- find_NAs(vcf_noinvar[i,])
#how many NAs in each row
N_NA_i <- length(i_NA)
N_NA[i] <- N_NA_i
}###Evaluating NA information Check if any row has >50% NAs
cutoff50 <- N_SNPs*0.5
percent_NA <- N_NA/N_SNPs*100
any(percent_NA > 50)## [1] FALSEAverage number of NAs per row (0)
mean(percent_NA)## [1] 0Save the mean percent of NAs per row
my_meta_N_meanNA_rows <- mean(percent_NA)###Mean imputation for NAs Load mean imputation function
mean_imputation <- function(df){
n_cols <- ncol(df)
for(i in 1:n_cols){
#current column
column_i <- df[, i]
#mean of that column
mean_i <- mean(column_i, na.rm = TRUE)
#NAs in current column
NAs_i <- which(is.na(column_i))
#number of NAs
N_NAs <- length(NAs_i)
#replace NAs in current column
column_i[NAs_i] <- mean_i
#replace column with updated column
df[, i] <- column_i
}
return(df)
}Check to make sure you are running mean_imputation() on numeric columns
names(vcf_noinvar)[1:10]## [1] "sample" "pop" "super_pop" "sex" "lat" "lng"
## [7] "X1" "X2" "X4" "X5"Perform mean imputation on a new copy of the data
vcf_noNA <- vcf_noinvar
vcf_noNA[, -c(1:6)] <- mean_imputation(vcf_noinvar[, -c(1:6)])###Prepare for PCA Scale the data (center around 0 and scale by standard deviation) Only perform on SNPs columns
#new copy of data
vcf_scaled <- vcf_noNA
#scale
vcf_scaled[, -c(1:6)] <- scale(vcf_noNA[, -c(1:6)])
write.csv(vcf_scaled, file = "vcf_scaled.csv", row.names = F)###Run PCA
vcf_pca <- prcomp(vcf_scaled[, -c(1:6)])###PCA diagnostics Default scree plot
screeplot(vcf_pca)
Calculate explained 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)
}PCA summary information
vcf_pca_summary <- summary(vcf_pca)Extract raw variation data
var_out <- PCA_variation(vcf_pca_summary, PCs = 500)Calculate the cut off
#number of dimensions in the data
N_columns <- ncol(vcf_scaled)
#value of the cut off using rull of thumb
cut_off <- 1/N_columns*100
#PCs below the cut off
i_cut_off <- which(var_out < cut_off)
#first value below the cut off
i_cut_off <- min(i_cut_off)## Warning in min(i_cut_off): no non-missing arguments to min; returning Inf#first value before the cut off
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
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 cutoff"))
###Plot cumulative percentage variation Calculate the cumulative amount
of variation explained
cumulative_variation <-cumsum(var_out)
plot(cumulative_variation, type = "l")
###Plot PCA results Get the scores and combine them with the species
information into a dataframe
vcf_pca_scores <- vegan::scores(vcf_pca)
vcf_pca_scores2 <- data.frame(super_pop = vcf_noNA$super_pop, vcf_pca_scores)Information on the variation explained by the PCs
my_meta_var_PC123[1]## PC1
## 1.927my_meta_var_PC123[2]## PC2
## 1.833my_meta_var_PC123[3]## PC3
## 1.671###Plot the results Plot the scores with super population color-coded PC1 versus PC2
ggpubr::ggscatter(data = vcf_pca_scores2,
y = "PC2",
x = "PC1",
color = "super_pop",
shape = "super_pop",
main = "PCA Scatterplot",
xlab = "PC1 (2.4% of variation)",
ylab = "PC2 (1.9% of variation)")
PC2 versus PC3
ggpubr::ggscatter(data = vcf_pca_scores2,
y = "PC3",
x = "PC2",
color = "super_pop",
shape = "super_pop",
main = "PCA Scatterplot",
xlab = "PC2 (1.9% of variation)",
ylab = "PC3 (1.8% of variation)")
PC1 versus PC3
ggpubr::ggscatter(data = vcf_pca_scores2,
y = "PC3",
x = "PC1",
color = "super_pop",
shape = "super_pop",
main = "PCA Scatterplot",
xlab = "PC1 (2.4% of variation)",
ylab = "PC3 (1.8% of variation)")