library(vcfR)##
## ***** *** vcfR *** *****
## This is vcfR 1.13.0
## browseVignettes('vcfR') # Documentation
## citation('vcfR') # Citation
## ***** ***** ***** *****#Download the needed packages
library(vegan)## Loading required package: permute## Loading required package: lattice## This is vegan 2.6-4library(ggplot2)
library(ggpubr)#Check working directory
getwd()## [1] "/Users/ishashah/Downloads/comp_bio_final_project"#Confirm the files in the list
list.files()## [1] "1000genomes_people_info2-1.csv"
## [2] "2.136483646-136733646.ALL.chr2_GRCh38.genotypes.20170504 2.vcf"
## [3] "rsconnect"
## [4] "Shah_FinalProject_Workflow.html"
## [5] "Shah_FinalProject_Workflow.Rmd"
## [6] "vcf_num_df.csv"
## [7] "vcf_num_df2.csv"
## [8] "vcf_num.csv"#Load VCF data
my_vcf <- "2.136483646-136733646.ALL.chr2_GRCh38.genotypes.20170504 2.vcf"
vcf <- vcfR::read.vcfR(my_vcf, convertNA = T)## Scanning file to determine attributes.
## File attributes:
## meta lines: 130
## header_line: 131
## variant count: 7896
## column count: 2513
##
Meta line 130 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
## Character matrix gt rows: 7896
## Character matrix gt cols: 2513
## skip: 0
## nrows: 7896
## 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: 7896
## All variants processed#Convert raw VCF file to genotype scores
vcf_num <- vcfR::extract.gt(vcf,
element = "GT",
IDtoRowNames = F,
as.numeric = T,
convertNA = T)Save the CSV
write.csv(vcf_num, file = "vcf_num.csv", row.names = F)#Confirm presence of file
list.files()## [1] "1000genomes_people_info2-1.csv"
## [2] "2.136483646-136733646.ALL.chr2_GRCh38.genotypes.20170504 2.vcf"
## [3] "rsconnect"
## [4] "Shah_FinalProject_Workflow.html"
## [5] "Shah_FinalProject_Workflow.Rmd"
## [6] "vcf_num_df.csv"
## [7] "vcf_num_df2.csv"
## [8] "vcf_num.csv"#transpose original VCF orientation to R df orientation
vcf_num_t <- t(vcf_num)
#make into a df
vcf_num_df <- data.frame(vcf_num_t)
#get sample name
sample <- row.names(vcf_num_df)
#add sample into df
vcf_num_df <- data.frame(sample,vcf_num_df)getwd## function ()
## .Internal(getwd())
## <bytecode: 0x7f9a892afd90>
## <environment: namespace:base>write.csv(vcf_num_df, file = "vcf_num_df.csv", row.names = F)list.files()## [1] "1000genomes_people_info2-1.csv"
## [2] "2.136483646-136733646.ALL.chr2_GRCh38.genotypes.20170504 2.vcf"
## [3] "rsconnect"
## [4] "Shah_FinalProject_Workflow.html"
## [5] "Shah_FinalProject_Workflow.Rmd"
## [6] "vcf_num_df.csv"
## [7] "vcf_num_df2.csv"
## [8] "vcf_num.csv"load population data
pop_meta <- read.csv(file = "1000genomes_people_info2-1.csv")merge data with 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 data sets
vcf_num_df2 <- merge(pop_meta, vcf_num_df, by = "sample")#check dimensions before and after merge
nrow(vcf_num_df)== nrow(vcf_num_df2)## [1] TRUE#check names of 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 wd and also files
getwd()## [1] "/Users/ishashah/Downloads/comp_bio_final_project"write.csv(vcf_num_df2, file = "vcf_num_df2.csv", row.names = F)
list.files()## [1] "1000genomes_people_info2-1.csv"
## [2] "2.136483646-136733646.ALL.chr2_GRCh38.genotypes.20170504 2.vcf"
## [3] "rsconnect"
## [4] "Shah_FinalProject_Workflow.html"
## [5] "Shah_FinalProject_Workflow.Rmd"
## [6] "vcf_num_df.csv"
## [7] "vcf_num_df2.csv"
## [8] "vcf_num.csv"#omit invariant features
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 character data
names(vcf_num_df2)[1:10]## [1] "sample" "pop" "super_pop" "sex" "lat" "lng"
## [7] "X1" "X2" "X3" "X4"#new df to store output
vcf_noinvar <- vcf_num_df2
vcf_noinvar[, -c(1:6)] <- invar_omit(vcf_noinvar[, -c(1:6)])## Dataframe of dim 2504 7896 processed...
## 2117 columns removed#storing number of invar columns removed
my_meta_N_invar_cols <- 2117#removing low quality data
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’s
N_rows <- nrow(vcf_noinvar)
#number of rows (individuals)
N_NA <- rep(x = 0, times = N_rows)
#vector to hold output (number of NA)
N_SNPs <- ncol(vcf_noinvar)
#total number of columns (SNP)
for(i in 1:N_rows){
#find location of NAs
i_NA <- find_NAs(vcf_noinvar[i,])
#determine NA's with length
N_NA_i <- length(i_NA)
#save output to storage vector
N_NA[i] <- N_NA_i
}#check for less than 50% NA’s
cutoff50 <- N_SNPs*0.5
percent_NA <- N_NA/N_SNPs*100
any(percent_NA>50)## [1] FALSE#average number of NA’s per row
mean(percent_NA)## [1] 0#save mean percent of NA’s per row
my_meta_N_meanNA_rows <- mean(percent_NA)#load imputation function for NA’s in case
mean_imputation <- function(df){
n_cols <- ncol(df)
for(i in 1:n_cols){
#get current column
column_i <- df[,i]
#get mean of current column
mean_i <- mean(column_i, na.rm = TRUE)
#get NA's in the current column
NAs_i <- which(is.na(column_i))
#report number of NA's
N_NAs <- length(NAs_i)
#replace NA's in current column
column_i[NAs_i] <- mean_i
#replace original column with updated
df[,i <- column_i]
}
return(df)
}#run on numeric columns
names(vcf_noinvar)[1:10]## [1] "sample" "pop" "super_pop" "sex" "lat" "lng"
## [7] "X1" "X2" "X3" "X4"#new copy of data
vcf_noNA <- vcf_noinvar
vcf_noNA[,-c(1:6)] <- mean_imputation(vcf_noinvar[, -c(1:6)])##Prepare for PCA # Scale the data
vcf_scaled <- vcf_noNA
#scale
vcf_scaled[,-c(1:6)] <- scale(vcf_noNA[,-c(1:6)])#run PCA
vcf_pca <- prcomp(vcf_scaled[,-c(1:6)])##PCA diagnostics # Examine defauly 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)}#get summary info
vcf_pca_summary <- summary(vcf_pca)#extract raw variation data
var_out <- PCA_variation(vcf_pca_summary, PCs = 500)#caclulate cut off
#number of dimensions in the data
N_columns <- ncol(vcf_scaled)
#value of cutoff
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 Inf#extract amount of variation explained by first 3 PC’s
N_meanNA_rowsPCs <- i_cut_off
var_PC123 <- var_out[c(1,2,3)]plot percent variation
barplot(var_out,
main = "Percent variation (%) Scree Plot",
ylab = "Percent variation (%) explained",
xlab = "PCs",
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 cummulative percentage variation of PCs
cumulative_variation <- cumsum(var_out)
plot(cumulative_variation,
type = "l",
main = "Cumulative variation (%) Plot",
ylab = "Cumulative variation (%) explained",
xlab = "PCs")
Plot PCA results
Use the vgan package to get PCA scores & check variation by first 2 PCs
vcf_pca_scores <- vegan::scores(vcf_pca)
vcf_pca_scores2 <- data.frame(population = vcf_noNA$super_pop, vcf_pca_scores)
var_PC123[1]## PC1
## 2.959var_PC123[2]## PC2
## 2.051var_PC123[3]## PC3
## 1.949#compare PC1 & PC2 in a scatterplot color coded by super population
ggpubr::ggscatter(data = vcf_pca_scores2,
y = "PC2",
x = "PC1",
color = "population",
shape = "population",
main = "PCA Scatterplot",
xlab = "PC1 (3.0% of variation)",
ylab = "PC2 (2.05% of variation)")
#Compare PC2 & PC3 in a scatterplot color coded by super population
ggpubr::ggscatter(data = vcf_pca_scores2,
y = "PC3",
x = "PC2",
color = "population",
shape = "population",
main = "PCA Scatterplot",
xlab = "PC2 (2.05% of variation)",
ylab = "PC3 (1.95% of variation)")
#Compare PC1 & PC3 in a scatterplot color coded by super population
ggpubr::ggscatter(data = vcf_pca_scores2,
y = "PC3",
x = "PC1",
color = "population",
shape = "population",
main = "PCA Scatterplot",
xlab = "PC1 (3.0% of variation)",
ylab = "PC3 (1.95% of variation)")
#3D scatter plot
colors_use <- as.numeric(factor(vcf_pca_scores2$population))
scatterplot3d::scatterplot3d(
x = vcf_pca_scores2$PC1,
y = vcf_pca_scores2$PC2,
z = vcf_pca_scores2$PC3,
color = colors_use,
xlab = "PC1 (3.0%)",
ylab = "PC2 (2.05%)",
zlab = "PC3 (1.95%)")