R Markdown

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this: ## Load Necessary Packages

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(ggpubr)
## Loading required package: ggplot2
library(ggplot2)

Confirm Working Directory and File Location

getwd()
## [1] "/Users/sumedhasripada/Downloads/Final Project CB"
list.files(pattern = "vcf")
## [1] "sumedha_snps.vcf.gz" "vcf_num_df.csv"      "vcf_num_df2.csv"    
## [4] "vcf_num.csv"

Load VCF Data

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

Convert Raw VCF File to Genotyope Scores

vcf_num <- vcfR::extract.gt(vcf,
                            element = "GT",
                            IDtoRowNames = F,
                            as.numeric = T,
                            convertNA = T)

Save the CSV

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

Confirm Presence of File

list.files()
##  [1] "1000genomes_people_info2-1.csv" "Final Project CB.Rproj"        
##  [3] "final_report_template.Rmd"      "rsconnect"                     
##  [5] "scaled_data.csv"                "sumedha_snps.vcf.gz"           
##  [7] "vcf_num_df.csv"                 "vcf_num_df2.csv"               
##  [9] "vcf_num.csv"                    "Workflow for PCA Analysis.Rmd" 
## [11] "Workflow-for-PCA-Analysis.html" "Workflow-for-PCA-Analysis.Rmd"

Transpose Original VCF Orientation to R Dataframe Orientation

vcf_num_t <- t(vcf_num)

Make into a Data Frame

vcf_num_df <- data.frame(vcf_num_t)

sample <- row.names(vcf_num_df)

vcf_num_df <- data.frame(vcf_num_df, 
                         sample)
getwd()
## [1] "/Users/sumedhasripada/Downloads/Final Project CB"
write.csv(vcf_num_df,
          file = "vcf_num_df.csv",
          row.names = F)

Clean 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] "X1"  "X2"  "X3"  "X4"  "X5"  "X6"  "X7"  "X8"  "X9"  "X10"

Merge Datasets

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

Check Names of the new DF

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

getwd()
## [1] "/Users/sumedhasripada/Downloads/Final Project CB"

Save CSV

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

Confirm Presence of File

list.files()
##  [1] "1000genomes_people_info2-1.csv" "Final Project CB.Rproj"        
##  [3] "final_report_template.Rmd"      "rsconnect"                     
##  [5] "scaled_data.csv"                "sumedha_snps.vcf.gz"           
##  [7] "vcf_num_df.csv"                 "vcf_num_df2.csv"               
##  [9] "vcf_num.csv"                    "Workflow for PCA Analysis.Rmd" 
## [11] "Workflow-for-PCA-Analysis.html" "Workflow-for-PCA-Analysis.Rmd"

Omit Invariant Features

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 removes\n")
  if(length(i_var0)>0){
    x <- x[, -i_var0]
  }
  
  return(x)
}
names(vcf_num_df2)[1:10]
##  [1] "sample"    "pop"       "super_pop" "sex"       "lat"       "lng"      
##  [7] "X1"        "X2"        "X3"        "X4"

Create New Data Frame

vcf_noinvar <- vcf_num_df2
vcf_noinvar[, -c(1:6)] <- invar_omit(vcf_noinvar[, -c(1:6)])
## Dataframe of dim 2504 7227 processed...
## 1805 columns removes

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

cutoff50 <- N_SNPs*.5
percent_NA <- N_NA/N_SNPs*100
any(percent_NA > 50)
## [1] FALSE
mean(percent_NA)
## [1] 0

Imputation of NAs

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

Scale Data to Prepare for PCA

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

Write CSV File

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