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/wyattkriebel/Desktop/data sclicer/Final Project"
list.files(pattern="vcf")
## [1] "13.vcf.gz"       "vcf_num_df.csv"  "vcf_num_df2.csv" "vcf_num.csv"
chromosome_13 <- "13.vcf.gz"
  
vcf <- vcfR::read.vcfR(chromosome_13, convertNA = T)
## Scanning file to determine attributes.
## File attributes:
##   meta lines: 130
##   header_line: 131
##   variant count: 6932
##   column count: 2513
## 
Meta line 130 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
##   Character matrix gt rows: 6932
##   Character matrix gt cols: 2513
##   skip: 0
##   nrows: 6932
##   row_num: 0
## 
Processed variant 1000
Processed variant 2000
Processed variant 3000
Processed variant 4000
Processed variant 5000
Processed variant 6000
Processed variant: 6932
## All variants processed
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)

list.files(pattern = "csv")
## [1] "1000genomes_people_info2.csv" "Data_cleaned.csv"            
## [3] "vcf_num_df.csv"               "vcf_num_df2.csv"             
## [5] "vcf_num.csv"

##pop_meta <- read.csv(file = “1000 genomes people info”)

vcf_num_transpose <- t(vcf_num)

vcf_df <- data.frame(vcf_num_transpose)

sample <- row.names(vcf_df)

vcf_df <- data.frame(sample, vcf_df)
getwd()
## [1] "/Users/wyattkriebel/Desktop/data sclicer/Final Project"
list.files()
##  [1] "1000genomes_people_info2.csv" "13.vcf.gz"                   
##  [3] "Data_cleaned.csv"             "final project.Rmd"           
##  [5] "Final Project.Rproj"          "final_report_template.Rmd"   
##  [7] "final-project.html"           "final-project.Rmd"           
##  [9] "vcf_num_df.csv"               "vcf_num_df2.csv"             
## [11] "vcf_num.csv"
write.csv(vcf_df, file = "vcf_num_df.csv", row.names = F )
list.files(pattern="csv")
## [1] "1000genomes_people_info2.csv" "Data_cleaned.csv"            
## [3] "vcf_num_df.csv"               "vcf_num_df2.csv"             
## [5] "vcf_num.csv"

##Time to clean da data

list.files()
##  [1] "1000genomes_people_info2.csv" "13.vcf.gz"                   
##  [3] "Data_cleaned.csv"             "final project.Rmd"           
##  [5] "Final Project.Rproj"          "final_report_template.Rmd"   
##  [7] "final-project.html"           "final-project.Rmd"           
##  [9] "vcf_num_df.csv"               "vcf_num_df2.csv"             
## [11] "vcf_num.csv"
pop_meta <- read.csv(file = "1000genomes_people_info2.csv")
names(pop_meta)
## [1] "pop"       "super_pop" "sample"    "sex"       "lat"       "lng"
names(vcf_df)[1:10]
##  [1] "sample" "X1"     "X2"     "X3"     "X4"     "X5"     "X6"     "X7"    
##  [9] "X8"     "X9"
vcf_df2 <- merge(pop_meta, vcf_df, by = "sample")

nrow(vcf_df) == nrow(vcf_df2)
## [1] TRUE
names(vcf_df2)[1:15]
##  [1] "sample"    "pop"       "super_pop" "sex"       "lat"       "lng"      
##  [7] "X1"        "X2"        "X3"        "X4"        "X5"        "X6"       
## [13] "X7"        "X8"        "X9"
getwd()
## [1] "/Users/wyattkriebel/Desktop/data sclicer/Final Project"
write.csv(vcf_df2, file = "vcf_num_df2.csv", row.names = F)
list.files()
##  [1] "1000genomes_people_info2.csv" "13.vcf.gz"                   
##  [3] "Data_cleaned.csv"             "final project.Rmd"           
##  [5] "Final Project.Rproj"          "final_report_template.Rmd"   
##  [7] "final-project.html"           "final-project.Rmd"           
##  [9] "vcf_num_df.csv"               "vcf_num_df2.csv"             
## [11] "vcf_num.csv"

##Omit invar function

invar_omit <- function(x){
  cat("Datafram dim", dim(x), "processed...\n")
  sds <- apply(x, 2, sd, na.rm = TRUE)
  i_var0 <- which(sds ==0)
  
  cat(length(i_var0), "Columns remove\n")
  
  if(length(i_var0) > 0){
    x <- x[,-i_var0]
  }
  
  return(x)
}
names(vcf_df2)[1:10]
##  [1] "sample"    "pop"       "super_pop" "sex"       "lat"       "lng"      
##  [7] "X1"        "X2"        "X3"        "X4"

#NOTES: hyperparameter –> the user puts in the data that they want

vcf_no_invariants <- vcf_df2
vcf_no_invariants[,-c(1:6)] <- invar_omit(vcf_no_invariants[,-c(1:6)])
## Datafram dim 2504 6932 processed...
## 1699 Columns remove
my_meta_N_invar_columns <- 1699
find_NAs <- function(x){
  NAs_TF <- is.na(x)
  i_NA <- which(NAs_TF == TRUE)
  N_NA <- length(i_NA)
  
  return(i_NA)
}
N_rows <- nrow(vcf_no_invariants)

N_NA <- rep(x=0, times = N_rows)

N_SNPs <- ncol(vcf_no_invariants)

cat("This may take a minute...")
## This may take a minute...
for(i in 1:N_rows){
  i_NA <- find_NAs(vcf_no_invariants[i,])
  N_NA_i <- length(i_NA)
  N_NA[i] <- N_NA_i
}
cutoff50percent <- N_SNPs*0.5
percent_NA <- N_NA/N_SNPs*100
any(percent_NA > 50)
## [1] FALSE
mean(percent_NA)
## [1] 0.0003799053
my_mean_N_meanNA_rows <- mean(percent_NA)
mean_imputation <- function(df){
  cat("this may take some time, go watch a youtube video...")
  
  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) == TRUE)
    N_NAs <- length(NAs_i)
    
    df[NAs_i,i] <- mean_i
    
  }
  return(df)
}
names(vcf_no_invariants)[1:10]
##  [1] "sample"    "pop"       "super_pop" "sex"       "lat"       "lng"      
##  [7] "X1"        "X2"        "X3"        "X4"
vcf_noNA <- vcf_no_invariants
vcf_noNA[,-c(1:6)] <- mean_imputation(vcf_no_invariants[,-c(1:6)])
## this may take some time, go watch a youtube video...
write.csv(vcf_noNA, file = "Data_cleaned.csv", row.names = F)