Loading Necessary Packages

library(vcfR)
## Warning: package 'vcfR' was built under R version 4.1.2
## 
##    *****       ***   vcfR   ***       *****
##    This is vcfR 1.13.0 
##      browseVignettes('vcfR') # Documentation
##      citation('vcfR') # Citation
##    *****       *****      *****       *****
library(vegan)
## Warning: package 'vegan' was built under R version 4.1.2
## Loading required package: permute
## Warning: package 'permute' was built under R version 4.1.2
## Loading required package: lattice
## This is vegan 2.6-4

Loading Data Set

my_vcf1 <- "5.24000-264000.ALL.chr5_GRCh38.genotypes.20170504.vcf.gz"

Chrom20 <- vcfR::read.vcfR(my_vcf1, convertNA = T)
## Scanning file to determine attributes.
## File attributes:
##   meta lines: 130
##   header_line: 131
##   variant count: 8412
##   column count: 2513
## 
Meta line 130 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
##   Character matrix gt rows: 8412
##   Character matrix gt cols: 2513
##   skip: 0
##   nrows: 8412
##   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: 8412
## All variants processed

Extract genotype scores, transpose, and create data frame

snp_gt <- vcfR::extract.gt(Chrom20, # TODO
           element = "GT",
           IDtoRowNames  = F,
           as.numeric = T,
           convertNA = T)
#Transpose VCF orientation into R dataframe orientation
snp_gt_t <- t(snp_gt)

#Create dataframe 
snp_gt_df <- data.frame(snp_gt_t)

#Get person (sample) name
sample <- row.names(snp_gt_df)

#Add sample infor into dataframe
snp_sample_df <- data.frame(sample, snp_gt_df)

Data Clean-Up

Load population Data and merge existing data frame

pop_meta <- read.csv("1000genomes_people_info2.csv")

names(pop_meta)
## [1] "pop"       "super_pop" "sample"    "sex"       "lat"       "lng"
names(snp_sample_df)[1:10]
##  [1] "sample" "X1"     "X2"     "X3"     "X4"     "X5"     "X6"     "X7"    
##  [9] "X8"     "X9"
vcf_df2 <- merge(pop_meta, snp_sample_df, by = "sample")

Check dimensions after the merge

nrow(vcf_df2) == nrow(snp_sample_df)
## [1] TRUE

create a function for invariate data

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 names of df
names(vcf_df2)[1:10]
##  [1] "sample"    "pop"       "super_pop" "sex"       "lat"       "lng"      
##  [7] "X1"        "X2"        "X3"        "X4"
vcf_noinvar <- vcf_df2

vcf_noinvar[, -c(1:6)] <- invar_omit(vcf_noinvar[, -c(1:6)])
## Dataframe of dim 2504 8412 processed...
## 2037 columns removed
my_meta_N_invar_cols <- 2037

Remove 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

N_rows <- nrow(vcf_noinvar)
N_NA <- rep(x = 0, times = N_rows)
N_SNPs <- ncol(vcf_noinvar)

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 greater than 50% NAs

cutoff50 <- N_SNPs*0.5
percent_NA <- N_NA/N_SNPs*100
any(percent_NA > 50)
## [1] FALSE

Mean imputation

mean_imputation <- function(df){
  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)
}

Run Mean-imputation only on numeric columns (not 1-6)

vcf_noNA <- vcf_noinvar

vcf_noNA[, -c(1:6)] <- mean_imputation(vcf_noinvar[, -c(1:6)])

Prepare for PCA

#New Copy
vcf_scaled <- vcf_noNA

#Scale the data 
vcf_scaled[,-c(1:6)] <- scale(vcf_noNA[, -c(1:6)])

Write the Data

write.csv(vcf_scaled, file = "vcf_scaled.csv", row.names = F)