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)getwd()## [1] "/Users/allenzhang/Desktop/BIOSC 1540 COMP BIO/FINAL PROJECT"list.files(pattern = "vcf")## [1] "17.45622108-45862108.ALL.chr17_GRCh38.genotypes.20170504.vcf.gz"
## [2] "ALL.chr17_GRCh38.genotypes.20170504 (1).vcf.gz"
## [3] "vcf_num_df2.csv"
## [4] "vcf_num.csv"#load the vcf data begin analyzing and cleaning the data
my_vcf <- "17.45622108-45862108.ALL.chr17_GRCh38.genotypes.20170504.vcf.gz"#load the vcf file
vcf <- vcfR::read.vcfR(my_vcf,
convertNA = T)## Scanning file to determine attributes.
## File attributes:
## meta lines: 130
## header_line: 131
## variant count: 7320
## column count: 2513
##
Meta line 130 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
## Character matrix gt rows: 7320
## Character matrix gt cols: 2513
## skip: 0
## nrows: 7320
## 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: 7320
## All variants processed#Convert Raw VCF file to genotype scores (allele counts) convert SNP data from categorical to numeric data dimensional reduction #get genotype score
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] "17.45622108-45862108.ALL.chr17_GRCh38.genotypes.20170504.vcf.gz"
## [3] "ALL.chr17_GRCh38.genotypes.20170504 (1).vcf.gz"
## [4] "FINAL PROJECT.Rproj"
## [5] "Final_Project_Analysis_of_1000_Genomes_Data_with_PCA.Rmd"
## [6] "rsconnect"
## [7] "scaled_data_for_PCA.csv"
## [8] "vcf_num_df2.csv"
## [9] "vcf_num.csv"
## [10] "vsc_num_df.csv"
## [11] "Zhang_FinalProject_Workflow.Rmd"it is now SNPs in rows and individuals in columns, so transpose the matrix and convert our data to a dataframe
vcf_num_t <- t(vcf_num)#make in to a dataframe
vcf_num_df <- data.frame(vcf_num_t)#get person names #add into dataframe add in the row names from the dataframe as a column that can be removed
sample <- row.names(vcf_num_df)
vcf_num_df <- data.frame(sample,
vcf_num_df)#check working directory #save the csv
getwd()## [1] "/Users/allenzhang/Desktop/BIOSC 1540 COMP BIO/FINAL PROJECT"write.csv(vcf_num_df,
file = "vsc_num_df.csv",
row.names = F)#confirm presence of file
list.files()## [1] "1000genomes_people_info2-1.csv"
## [2] "17.45622108-45862108.ALL.chr17_GRCh38.genotypes.20170504.vcf.gz"
## [3] "ALL.chr17_GRCh38.genotypes.20170504 (1).vcf.gz"
## [4] "FINAL PROJECT.Rproj"
## [5] "Final_Project_Analysis_of_1000_Genomes_Data_with_PCA.Rmd"
## [6] "rsconnect"
## [7] "scaled_data_for_PCA.csv"
## [8] "vcf_num_df2.csv"
## [9] "vcf_num.csv"
## [10] "vsc_num_df.csv"
## [11] "Zhang_FinalProject_Workflow.Rmd"#load population meta data the newly added sample
pop_meta <- read.csv(file = "1000genomes_people_info2-1.csv")#check column confirm both dataframes contain a column which are representative of the sample Merge data with population meta 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 set of data
vcf_num_df2 <- merge(pop_meta,
vcf_num_df,
by = "sample")#check the dimentsion check the simension before and after merge
nrow(vcf_num_df) == nrow(vcf_num_df2)## [1] TRUE#check the names of new dataframe check the names of new merged dataframe columns
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 #save the csv #confirm presence of file
getwd()## [1] "/Users/allenzhang/Desktop/BIOSC 1540 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] "17.45622108-45862108.ALL.chr17_GRCh38.genotypes.20170504.vcf.gz"
## [3] "ALL.chr17_GRCh38.genotypes.20170504 (1).vcf.gz"
## [4] "FINAL PROJECT.Rproj"
## [5] "Final_Project_Analysis_of_1000_Genomes_Data_with_PCA.Rmd"
## [6] "rsconnect"
## [7] "scaled_data_for_PCA.csv"
## [8] "vcf_num_df2.csv"
## [9] "vcf_num.csv"
## [10] "vsc_num_df.csv"
## [11] "Zhang_FinalProject_Workflow.Rmd"omit any invariant data to avoid issues when scaling data for PCA later
#load invar_omit() function
invar_omit <- function(x){
cat("Dataframe of dim", dim(x), "procedded...\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"#make sure we run this function only on columns with numeric data. keep using negative indexing
vcf_noinvar <- vcf_num_df2
#run invar_omit() on numeric data
#the first 6 are character data
vcf_noinvar[,-c(1:6)] <- invar_omit(vcf_noinvar[,-c(1:6)])## Dataframe of dim 2504 7320 procedded...
## 1725 columns removed#new object store the number
my_meta_N_invar_cols <- 1725
#1725 columns are removedpoor quality of sample may skew the results #load find_NAs() create a function to determine where NAs occur
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 find NAs present in each row of our dataframe
#number of rows (individuals)
N_rows <- nrow(vcf_noinvar)
#vector to hold output (number of NAs)
N_NA <- rep(x = 0, times = N_rows)
#total number of columns (SNPs)
N_SNPs <- ncol(vcf_noinvar)
cat("This may take a minute...")## This may take a minute...#for() loop
for(i in 1:N_rows){
#find the location of NAs with vcf_noinvar()
i_NA <- find_NAs(vcf_noinvar[i,])
#determine how many NAs with length()
N_NA_i <- length(i_NA)
#save the output to storage vector
N_NA[i] <- N_NA_i
}#check if any rwo has >50% NAs
cutoff50 <- N_SNPs * 0.5
percent_NA <- N_NA / N_SNPs*100
any(percent_NA > 50)## [1] FALSE#check the average number of NAs per row #save the mean
mean(percent_NA)## [1] 0#save the mean percent of NAs per row
my_meta_N_meanNA_rows <- mean(percent_NA)although it is rare to have NAs, we still need to do a mean imputation to make sure all NAs are removed #mean imputation
mean_imputation <- function(df){
cat("This may take some time ...")
n_cols <- ncol(df)
for(i in 1:n_cols){
#get current column number
column_i <-df[,i]
#get the mean of the column
mean_i <- mean(column_i, na.rm = TRUE)
#get indices of NAs in column
NAs_i <- which(is.na(column_i))
#report number of NAs
N_NAs <- length(NAs_i)
#replace the NAs in the column with the mean
column_i[NAs_i] <- mean_i
#replace original columns with that of the df
df[,i] <- column_i
}
return(df)
}only run this function on numeric columns #check for character data
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)])## This may take some time ...we have cleaned and processed the data, then scale the data so it can be used for PCA
#new copy of data
vcf_scaled <- vcf_noNA
#scale the data
vcf_scaled[,-c(1:6)] <- scale(vcf_noNA[,-c(1:6)])#Write Data
write.csv(vcf_scaled, file = "scaled_data_for_PCA.csv", row.names = F)