Load necessary R packages

Load and/or install all packages required for analysis.

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-2
library(ggplot2)
library(ggpubr)

Confirm location of working directory and presence of desired vcf file.

getwd()
## [1] "C:/Users/willi/Desktop/RStudio"
list.files(pattern = "vcf")
## [1] "1.159051856-159301856.ALL.chr1_GRCh38.genotypes.20170504.vcf.gz"
## [2] "12.12000-252000.ALL.chr12_GRCh38.genotypes.20170504.vcf.gz"     
## [3] "all_loci-1.vcf"                                                 
## [4] "all_loci.vcf"                                                   
## [5] "code_checkpoint_vcfR.Rmd"                                       
## [6] "vcf_for_PCA.csv"                                                
## [7] "vcf_num.csv"                                                    
## [8] "vcf_num_df.csv"                                                 
## [9] "vcf_num_df2.csv"
my_vcf <- "12.12000-252000.ALL.chr12_GRCh38.genotypes.20170504.vcf.gz"

Data preparation

Load 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: 8084
##   column count: 2513
## 
Meta line 130 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
##   Character matrix gt rows: 8084
##   Character matrix gt cols: 2513
##   skip: 0
##   nrows: 8084
##   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: 8084
## All variants processed

Convert raw VCF file to genotype scores and save as a csv

vcf_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)

Transpose VCF, create dataframe, and add row names information into the dataframe created

vcf_num_t <- t(vcf_num)

vcf_num_df <- data.frame(vcf_num_t)

sample<- row.names(vcf_num_df)

vcf_num_df <- data.frame(sample, 
                         vcf_num_df)

Save genotypic scores as a csv file and confirm its presence in working directory.

write.csv(vcf_num_df,
          file = "vcf_num_df.csv",
          row.names = F)
list.files(pattern = "csv")
## [1] "1000genomes_people_info2-1.csv" "gwas_pheno_env.csv"            
## [3] "SNPs_cleaned.csv"               "vcf_for_PCA.csv"               
## [5] "vcf_num.csv"                    "vcf_num_df.csv"                
## [7] "vcf_num_df2.csv"                "walsh2017morphology.csv"

Data subsetting and meta data

Load population meta data, and check that “sample” appears in the column names for both meta AND SNP 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"

Merge the two data sets and check the dimensions to see if the row numbers are the same before and after merging

vcf_num_df2 <- merge(pop_meta, vcf_num_df, by= "sample")
nrow(vcf_num_df) == nrow(vcf_num_df2)
## [1] TRUE

Check the names of the new dataframe created

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"

Make a csv file for vcf_num_df2

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

Omitting invariants 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)
}

Run invar_omit() on non-character columns and store to a new dataframe object.

vcf_noinvar<- vcf_num_df2
vcf_noinvar[, -c(1:6)] <- invar_omit(vcf_noinvar[, -c(1:6)])
## Dataframe of dim 2504 8084 processed... 
## 2007 columns removed
N_invar_cols <- 2007

Removing data with too many NAs Create function to find NAs

find_NAs <- function(x){
  NAs_TF <- is.na(x)
  i_NA <- which(NAs_TF ==T)
  N_NA <- length(i_NA)
  
  return(i_NA)
}

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

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

Mean imputation of NAs

NAs are rare. This is doing an imputation just in case.

mean_imputation <- function(df){
  n_cols <- ncol(df)
  
  for(i in 1:n_cols){
    column_i <- df[, i]
    mean_i <- mean(column_i, na.rm = TRUE)
    NAs_i <- which(is.na(column_i))
    N_NAs <- length(NAs_i)
    column_i[NAs_i] <- mean_i
    df[, i] <- column_i  }
  
  return(df)  }

Apply mean imputation to non-character columns.

names(vcf_noinvar)[1:10] 
##  [1] "sample"    "pop"       "super_pop" "sex"       "lat"       "lng"      
##  [7] "X1"        "X2"        "X3"        "X4"
vcf_noNA <- vcf_noinvar
vcf_noNA[,-c(1:6)] <- mean_imputation(vcf_noinvar[,-c(1:6)])

Preparing for PCA Scaling data

Prepare data by centering mean on 0 and scaling by standard deviation.

#new copy of data
vcf_scaled <- vcf_noNA

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

Write a csv file for PCA analysis

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