Final Project Workup (Pt1)

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

getwd()

## [1] "/Users/eeshamukherjee/Downloads/Computaional Bio 1/Final Project"

list.files(pattern="vcf")

## [1] "all_loci.vcf"  "My_snp.vcf.gz"

snps <- vcfR::read.vcfR("My_snp.vcf.gz", convertNA= T)

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

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

snps_num_t <- t(snps_num) 

snps_num_df <- data.frame(snps_num_t) 
sample <- row.names(snps_num_df)
snps_num_df <- data.frame(sample, snps_num_df)

pop_meta <- read.csv("1000genomes_people_info2-1.csv")
snps_num_df2 <- merge(pop_meta,
                      snps_num_df,
                      by = "sample",)

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)                      
}
snps_no_invar <- data.frame(snps_num_df2[,c(1:6)],
                            invar_omit(snps_num_df2[,-c(1:6)]),
                            row.names = NULL)

## Dataframe of dim 2504 7109 processed...
## 1844 columns removed

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

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

for (i in 1:N_rows) {
  i_NA <- find_NAs(snps_no_invar[i,])
  N_NA_i <- length(i_NA)
  N_NA[i] <- N_NA_i
}

cutoff50 <- N_SNPs*0.5
percent_NA <- N_NA/N_SNPs*100
any(percent_NA > 50)

## [1] FALSE

mean(percent_NA)

## [1] 0

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 = T)
    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(snps_no_invar)[1:10] 

##  [1] "sample"    "pop"       "super_pop" "sex"       "lat"       "lng"      
##  [7] "X2"        "X3"        "X7"        "X8"

snp_noNA <- snps_no_invar
snp_noNA[, -c(1:6)]<- mean_imputation(snps_no_invar[, -c(1:6)])

snp_scaled <- snp_noNA
snp_scaled[, -c(1:6)]<- scale(snp_noNA[, -c(1:6)])

getwd()

## [1] "/Users/eeshamukherjee/Downloads/Computaional Bio 1/Final Project"

write.csv(snp_scaled,
          file = "snp_scaled.csv",
          row.names = F)
list.files(pattern = "csv")

## [1] "1000genomes_people_info2-1.csv" "snp_scaled.csv"                
## [3] "SNPs_cleaned.csv"               "walsh2017morphology.csv"

```