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)Check working directory
getwd()## [1] "/Users/ekoneduncan/Desktop/compbio"list.files()## [1] "07-mean_imputation.html"
## [2] "07-mean_imputation.Rmd"
## [3] "08-PCA_worked.html"
## [4] "08-PCA_worked.Rmd"
## [5] "09-PCA_worked_example-SNPs-part1.html"
## [6] "09-PCA_worked_example-SNPs-part1.Rmd"
## [7] "1.159051856-159301856.ALL.chr1_GRCh38.genotypes.20170504.vcf.gz"
## [8] "10-PCA_worked_example-SNPs-part2.html"
## [9] "10-PCA_worked_example-SNPs-part2.Rmd"
## [10] "1000genomes_people_info2-1.csv"
## [11] "1540_120122-1.pdf"
## [12] "1540_final_project_Final_Report_template.pdf"
## [13] "1540_final_report_flowchart.pdf"
## [14] "1540_week14_PCA_SNP_workflow.pdf"
## [15] "17.12071392-12311392.ALL.chr17_GRCh38.genotypes.20170504.vcf"
## [16] "17.12071392-12311392.ALL.chr17_GRCh38.genotypes.20170504.vcf.gz"
## [17] "2.136483646-136733646.ALL.chr2_GRCh38.genotypes.20170504.vcf.gz"
## [18] "all_loci-1.vcf"
## [19] "all_loci.vcf"
## [20] "center_function.R"
## [21] "code_checkpoint_vcfR.html"
## [22] "code_checkpoint_vcfR.Rmd"
## [23] "data_prep.html"
## [24] "data_prep.Rmd"
## [25] "feature_engineering_intro_2_functions-part2.Rmd"
## [26] "feature_engineering.Rmd"
## [27] "final_report_template.Rmd"
## [28] "for_pca.csv"
## [29] "fst_exploration_in_class-STUDENT.html"
## [30] "fst_exploration_in_class-STUDENT.Rmd"
## [31] "fst_exploration_in_class.Rmd"
## [32] "gwas_pheno_env.csv"
## [33] "lecture-introd2RStudio-with_scripts.pdf"
## [34] "line_of_best_fit_example-tibet_allele_freq.pdf"
## [35] "my_snps"
## [36] "Navarro_regression_part01.pdf"
## [37] "PCA_analysis_in_class_work-for_students.pdf"
## [38] "PCA_with_SNPs_handout_worksheet_031122-1.docx"
## [39] "PCA_with_SNPs_handout_worksheet_031122-1.pdf"
## [40] "PCA-missing_data-KEY.Rmd"
## [41] "PCA-missing_data.Rmd"
## [42] "pheno.csv"
## [43] "R_data_structures_vectors_intro.pdf"
## [44] "R_Directory"
## [45] "r_help_hclust_intro-vs2.pdf"
## [46] "removing_fixed_alleles.html"
## [47] "removing_fixed_alleles.Rmd"
## [48] "rsconnect"
## [49] "SNPs_cleaned.csv"
## [50] "summary_stats.pdf"
## [51] "test.docx"
## [52] "test.html"
## [53] "test.Rmd"
## [54] "test2.Rmd"
## [55] "transpose_1000_genomes.html"
## [56] "transpose_1000_genomes.Rmd"
## [57] "vcf_data.csv"
## [58] "vcf_num_df.csv"
## [59] "vcf_num_df2.csv"
## [60] "vcf_num.csv"
## [61] "vcfR_test.vcf"
## [62] "vcfR_test.vcf.gz"
## [63] "vegan_PCA_amino_acids-STUDENT.html"
## [64] "vegan_PCA_amino_acids-STUDENT.Rmd"
## [65] "vegan_pca_with_msleep-STUDENT.html"
## [66] "vegan_pca_with_msleep-STUDENT.Rmd"
## [67] "walsh2017morphology.csv"
## [68] "walsh2017morphology.RData"
## [69] "week08_cluster_analysis-1.pdf"
## [70] "What is computational biology_exert.pdf"list.files(pattern = "vcf")## [1] "1.159051856-159301856.ALL.chr1_GRCh38.genotypes.20170504.vcf.gz"
## [2] "17.12071392-12311392.ALL.chr17_GRCh38.genotypes.20170504.vcf"
## [3] "17.12071392-12311392.ALL.chr17_GRCh38.genotypes.20170504.vcf.gz"
## [4] "2.136483646-136733646.ALL.chr2_GRCh38.genotypes.20170504.vcf.gz"
## [5] "all_loci-1.vcf"
## [6] "all_loci.vcf"
## [7] "code_checkpoint_vcfR.html"
## [8] "code_checkpoint_vcfR.Rmd"
## [9] "vcf_data.csv"
## [10] "vcf_num_df.csv"
## [11] "vcf_num_df2.csv"
## [12] "vcf_num.csv"
## [13] "vcfR_test.vcf"
## [14] "vcfR_test.vcf.gz"##Load VCF
Load vcf data
my_vcf<- "17.12071392-12311392.ALL.chr17_GRCh38.genotypes.20170504.vcf.gz"Load the vcf file
vcf<- vcfR::read.vcfR(my_vcf, convertNA = TRUE)## Scanning file to determine attributes.
## File attributes:
## meta lines: 130
## header_line: 131
## variant count: 7493
## column count: 2513
##
Meta line 130 read in.
## All meta lines processed.
## gt matrix initialized.
## Character matrix gt created.
## Character matrix gt rows: 7493
## Character matrix gt cols: 2513
## skip: 0
## nrows: 7493
## 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: 7493
## All variants processed##Convert raw VCF file to genotype scores
Get genotype score
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 the presence of file
list.files()## [1] "07-mean_imputation.html"
## [2] "07-mean_imputation.Rmd"
## [3] "08-PCA_worked.html"
## [4] "08-PCA_worked.Rmd"
## [5] "09-PCA_worked_example-SNPs-part1.html"
## [6] "09-PCA_worked_example-SNPs-part1.Rmd"
## [7] "1.159051856-159301856.ALL.chr1_GRCh38.genotypes.20170504.vcf.gz"
## [8] "10-PCA_worked_example-SNPs-part2.html"
## [9] "10-PCA_worked_example-SNPs-part2.Rmd"
## [10] "1000genomes_people_info2-1.csv"
## [11] "1540_120122-1.pdf"
## [12] "1540_final_project_Final_Report_template.pdf"
## [13] "1540_final_report_flowchart.pdf"
## [14] "1540_week14_PCA_SNP_workflow.pdf"
## [15] "17.12071392-12311392.ALL.chr17_GRCh38.genotypes.20170504.vcf"
## [16] "17.12071392-12311392.ALL.chr17_GRCh38.genotypes.20170504.vcf.gz"
## [17] "2.136483646-136733646.ALL.chr2_GRCh38.genotypes.20170504.vcf.gz"
## [18] "all_loci-1.vcf"
## [19] "all_loci.vcf"
## [20] "center_function.R"
## [21] "code_checkpoint_vcfR.html"
## [22] "code_checkpoint_vcfR.Rmd"
## [23] "data_prep.html"
## [24] "data_prep.Rmd"
## [25] "feature_engineering_intro_2_functions-part2.Rmd"
## [26] "feature_engineering.Rmd"
## [27] "final_report_template.Rmd"
## [28] "for_pca.csv"
## [29] "fst_exploration_in_class-STUDENT.html"
## [30] "fst_exploration_in_class-STUDENT.Rmd"
## [31] "fst_exploration_in_class.Rmd"
## [32] "gwas_pheno_env.csv"
## [33] "lecture-introd2RStudio-with_scripts.pdf"
## [34] "line_of_best_fit_example-tibet_allele_freq.pdf"
## [35] "my_snps"
## [36] "Navarro_regression_part01.pdf"
## [37] "PCA_analysis_in_class_work-for_students.pdf"
## [38] "PCA_with_SNPs_handout_worksheet_031122-1.docx"
## [39] "PCA_with_SNPs_handout_worksheet_031122-1.pdf"
## [40] "PCA-missing_data-KEY.Rmd"
## [41] "PCA-missing_data.Rmd"
## [42] "pheno.csv"
## [43] "R_data_structures_vectors_intro.pdf"
## [44] "R_Directory"
## [45] "r_help_hclust_intro-vs2.pdf"
## [46] "removing_fixed_alleles.html"
## [47] "removing_fixed_alleles.Rmd"
## [48] "rsconnect"
## [49] "SNPs_cleaned.csv"
## [50] "summary_stats.pdf"
## [51] "test.docx"
## [52] "test.html"
## [53] "test.Rmd"
## [54] "test2.Rmd"
## [55] "transpose_1000_genomes.html"
## [56] "transpose_1000_genomes.Rmd"
## [57] "vcf_data.csv"
## [58] "vcf_num_df.csv"
## [59] "vcf_num_df2.csv"
## [60] "vcf_num.csv"
## [61] "vcfR_test.vcf"
## [62] "vcfR_test.vcf.gz"
## [63] "vegan_PCA_amino_acids-STUDENT.html"
## [64] "vegan_PCA_amino_acids-STUDENT.Rmd"
## [65] "vegan_pca_with_msleep-STUDENT.html"
## [66] "vegan_pca_with_msleep-STUDENT.Rmd"
## [67] "walsh2017morphology.csv"
## [68] "walsh2017morphology.RData"
## [69] "week08_cluster_analysis-1.pdf"
## [70] "What is computational biology_exert.pdf"##Transpose data
Transpose original vcf orientation
vcf_num_t<- t(vcf_num)Make into dataframe
vcf_num_df<- data.frame(vcf_num_t)Get person(sample) names
sample <- row.names(vcf_num_df)Add sample info into dataframe
vcf_num_df<- data.frame(sample, vcf_num_df)save the csv
write.csv(vcf_num_df,
file = "vcf_num_df.csv",
row.names = F)list.files()## [1] "07-mean_imputation.html"
## [2] "07-mean_imputation.Rmd"
## [3] "08-PCA_worked.html"
## [4] "08-PCA_worked.Rmd"
## [5] "09-PCA_worked_example-SNPs-part1.html"
## [6] "09-PCA_worked_example-SNPs-part1.Rmd"
## [7] "1.159051856-159301856.ALL.chr1_GRCh38.genotypes.20170504.vcf.gz"
## [8] "10-PCA_worked_example-SNPs-part2.html"
## [9] "10-PCA_worked_example-SNPs-part2.Rmd"
## [10] "1000genomes_people_info2-1.csv"
## [11] "1540_120122-1.pdf"
## [12] "1540_final_project_Final_Report_template.pdf"
## [13] "1540_final_report_flowchart.pdf"
## [14] "1540_week14_PCA_SNP_workflow.pdf"
## [15] "17.12071392-12311392.ALL.chr17_GRCh38.genotypes.20170504.vcf"
## [16] "17.12071392-12311392.ALL.chr17_GRCh38.genotypes.20170504.vcf.gz"
## [17] "2.136483646-136733646.ALL.chr2_GRCh38.genotypes.20170504.vcf.gz"
## [18] "all_loci-1.vcf"
## [19] "all_loci.vcf"
## [20] "center_function.R"
## [21] "code_checkpoint_vcfR.html"
## [22] "code_checkpoint_vcfR.Rmd"
## [23] "data_prep.html"
## [24] "data_prep.Rmd"
## [25] "feature_engineering_intro_2_functions-part2.Rmd"
## [26] "feature_engineering.Rmd"
## [27] "final_report_template.Rmd"
## [28] "for_pca.csv"
## [29] "fst_exploration_in_class-STUDENT.html"
## [30] "fst_exploration_in_class-STUDENT.Rmd"
## [31] "fst_exploration_in_class.Rmd"
## [32] "gwas_pheno_env.csv"
## [33] "lecture-introd2RStudio-with_scripts.pdf"
## [34] "line_of_best_fit_example-tibet_allele_freq.pdf"
## [35] "my_snps"
## [36] "Navarro_regression_part01.pdf"
## [37] "PCA_analysis_in_class_work-for_students.pdf"
## [38] "PCA_with_SNPs_handout_worksheet_031122-1.docx"
## [39] "PCA_with_SNPs_handout_worksheet_031122-1.pdf"
## [40] "PCA-missing_data-KEY.Rmd"
## [41] "PCA-missing_data.Rmd"
## [42] "pheno.csv"
## [43] "R_data_structures_vectors_intro.pdf"
## [44] "R_Directory"
## [45] "r_help_hclust_intro-vs2.pdf"
## [46] "removing_fixed_alleles.html"
## [47] "removing_fixed_alleles.Rmd"
## [48] "rsconnect"
## [49] "SNPs_cleaned.csv"
## [50] "summary_stats.pdf"
## [51] "test.docx"
## [52] "test.html"
## [53] "test.Rmd"
## [54] "test2.Rmd"
## [55] "transpose_1000_genomes.html"
## [56] "transpose_1000_genomes.Rmd"
## [57] "vcf_data.csv"
## [58] "vcf_num_df.csv"
## [59] "vcf_num_df2.csv"
## [60] "vcf_num.csv"
## [61] "vcfR_test.vcf"
## [62] "vcfR_test.vcf.gz"
## [63] "vegan_PCA_amino_acids-STUDENT.html"
## [64] "vegan_PCA_amino_acids-STUDENT.Rmd"
## [65] "vegan_pca_with_msleep-STUDENT.html"
## [66] "vegan_pca_with_msleep-STUDENT.Rmd"
## [67] "walsh2017morphology.csv"
## [68] "walsh2017morphology.RData"
## [69] "week08_cluster_analysis-1.pdf"
## [70] "What is computational biology_exert.pdf"##Clean data
Load population metadata
meta_pop<- read.csv(file = "1000genomes_people_info2-1.csv")Merge data with SNP data
names(meta_pop)## [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(meta_pop,
vcf_num_df,
by = 'sample')Check the dimensions before and after merge
nrow(vcf_num_df) == nrow(vcf_num_df2)## [1] TRUEcheck names of new dataframe
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"save the csv
getwd()## [1] "/Users/ekoneduncan/Desktop/compbio"write.csv(vcf_num_df2, file = "vcf_num_df2.csv", row.names = F)
list.files(pattern = "csv")## [1] "1000genomes_people_info2-1.csv" "for_pca.csv"
## [3] "gwas_pheno_env.csv" "pheno.csv"
## [5] "SNPs_cleaned.csv" "vcf_data.csv"
## [7] "vcf_num_df.csv" "vcf_num_df2.csv"
## [9] "vcf_num.csv" "walsh2017morphology.csv"##Omit invariant features
Load invar_omit() function
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)
}Check which columns have character data
names(vcf_num_df2)[1:10]## [1] "sample" "pop" "super_pop" "sex" "lat" "lng"
## [7] "X1" "X2" "X3" "X4"##vcf_noinvar <- vcf_num_df2run invar_omit on numeric data and store results in new object
# the "patch"
# create a vector of just the numeric data, no columns
vcf_noinvar <- vcf_num_df2[, -c(1:7)]
# run invar omit on the data
vcf_noinvar <- invar_omit(vcf_noinvar)## Dataframe of dim 2504 7492 processed ...
## 1794 columns removed#put the metadata back with the numeric invar omitted data
vcf_noinvar <- data.frame(vcf_num_df2[, c("sample","pop","super_pop","sex","lat","lng")],
vcf_noinvar)
my_meta_N_invar_cols <- 1794##Remove low-quality data
Load find_NAs()
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
# number of rows (individuals)
N_rows <- nrow(vcf_noinvar)
#N_Na
# vector to hold output (number of NAs)
N_NA <- rep(x = 0, times = N_rows)
# N_SNPs
# total number of columns (SNPs)
N_SNPs <- ncol(vcf_noinvar)
cat("This may take a minute ... ")## This may take a minute ...## This may take a minute ...
# the for() Loop
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*0.5
percent_NA <- N_NA/N_SNPs*100
any(percent_NA > 50)## [1] FALSEmean(percent_NA)## [1] 0.0004410895my_meta_N_meanNA_rows <- mean(percent_NA)Load imputation function
mean_imputation <- function(df){
cat("This make take some time ... ")
n_cols <- ncol(df)
for(i in 1:n_cols){
# get the current column
column_i <- df[, i]
# get the mean of the current column
mean_i <- mean(column_i, na.rm = TRUE)
# get the NAs in the current column
NAs_i <- which(is.na(column_i))
# report the number of NAs
N_NAs <- length(NAs_i)
# replace the NAs in the current column
column_i[NAs_i] <- mean_i
# replace the original column with the
## updated columns
df[, i] <- column_i
}
return(df)
}Run the function
names(vcf_noinvar)[1:10]## [1] "sample" "pop" "super_pop" "sex" "lat" "lng"
## [7] "X3" "X9" "X10" "X11"vcf_noNA <- vcf_noinvar
vcf_noNA[,-c(1:6)] <- mean_imputation(vcf_noinvar[,-c(1:6)])## This make take some time ...##Prepare for PCA
Scale data
# new copy of data
vcf_scaled <- vcf_noNA
# scale
vcf_scaled[,-c(1:6)] <- scale(vcf_noNA[,-c(1:6)])
dim(vcf_scaled)## [1] 2504 5704write to csv
write.csv(vcf_scaled, file = "vcf_data.csv",
row.names = F)
list.files(pattern = "csv")## [1] "1000genomes_people_info2-1.csv" "for_pca.csv"
## [3] "gwas_pheno_env.csv" "pheno.csv"
## [5] "SNPs_cleaned.csv" "vcf_data.csv"
## [7] "vcf_num_df.csv" "vcf_num_df2.csv"
## [9] "vcf_num.csv" "walsh2017morphology.csv"