Prepare workspace

# Load dependencies
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-4
library(ggplot2)
library(ggpubr)
# Verify Working directory
getwd()
## [1] "C:/Users/Ian/OneDrive - University of Pittsburgh/Academic/Fall 2022/Comp bio/Final Project"
list.files()
## [1] "1000genomes_people_info2-1.csv" "Final Project.Rproj"           
## [3] "final_report_template.Rmd"      "ity1_clean_data.csv"           
## [5] "myData.vcf.gz"                  "rsconnect"                     
## [7] "Workflow.html"                  "Workflow.Rmd"

Load VCF File -> Dataframe

vcf <- read.vcfR("myData.vcf.gz",convertNA = T)
## Scanning file to determine attributes.
## File attributes:
##   meta lines: 130
##   header_line: 131
##   variant count: 6988
##   column count: 2513
## 
Meta line 130 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
##   Character matrix gt rows: 6988
##   Character matrix gt cols: 2513
##   skip: 0
##   nrows: 6988
##   row_num: 0
## 
Processed variant 1000
Processed variant 2000
Processed variant 3000
Processed variant 4000
Processed variant 5000
Processed variant 6000
Processed variant: 6988
## All variants processed
# Get genotype section and convert to numeric scores
vcf_num <- extract.gt(vcf, 
                      element = "GT", 
                      IDtoRowNames = F,
                      as.numeric = T,
                      convertNA = T)
#transpose and convert to dataframe
vcf_num_t <- t(vcf_num)
vcf_df <- data.frame(vcf_num_t)
# Add sample infor to data frame
sample <- row.names(vcf_df)
vcf_df <- data.frame(sample, vcf_df)

Clean Data

Merge with population data

pop_meta <- read.csv(file = "1000genomes_people_info2-1.csv")

#verify pop and vcf have a sample field
names(pop_meta)
## [1] "pop"       "super_pop" "sample"    "sex"       "lat"       "lng"
names(vcf_df[1:10])
##  [1] "sample" "X1"     "X2"     "X3"     "X4"     "X5"     "X6"     "X7"    
##  [9] "X8"     "X9"
#Merge!
vcf_full_df <- merge(pop_meta, vcf_df, by = "sample")

#Verify no data loss
nrow(vcf_full_df) == nrow(vcf_df)
## [1] TRUE
#Check fields for new df
names(vcf_full_df)[1:15]
##  [1] "sample"    "pop"       "super_pop" "sex"       "lat"       "lng"      
##  [7] "X1"        "X2"        "X3"        "X4"        "X5"        "X6"       
## [13] "X7"        "X8"        "X9"

Omit Invariant SNPs

#find and omit invariant columns
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)                      
}
#find character columns
names(vcf_full_df)[1:10] #first 6 are 
##  [1] "sample"    "pop"       "super_pop" "sex"       "lat"       "lng"      
##  [7] "X1"        "X2"        "X3"        "X4"
#run omit function of the df
#NOTE: the example in the pdf resulted in no cols being removed, I have instead used cbind to get the desired result
vcf_noinvar <- invar_omit(vcf_full_df[,-c(1:6)])
## Dataframe of dim 2504 6988 processed...
## 1721 columns removed
vcf_noinvar <- cbind(vcf_full_df[1:6],vcf_noinvar)
#store num removed
num_invar <- 1721
total_snps <- 6988
num_noninvar_snps <- total_snps - num_invar
num_noninvar_snps
## [1] 5267

Remove NA cols

NOTE: The for loop code takes about 5ish minutes to run, I recommend against running it. I did it once and found this data has 0 NAs in any SNP col.

find_NAs <- function(x){
  NAs_TF <- is.na(x)
  i_NA <- which(NAs_TF == TRUE)
  N_NA <- length(i_NA)
  
  cat("Results:",N_NA, "NAs present\n.")
  return(i_NA)
}
#N_rows <- nrow(vcf_noinvar)

# N_NA
# vector to hold output (number of NAs)
#N_NA   <- rep(x = 0, times = N_rows)

# N_SNPs
# total number of columns (SNPs)
#N_SNPs <- ncol(vcf_noinvar)

# the for() loop
#for(i in 1:N_rows){
  
  # for each row, find the location of
  ## NAs with snps_num_t()
  #i_NA <- find_NAs(vcf_noinvar[i,]) 
  
  # then determine how many NAs
  ## with length()
  #N_NA_i <- length(i_NA)
  
  # then save the output to 
  ## our storage vector
  #N_NA[i] <- N_NA_i
#}
# No NAs, so mean imputaion not necessary
N_NAs <- 0

Scale data

vcf_scaled <- vcf_noinvar
vcf_scaled[,-c(1:6)] <- scale(vcf_noinvar[,-c(1:6)])
names
## function (x)  .Primitive("names")

Save cleaned data to .csv

write.csv(vcf_scaled, file = "ity1_clean_data.csv", row.names = F)
#Overview of cleaned data
cat("dimentions of clean data:", dim(vcf_scaled),"\nMetadata: ", names(vcf_scaled)[1:6],"\nSNPs:\n\tTotal SNPS:",total_snps,"\n\tInvariate SNPS:", num_invar,"\n\tSNPS in clean:", num_noninvar_snps,"\n\nSamples:",nrow(vcf_scaled))
## dimentions of clean data: 2504 5273 
## Metadata:  sample pop super_pop sex lat lng 
## SNPs:
##  Total SNPS: 6988 
##  Invariate SNPS: 1721 
##  SNPS in clean: 5267 
## 
## Samples: 2504