Introduction

In this worked example you will replicate a PCA on a published dataset.

The example is split into 2 Parts:

In this Data Preparation phase, you will do the following things:

  1. Load the SNP genotypes in .vcf format (vcfR::read.vcfR())
  2. Extract the genotypes into an R-compatible format (vcfR::extract.gt())
  3. Rotate the data into the standard R analysis format (t())
  4. Remove individuals (rows) from the data set that have >50% NAs (using a function I wrote)
  5. Remove SNPs (columns) that are fixed
  6. Impute remaining NAs (using a for() loop)
  7. Save the prepared data as a .csv file for the next step (write.csv())

Biological background

This worked example is based on a paper in the journal Molecular Ecology from 2017 by Jennifer Walsh titled Subspecies delineation amid phenotypic, geographic and genetic discordance in a songbird.

The study investigated variation between two bird species in the genus Ammodramus: A. nenlsoni and A. caudacutus.

The species A. nenlsoni has been divided into 3 sub-species: A. n. nenlsoni, A.n. alterus, and A n. subvirgatus. The other species, A. caudacutus, has been divided into two subspecies, A.c. caudacutus and A.c. diversus.

The purpose of this study was to investigate to what extent these five subspecies recognized by taxonomists are supported by genetic data. The author’s collected DNA from 75 birds (15 per subspecies) and genotyped 1929 SNPs. They then analyzed the data with Principal Components Analysis (PCA), among other genetic analyzes.

This tutorial will work through all of the steps necessary to re-analyze Walsh et al.s data

Tasks

In the code below all code is provided. Your tasks will be to do 2 things:

  1. Give a meaningful title to all sections marked “TODO: TITLE”
  2. Write 1 to 2 sentences describing what is being done and why in all sections marked “TODO: EXPLAIN”

Preliminaries

Load the vcfR and other packages with library().

library(vcfR)    
## Warning: package 'vcfR' was built under R version 4.2.2
## 
##    *****       ***   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)

Make sure that your working directory is set to the location of the file all_loci.vcf.

getwd()
## [1] "C:/Users/afhar/Desktop/Comp Bio"
list.files()
##  [1] "07-mean_imputation.docx"                        
##  [2] "07-mean_imputation.html"                        
##  [3] "07-mean_imputation.Rmd"                         
##  [4] "08-PCA_worked.html"                             
##  [5] "08-PCA_worked.Rmd"                              
##  [6] "09-PCA_worked_example-SNPs-part1.Rmd"           
##  [7] "all_loci-1.vcf"                                 
##  [8] "all_loci.vcf"                                   
##  [9] "bird_snps_remove_NAs.html"                      
## [10] "bird_snps_remove_NAs.Rmd"                       
## [11] "blue snp.png"                                   
## [12] "center_function (1).R"                          
## [13] "center_function.R"                              
## [14] "code_checkpoint_vcfR.Rmd"                       
## [15] "Comp Bio.Rproj"                                 
## [16] "data slicer preview results.png"                
## [17] "feature_engineering.Rmd"                        
## [18] "feature_engineering_intro_2_functions-part2.Rmd"
## [19] "graphs"                                         
## [20] "html"                                           
## [21] "in class 10-13.R"                               
## [22] "install screenshots"                            
## [23] "misc"                                           
## [24] "practice"                                       
## [25] "R 4.2.1.lnk"                                    
## [26] "removing_fixed_alleles.html"                    
## [27] "removing_fixed_alleles.Rmd"                     
## [28] "rsconnect"                                      
## [29] "SNPs_cleaned.csv"                               
## [30] "temp.Rmd"                                       
## [31] "test.Rmd"                                       
## [32] "transpose_VCF_data.html"                        
## [33] "transpose_VCF_data.Rmd"                         
## [34] "vcfR_test.vcf"                                  
## [35] "vcfR_test.vcf.gz"                               
## [36] "walsh2017morphology.csv"                        
## [37] "word docs"                                      
## [38] "working_directory_practice.html"                
## [39] "working_directory_practice.Rmd"
list.files(pattern = "vcf")
## [1] "all_loci-1.vcf"           "all_loci.vcf"            
## [3] "code_checkpoint_vcfR.Rmd" "vcfR_test.vcf"           
## [5] "vcfR_test.vcf.gz"

Data preparation

LOAD DATA

VCF File must be loaded and read in order to proceed with analysis of the data in said file. To do this must use vcfR, which was a previously loaded package.

snps <- vcfR::read.vcfR("all_loci.vcf", convertNA  = TRUE)
## Scanning file to determine attributes.
## File attributes:
##   meta lines: 8
##   header_line: 9
##   variant count: 1929
##   column count: 81
## 
Meta line 8 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
##   Character matrix gt rows: 1929
##   Character matrix gt cols: 81
##   skip: 0
##   nrows: 1929
##   row_num: 0
## 
Processed variant 1000
Processed variant: 1929
## All variants processed

EXTRACT DATA

After data is loaded it must be extracted in a way it can be utilized and transformed in the further analysis. The genotype of the SNPs in question is what is being extracted and analysed in subsequent steps.

snps_num <- vcfR::extract.gt(snps, 
           element = "GT",
           IDtoRowNames  = F,
           as.numeric = T,
           convertNA = T,
           return.alleles = F)

TRANSPOSE DATA

The extracted data must be transposed, flipped, so the samples become rows instead of columns for the rest of analysis.

snps_num_t <- t(snps_num)

CREATE DATAFRAME

The transposed data must be turned into a dataframe to proceed with further analysis.

snps_num_df <- data.frame(snps_num_t)

FIND LOCATION OF NAs IN DATASET

Create function that finds NA values within the data and returns how many NAs are present from a sample. If too many NAs occur from a single sample it may be better to omit the whole sample rather than impute/omit every singular value.

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

Creating a for loop to go through every sample in the data decreases the amount of code written to produce the same result. This cycles through every sample in the data and returns the number of NAs present for each. This number will help decide which samples to omit from analysis (if any)

# N_rows
# number of rows (individuals)
N_rows <- nrow(snps_num_t)

# 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(snps_num_t)

# 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(snps_num_t[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
}
## Results: 28 NAs present
## Results: 20 NAs present
## Results: 28 NAs present
## Results: 24 NAs present
## Results: 23 NAs present
## Results: 63 NAs present
## Results: 51 NAs present
## Results: 38 NAs present
## Results: 34 NAs present
## Results: 24 NAs present
## Results: 48 NAs present
## Results: 21 NAs present
## Results: 42 NAs present
## Results: 78 NAs present
## Results: 45 NAs present
## Results: 21 NAs present
## Results: 42 NAs present
## Results: 34 NAs present
## Results: 66 NAs present
## Results: 54 NAs present
## Results: 59 NAs present
## Results: 52 NAs present
## Results: 47 NAs present
## Results: 31 NAs present
## Results: 63 NAs present
## Results: 40 NAs present
## Results: 40 NAs present
## Results: 22 NAs present
## Results: 60 NAs present
## Results: 48 NAs present
## Results: 961 NAs present
## Results: 478 NAs present
## Results: 59 NAs present
## Results: 26 NAs present
## Results: 285 NAs present
## Results: 409 NAs present
## Results: 1140 NAs present
## Results: 600 NAs present
## Results: 1905 NAs present
## Results: 25 NAs present
## Results: 1247 NAs present
## Results: 23 NAs present
## Results: 750 NAs present
## Results: 179 NAs present
## Results: 433 NAs present
## Results: 123 NAs present
## Results: 65 NAs present
## Results: 49 NAs present
## Results: 192 NAs present
## Results: 433 NAs present
## Results: 66 NAs present
## Results: 597 NAs present
## Results: 1891 NAs present
## Results: 207 NAs present
## Results: 41 NAs present
## Results: 268 NAs present
## Results: 43 NAs present
## Results: 110 NAs present
## Results: 130 NAs present
## Results: 90 NAs present
## Results: 271 NAs present
## Results: 92 NAs present
## Results: 103 NAs present
## Results: 175 NAs present
## Results: 31 NAs present
## Results: 66 NAs present
## Results: 64 NAs present
## Results: 400 NAs present
## Results: 192 NAs present
## Results: 251 NAs present
## Results: 69 NAs present
## Results: 58 NAs present

The decided cut off for omition of a sample is if 50% of the SNPs are an NA value. This takes the information of each sample and shows how many samples fall above/below this cut off line.

# 50% of N_SNPs
cutoff50 <- N_SNPs*0.5

hist(N_NA)            
abline(v = cutoff50, 
       col = 2, 
       lwd = 2, 
       lty = 2)

To find the percentage of NAs present the number of NAs in the sample must be divided by the total number of SNPs in said sample, then multiplied by 100. This percentage will then be compared to the cut off of 50% to decide whether the sample is used in further analysis or omitted.

percent_NA <- N_NA/N_SNPs*100

# Call which() on percent_NA
i_NA_50percent <- which(percent_NA > 50) 

snps_num_t02 <- snps_num_t[-i_NA_50percent, ]

ALTERING ROW NAMES

Regular expressions are used to edit the text in the row names, to provide clearer information about the samples we are focused on in this step.

row_names <- row.names(snps_num_t02) # Key

row_names02 <- gsub("sample_","",row_names)

sample_id <- gsub("^([ATCG]*)(_)(.*)",
                  "\\3",
                  row_names02)
pop_id <- gsub("[01-9]*",    
               "",
               sample_id)

table(pop_id)  
## pop_id
## Alt Cau Div Nel Sub 
##  15  12  15  15  11

REMOVING INVARIANT DATA

Some SNPs may be invariant, meaning it is a waste of computational power to continue analysis with that data present, as it takes up space and time but provides no new means of analysis, so it is better to remove it.

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


snps_no_invar <- invar_omit(snps_num_t02) 
## Dataframe of dim 68 1929 processed...
## 591 columns removed

MEAN IMPUTATION

After removing samples that have too many NAs and invariant data, there is still likely NA values present in the data deemed sufficient for analysis. These values are often treated with mean imputation, where the mean of the data is put in place of the NA values.

snps_noNAs <- snps_no_invar

N_col <- ncol(snps_no_invar)
for(i in 1:N_col){
  
  # get the current column
  column_i <- snps_noNAs[, 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))
  
  # record 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
  snps_noNAs[, i] <- column_i
  
}

Save the data

Save the data as a .csv file which can be loaded again later.

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

Check for the presence of the file with list.files()

list.files(pattern = ".csv")
## [1] "SNPs_cleaned.csv"        "walsh2017morphology.csv"

Next steps:

In Part 2, we will re-load the SNPs_cleaned.csv file and carry an an analysis with PCA.