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:
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-4library(ggplot2)
library(ggpubr)
getwd()## [1] "/Users/ziqili/Desktop/COMPBIO Project"list.files(pattern = "vcf")## [1] "4.13132463-13372463.ALL.chr4_GRCh38.genotypes.20170504.vcf.gz"
## [2] "ALL.chr4_GRCh38.genotypes.20170504.vcf.gz"
## [3] "vcf_num_df.csv"
## [4] "vcf_num_df2.csv"my_vcf <- "4.13132463-13372463.ALL.chr4_GRCh38.genotypes.20170504.vcf.gz"
vcf <- vcfR::read.vcfR(my_vcf, convertNA = T)## Scanning file to determine attributes.
## File attributes:
## meta lines: 130
## header_line: 131
## variant count: 7497
## column count: 2513
##
Meta line 130 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
## Character matrix gt rows: 7497
## Character matrix gt cols: 2513
## skip: 0
## nrows: 7497
## 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: 7497
## All variants processedvcf_num <- vcfR::extract.gt(vcf, element = "GT", IDtoRowNames = F, as.numeric = T, convertNA = T)
list.files() ## [1] "1000genomes_people_info2-1.csv"
## [2] "4.13132463-13372463.ALL.chr4_GRCh38.genotypes.20170504.vcf.gz"
## [3] "ALL.chr4_GRCh38.genotypes.20170504.vcf.gz"
## [4] "COMPBIO FInal Project Workflow.Rmd"
## [5] "COMPBIO FInal Project.Rmd"
## [6] "Compbio Final_report.Rmd"
## [7] "COMPBIO Project.Rproj"
## [8] "Compbio-Final_report.docx"
## [9] "Compbio-Final_report.html"
## [10] "COMPBIO-FInal-Project-Workflow.Rmd"
## [11] "datascaled.csv"
## [12] "rsconnect"
## [13] "vcf_num_df.csv"
## [14] "vcf_num_df2.csv"vcf_num_t <- t(vcf_num)
vcf_num_df <- data.frame(vcf_num_t)
sample <- row.names(vcf_num_df)
vcf_num_df <- data.frame(sample, vcf_num_df)getwd() ## [1] "/Users/ziqili/Desktop/COMPBIO Project"write.csv(vcf_num_df, file = "vcf_num_df.csv", row.names = F)
list.files()## [1] "1000genomes_people_info2-1.csv"
## [2] "4.13132463-13372463.ALL.chr4_GRCh38.genotypes.20170504.vcf.gz"
## [3] "ALL.chr4_GRCh38.genotypes.20170504.vcf.gz"
## [4] "COMPBIO FInal Project Workflow.Rmd"
## [5] "COMPBIO FInal Project.Rmd"
## [6] "Compbio Final_report.Rmd"
## [7] "COMPBIO Project.Rproj"
## [8] "Compbio-Final_report.docx"
## [9] "Compbio-Final_report.html"
## [10] "COMPBIO-FInal-Project-Workflow.Rmd"
## [11] "datascaled.csv"
## [12] "rsconnect"
## [13] "vcf_num_df.csv"
## [14] "vcf_num_df2.csv"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] "sample" "X1" "X2" "X3" "X4" "X5" "X6" "X7"
## [9] "X8" "X9"vcf_num_df2 <- merge(pop_meta, vcf_num_df, by = "sample")nrow(vcf_num_df) == nrow(vcf_num_df2)## [1] TRUEnames(vcf_num_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/ziqili/Desktop/COMPBIO Project"write.csv(vcf_num_df2, file = "vcf_num_df2.csv", row.names = F)
list.files()## [1] "1000genomes_people_info2-1.csv"
## [2] "4.13132463-13372463.ALL.chr4_GRCh38.genotypes.20170504.vcf.gz"
## [3] "ALL.chr4_GRCh38.genotypes.20170504.vcf.gz"
## [4] "COMPBIO FInal Project Workflow.Rmd"
## [5] "COMPBIO FInal Project.Rmd"
## [6] "Compbio Final_report.Rmd"
## [7] "COMPBIO Project.Rproj"
## [8] "Compbio-Final_report.docx"
## [9] "Compbio-Final_report.html"
## [10] "COMPBIO-FInal-Project-Workflow.Rmd"
## [11] "datascaled.csv"
## [12] "rsconnect"
## [13] "vcf_num_df.csv"
## [14] "vcf_num_df2.csv"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)
}names(vcf_num_df2)[1:10]## [1] "sample" "pop" "super_pop" "sex" "lat" "lng"
## [7] "X1" "X2" "X3" "X4"vcf_noinvar <- vcf_num_df2
vcf_noinvar <- invar_omit(vcf_noinvar[, -c(1:6)]) ## Dataframe of dim 2504 7497 processed...
## 1877 columns removedvcf_noinvar <- invar_omit(vcf_noinvar)## Dataframe of dim 2504 5620 processed...
## 0 columns removedvcf_noinvar <- data.frame(vcf_num_df2[, c(1:6)], vcf_noinvar)
dim(vcf_noinvar) ## [1] 2504 5626my_meta_N_invar_cols <- 1877find_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_noinvar)
N_NA <- rep(x=0, times = N_rows)
N_SNPs <- ncol(vcf_noinvar)
cat("This may take a minutes...")## This may take a minutes...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
}cutoff50 <- N_SNPs*0.5
percent_NA <- N_NA/N_SNPs*100
any(percent_NA > 50)## [1] FALSEmean(percent_NA)## [1] 0my_meta_N_meanNA_rows <- mean(percent_NA)mean_imputation <- function(df){
cat("This make take some time...")
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)])## This make take some time...vcf_scaled <- vcf_noNA
vcf_scaled[, -c(1:6)] <- scale(vcf_noNA[, -c(1:6)])
vcf_pca <- prcomp(vcf_scaled[,-c(1:6)])write.csv(vcf_scaled, file = "datascaled.csv", row.names = F)