This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:
#loading R packages needed
library(vcfR)## Warning: package 'vcfR' was built under R version 4.2.2##
## ***** *** 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-2library(ggplot2)
library(ggpubr)#Confirming working directory and location of files
getwd()## [1] "C:/Users/arhat/Desktop/Computatational Biology/Final Project"list.files()## [1] "1000genomes_people_info2-1.csv"
## [2] "5.360294-600294.ALL.chr5_GRCh38.genotypes.20170504.vcf.gz"
## [3] "Final Project.Rproj"
## [4] "final_report_template.Rmd"
## [5] "Pradhan_Arhat_Final_Project.Rmd"
## [6] "vcf_num.csv"
## [7] "vcf_num_df.csv"
## [8] "vcf_num_df2.csv"
## [9] "vcf_scaled.csv"#Loading my SNP data
my_vcf <-"5.360294-600294.ALL.chr5_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: 8669
## column count: 2513
##
Meta line 130 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
## Character matrix gt rows: 8669
## Character matrix gt cols: 2513
## skip: 0
## nrows: 8669
## 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: 8669
## All variants processed#Converting raw VCF to genotype scores
vcf_num <- vcfR::extract.gt(vcf,
element = "GT",
IDtoRowNames = F,
as.numeric = T,
convertNA = T)#Saving as a CSV
write.csv(vcf_num, file = "vcf_num.csv", row.names = F)#checking if file was made
list.files(pattern = "csv")## [1] "1000genomes_people_info2-1.csv" "vcf_num.csv"
## [3] "vcf_num_df.csv" "vcf_num_df2.csv"
## [5] "vcf_scaled.csv"#Transpose VCF for R
vcf_num_t <- t(vcf_num)
#make into data frame
vcf_num_df <- data.frame(vcf_num_t)#Get sample names
sample <- row.names(vcf_num_df)
#add sample into df
vcf_num_df <- data.frame(sample, vcf_num_df)#check wd and then save as csv
getwd()## [1] "C:/Users/arhat/Desktop/Computatational Biology/Final Project"write.csv(vcf_num_df,
file = "vcf_num_df.csv",
row.names = F)#check for df csv file
list.files()## [1] "1000genomes_people_info2-1.csv"
## [2] "5.360294-600294.ALL.chr5_GRCh38.genotypes.20170504.vcf.gz"
## [3] "Final Project.Rproj"
## [4] "final_report_template.Rmd"
## [5] "Pradhan_Arhat_Final_Project.Rmd"
## [6] "vcf_num.csv"
## [7] "vcf_num_df.csv"
## [8] "vcf_num_df2.csv"
## [9] "vcf_scaled.csv"#merge data with population meta data
pop_meta <- read.csv(file = "1000genomes_people_info2-1.csv")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"#Merging the 2 sets of data
vcf_num_df2 <- merge(pop_meta,
vcf_num_df,
by = "sample")#check the number of rows and col
nrow(vcf_num_df) == nrow(vcf_num_df2)## [1] TRUE#Checking the names of the new df
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"#Checking wd
getwd()## [1] "C:/Users/arhat/Desktop/Computatational Biology/Final Project"#save the CSV
write.csv(vcf_num_df2, file = "vcf_num_df2.csv", row.names = F)#Omiting invariants
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]
}
## add return() with x in it
return(x)
}#Check which cols have char data
names(vcf_num_df2)[1:10]## [1] "sample" "pop" "super_pop" "sex" "lat" "lng"
## [7] "X1" "X2" "X3" "X4"#creating the invar df
vcf_noinvar <- vcf_num_df2
vcf_noinvar[, -c(1:6)] <- invar_omit(vcf_noinvar[, -c(1:6)])## Dataframe of dim 2504 8669 processed...
## 2252 columns removed#store all invar col removed
my_meta_N_invar_cols <- 2252#Remove low-quality data
#Creating function to find NA's
find_NAs <- function(x){
NAs_TF <- is.na(x)
i_NA <- which(NAs_TF == TRUE)
N_NA <- length(i_NA)
return(i_NA)
}N_rows <- nrow(vcf_noinvar)
N_NA <- rep(x =0, times = N_rows)
N_SNPs <- ncol(vcf_noinvar)for(i in 1:N_rows){
#finds the NAs in each row
i_NA <- find_NAs(vcf_noinvar[i,])
#determine how many NAs
N_NA_i <- length(i_NA)
#save the output
N_NA[i] <- N_NA_i
}#Checking if any row has >50% of NA’s
cutoff50 <- N_SNPs * 0.5
percent_NA <- N_NA/N_SNPs*100
any(percent_NA > 50)## [1] FALSEmy_meta_N_meanNA_rows<- mean(percent_NA)#Imputation of NA’s
mean_imputation <- function(df){
n_cols <- ncol(df)
for(i in 1:n_cols){
#get the current col
column_i <- df[,i]
#get the mean of the current col
mean_i <- mean(column_i, na.rm = TRUE)
#get the Nas in the current col
NAs_i <- which(is.na(column_i))
#report num of Nas
N_NAs <- length(NAs_i)
#replace the Nas in the current col
column_i[NAs_i] <-mean_i
#replce the original column with the means
df[,i] <- column_i
}
return(df)
}names(vcf_noinvar)[1:10]## [1] "sample" "pop" "super_pop" "sex" "lat" "lng"
## [7] "X1" "X2" "X3" "X4"#Runing the Imputation function
vcf_noNA <- vcf_noinvar
vcf_noNA[,-c(1:6)] <- mean_imputation(vcf_noinvar[,-c(1:6)])#Prepping for PCA
#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)])#Examine the results
# screeplot(vcf_pca)#Calculate explained variation
#pca varation 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 PCA summary info
# vcf_pca_summary <- summary(vcf_pca)#Extract raw variation data
# var_out <- PCA_variation(vcf_pca_summary, PCs = 500)#calculate the cut off for the rule of thumb
# 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)
#
# my_meta_N_meanNA_rowsPCs <- i_cut_off# 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 cut off"))#Plot cumulative percentage variation
# cumulative_variation <- cumsum(var_out)
# plot(cumulative_variation, type = "l")#Plot PCA results
# vcf_pca_scores <- vegan::scores(vcf_pca)
#
# vcf_pca_scores2 <- data.frame(super_pop = vcf_noNA$super_pop,
# vcf_pca_scores)
# my_meta_var_PC123[1]
# my_meta_var_PC123[2]
# my_meta_var_PC123[3]
# ```
#Plot the PC results #Plot PC1 vs PC2
# ggpubr::ggscatter(data = vcf_pca_scores2,
# y = "PC2",
# x = "PC1",
# color = "super_pop",
# shape = "super_pop",
# main = "PCA Scatterplot",
# xlab = "PC1 (3.1% of variation)",
# ylab = "PC2 (1.8% of variation)")#Plot PC2 vs PC3
# ggpubr::ggscatter(data = vcf_pca_scores2,
# y = "PC3",
# x = "PC2",
# color = "super_pop",
# shape = "super_pop",
# main = "PCA Scatterplot",
# xlab = "PC2 (1.8% of variation)",
# ylab = "PC3 (1.5% of variation)")#Plotting PC1 vs PC3
##ggpubr::ggscatter(data = vcf_pca_scores2,
# y = "PC1",
# x = "PC3",
# color = "super_pop",
# shape = "super_pop",
# main = "PCA Scatterplot",
# xlab = "PC1 (3.1% of variation)",
# ylab = "PC3 (1.5% of variation)")
###K-means cluster analysis
#my_meta_N_meanNA_rowsPCs