LOAD NECESSARY R PACKAGES
library(vegan)## Loading required package: permute## Loading required package: lattice## This is vegan 2.6-4library(ggplot2)
library(ggpubr)getwd()## [1] "/Users/parthpatel/Library/Mobile Documents/com~apple~CloudDocs/Pitt/BIOSC 1540/Final Project"list.files(pattern = "vcf")## [1] "11.11000-251000.ALL.chr11_GRCh38.genotypes.20170504.vcf"
## [2] "18.18000-258000.ALL.chr18_GRCh38.genotypes.20170504.vcf.gz"
## [3] "ALL.chr14_GRCh38.genotypes.20170504.vcf"
## [4] "vcf_num_df.csv"
## [5] "vcf_num_df2.csv"
## [6] "vcf_num.csv"
## [7] "vcf_scaled.csv"my_vcf <- "18.18000-258000.ALL.chr18_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: 6491
## column count: 2513
##
Meta line 130 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
## Character matrix gt rows: 6491
## Character matrix gt cols: 2513
## skip: 0
## nrows: 6491
## row_num: 0
##
Processed variant 1000
Processed variant 2000
Processed variant 3000
Processed variant 4000
Processed variant 5000
Processed variant 6000
Processed variant: 6491
## All variants processedvcf_num <- vcfR::extract.gt(vcf, element = "GT",
IDtoRowNames = F,
as.numeric = T,
convertNA = T)write.csv(vcf_num, file = "vcf_num.csv", row.names = F)#Confirm the presence of file
list.files()## [1] "1000genomes_people_info2.csv"
## [2] "11.11000-251000.ALL.chr11_GRCh38.genotypes.20170504.vcf"
## [3] "1540_final_project_Final_Report_template.pdf"
## [4] "1540_final_report_flowchart.pdf"
## [5] "1540_week14_PCA_SNP_workflow.pdf"
## [6] "18.18000-258000.ALL.chr18_GRCh38.genotypes.20170504.vcf.gz"
## [7] "ALL.chr14_GRCh38.genotypes.20170504.vcf"
## [8] "extra.Rmd"
## [9] "Final Project.Rproj"
## [10] "final_report_template.docx"
## [11] "final_report_template.html"
## [12] "final_report_template.Rmd"
## [13] "Final_Workflow.html"
## [14] "Final_Workflow.Rmd"
## [15] "my_snps"
## [16] "rsconnect"
## [17] "Screenshot 2022-12-03 at 3.49.33 PM.png"
## [18] "vcf_num_df.csv"
## [19] "vcf_num_df2.csv"
## [20] "vcf_num.csv"
## [21] "vcf_scaled.csv"Transpose original VCF orientation to R dataframe orientation
vcf_num_t <- t(vcf_num)Make into a dataframe
vcf_num_df <- data.frame(vcf_num_t)Get person (sample) name
sample <- row.names(vcf_num_df)Add sample info to dataframe
vcf_num_df <- data.frame(sample,
vcf_num_df)Save the CSV
write.csv(vcf_num_df,
file = "vcf_num_df.csv",
row.names = F)Confirm working directory and confirm presence of file
getwd()## [1] "/Users/parthpatel/Library/Mobile Documents/com~apple~CloudDocs/Pitt/BIOSC 1540/Final Project"list.files()## [1] "1000genomes_people_info2.csv"
## [2] "11.11000-251000.ALL.chr11_GRCh38.genotypes.20170504.vcf"
## [3] "1540_final_project_Final_Report_template.pdf"
## [4] "1540_final_report_flowchart.pdf"
## [5] "1540_week14_PCA_SNP_workflow.pdf"
## [6] "18.18000-258000.ALL.chr18_GRCh38.genotypes.20170504.vcf.gz"
## [7] "ALL.chr14_GRCh38.genotypes.20170504.vcf"
## [8] "extra.Rmd"
## [9] "Final Project.Rproj"
## [10] "final_report_template.docx"
## [11] "final_report_template.html"
## [12] "final_report_template.Rmd"
## [13] "Final_Workflow.html"
## [14] "Final_Workflow.Rmd"
## [15] "my_snps"
## [16] "rsconnect"
## [17] "Screenshot 2022-12-03 at 3.49.33 PM.png"
## [18] "vcf_num_df.csv"
## [19] "vcf_num_df2.csv"
## [20] "vcf_num.csv"
## [21] "vcf_scaled.csv"Merge data with population meta data Data was obtained from these sources: Data obtained from: https://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/integrated_call_samples_v3.20130502.ALL.panel
Information about the population codes can be found here: http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/README_populations.md
Load population META DATA
pop_meta <- read.csv(file = "1000genomes_people_info2.csv")Merge meta 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 sets of data
vcf_num_df2 <- merge(pop_meta,
vcf_num_df,
by = "sample")QA/AC: Check the dimensions before and after the merge
nrow(vcf_num_df) == nrow(vcf_num_df2)## [1] TRUECheck the names of the 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 working directory again
getwd()## [1] "/Users/parthpatel/Library/Mobile Documents/com~apple~CloudDocs/Pitt/BIOSC 1540/Final Project"Save the csv
write.csv(vcf_num_df2, file = "vcf_num_df2.csv", row.names = F)Confirm presence of file
list.files()## [1] "1000genomes_people_info2.csv"
## [2] "11.11000-251000.ALL.chr11_GRCh38.genotypes.20170504.vcf"
## [3] "1540_final_project_Final_Report_template.pdf"
## [4] "1540_final_report_flowchart.pdf"
## [5] "1540_week14_PCA_SNP_workflow.pdf"
## [6] "18.18000-258000.ALL.chr18_GRCh38.genotypes.20170504.vcf.gz"
## [7] "ALL.chr14_GRCh38.genotypes.20170504.vcf"
## [8] "extra.Rmd"
## [9] "Final Project.Rproj"
## [10] "final_report_template.docx"
## [11] "final_report_template.html"
## [12] "final_report_template.Rmd"
## [13] "Final_Workflow.html"
## [14] "Final_Workflow.Rmd"
## [15] "my_snps"
## [16] "rsconnect"
## [17] "Screenshot 2022-12-03 at 3.49.33 PM.png"
## [18] "vcf_num_df.csv"
## [19] "vcf_num_df2.csv"
## [20] "vcf_num.csv"
## [21] "vcf_scaled.csv"load invar_omit() function
invar_omit <- function(x){
cat("Dataframe of dim", dim(x), "processed... \n")
sds <- apply(x,2,sd, na.rm=T)
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 dataframe to store output
vcf_noinvar <- vcf_num_df2Run invar_omit on numeric data
vcf_noinvar[, -c(1:6)] <- invar_omit(vcf_noinvar[, -c(1:6)])## Dataframe of dim 2504 6491 processed...
## 1713 columns removedCreate an object to store the number of invariant columns removed
my_meta_N_invar_cols <- 1713find_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
N_rows <- nrow(vcf_noinvar)
N_NA <- rep(x = 0, times = N_rows)
N_SNPs <- ncol(vcf_noinvar)
cat("This may take a minute...")## This may take a minute...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
}Check if any row has >50% NAs it will probably be 0 for 1000 genomes data.
cutoff50 <- N_SNPs*0.5
percent_NA <- N_NA/N_SNPs*100
any(percent_NA > 50)## [1] FALSEWhat is the average number of NAs per row? It is probably 0 or close to 0
mean(percent_NA)## [1] 0.009933314Save the mean percent of NAs per row
my_meta_N_meanNA_rows <- mean(percent_NA)NAs are rare. This is doing an imputation just in case.
mean_imputation <- function(df){
cat("This may take some time...")
n_cols <- ncol(df)
for(i in 1:n_cols){
#get the current column
column_i <- df[, i]
#get the mean of the current column
mean_i <- mean(column_i, na.rm = TRUE)
#get the NAs in the current column
NAs_i <- which(is.na(column_i))
#report the number of NAs
N_NAs <- length(NAs_i)
#replace the NAs in the current column
column_i[NAs_i] <- mean_i
#replace the original column with the updated columns
df[, i] <- column_i
}
return(df)
}We will only run this 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)])## This may take some time...Scale the data
Many studies use centering around 0 and scaling by the standard deviation.
We only want to run this on our SNP 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)