final project

LOAD DATA

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)

## Warning: package 'vegan' was built under R version 4.2.2

## Loading required package: permute

## Warning: package 'permute' was built under R version 4.2.2

## Loading required package: lattice

## This is vegan 2.6-4

library(ggplot2)
library(ggpubr)

# The VCF file was loaded into R (function `read.vcfR`)
vcf <- vcfR::read.vcfR("14.24353369-24593369.ALL.chr14_GRCh38.genotypes.20170504.vcf.gz", convertNA  = T)

## Scanning file to determine attributes.
## File attributes:
##   meta lines: 130
##   header_line: 131
##   variant count: 7710
##   column count: 2513
## 
Meta line 130 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
##   Character matrix gt rows: 7710
##   Character matrix gt cols: 2513
##   skip: 0
##   nrows: 7710
##   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: 7710
## All variants processed

# converted to counts of the minor allele
vcf_num <- vcfR::extract.gt(vcf, 
           element = "GT",
           IDtoRowNames = F,
           as.numeric = T,
           convertNA = T,
           return.alleles = F)

vcf_num_df <- data.frame(t(vcf_num))

# row names are added to the df as a new column
sample <- row.names(vcf_num_df)
vcf_num_df <- data.frame(sample, vcf_num_df)


# merge categorical with snp data
pop_meta <- read.csv(file = "1000genomes_people_info2-1.csv")
vcf_num_df <- merge(pop_meta, 
                 vcf_num_df,
                 by = 'sample')

# confirm successful merge and save data as '.csv'
write.csv(vcf_num_df, 
          file = "vcf_num_df.csv",
          row.names = F)

CLEAN DATA

# screen for any invariant SNPs using the `invar_omit()`
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)
}

vcf_no_invar <- vcf_num_df[,-c(1:6)]
vcf_no_invar <- invar_omit(vcf_no_invar)

## Dataframe of dim 2504 7710 processed...
## 1914 columns removed

vcf_no_invar <-  data.frame(vcf_num_df[,c("sample","pop","super_pop","sex","lat","lng")],
                            vcf_no_invar)

invars_removed <- 1914

# data is **centered and scaled** using R's `scale()` function.
vcf_scaled <- vcf_no_invar
vcf_scaled[,-c(1:6)] <- scale(vcf_no_invar[,-c(1:6)])

RUN PCA

vcf_pca <- prcomp(vcf_scaled[,-c(1:6)])
screeplot(vcf_pca)

PCA_variation <- function(pca_summary, PCs){
  var_explained <- pca_summary$importance[2,1:PCs]*100
  var_explained <- round(var_explained, 3)
  return(var_explained)
}

vcf_pca_summary <- summary(vcf_pca)

var_out <- PCA_variation(vcf_pca_summary, 
                         PCs = 100)

N_columns <- ncol(vcf_scaled)
cut_off <- 1/N_columns*100

i_cut_off <- which(var_out < cut_off)
i_cut_off <- min(i_cut_off)

## Warning in min(i_cut_off): no non-missing arguments to min; returning Inf

my_meta_N_meanNA_rowsPCs <- i_cut_off

my_meta_var_PC123 <- var_out[c(1,2,3)]

PLOT RESULTS

barplot(var_out,
        main = "Percent variation (%) Scree plot",
        ylab = "Percent variation (%) explained",
        names.arg = 1:length(var_out))

abline(h = cut_off, col = 2, lwd = 2)
abline(v = i_cut_off)

legend("topright",
       col = c(2,1),
       lty = c(1,1),
       legend = c("Vertical line: cutoff",
                  "Horizontal line: 1st value below cutoff"))

# Plotting Cumulative Percent Variation
cumulative_variation <- cumsum(var_out)
plot(cumulative_variation, type = "l")

# PLOT PCA
vcf_pca_scores <- vegan::scores(vcf_pca)
vcf_pca_scores2 <- data.frame(super_pop = vcf_no_invar$super_pop,
                              vcf_pca_scores)

# variation explained by the PCs
my_meta_var_PC123[1]

##   PC1 
## 4.684

my_meta_var_PC123[2]

##   PC2 
## 2.908

# PC1 vs PC2
ggpubr::ggscatter(data = vcf_pca_scores2,
                  x = "PC1",
                  y = "PC2",
                  color = "super_pop",
                  shape = "super_pop",
                  main = "PCA Scatterplot",
                  xlab = paste("PC1 (",var_out[1],"% variation)", sep = ""),
                  ylab = paste("PC2 (",var_out[2],"% variation)", sep = ""))

# PC2 vs PC3
ggpubr::ggscatter(data = vcf_pca_scores2,
                  x = "PC2",
                  y = "PC3",
                  color = "super_pop",
                  shape = "super_pop",
                  main = "PCA Scatterplot",
                  xlab = paste("PC2 (",var_out[2],"% variation)", sep = ""),
                  ylab = paste("PC3 (",var_out[3],"% variation)", sep = ""))

# PC1 vs PC3
ggpubr::ggscatter(data = vcf_pca_scores2,
                  x = "PC1",
                  y = "PC3",
                  color = "super_pop",
                  shape = "super_pop",
                  main = "PCA Scatterplot",
                  xlab = paste("PC1 (",var_out[1],"% variation)", sep = ""),
                  ylab = paste("PC3 (",var_out[3],"% variation)", sep = ""))