diff --git a/scripts/10_correction_batch_effects_unfiltered.Rmd b/scripts/10_correction_batch_effects_unfiltered.Rmd index fae0021..54ec777 100644 --- a/scripts/10_correction_batch_effects_unfiltered.Rmd +++ b/scripts/10_correction_batch_effects_unfiltered.Rmd @@ -51,7 +51,7 @@ knitr::opts_chunk$set(echo = FALSE) ```{r load data, include=FALSE} summary_table_preprocessing <- readRDS(paste0(user_choices$project_name, "/reports/summary_table_preprocessing.rds")) -load(paste0(user_choices$project_name, "/processed_data/Betas_clean_quantile_bmiq.Rdata")) +load(paste0(user_choices$project_name, "/processed_data/Betas_clean_unfiltered_quantile_bmiq.Rdata")) load(paste0(user_choices$project_name, "/processed_data/PhenoData_clean.Rdata")) load(paste0(user_choices$project_name, "/processed_data/RGSet_clean.Rdata")) phenotype_data <- readRDS(paste0(user_choices$project_name, "/processed_data/phenotype_data.rds")) @@ -60,7 +60,7 @@ phenotype_data <- readRDS(paste0(user_choices$project_name, "/processed_data/phe ## Correction of technical batch effects with ComBat ```{r, correction of first batch effect and PCA} -m_values <- apply(Betas_clean_quantile_bmiq, 2, function(x) log2((x)/(1-x))) # get M-values +m_values <- apply(Betas_clean_unfiltered_quantile_bmiq, 2, function(x) log2((x)/(1-x))) # get M-values if (correction_variable_1 == "PLEASE FILL IN"){ print("No first batch correction variable was specified by user, data remains unchanged") @@ -557,31 +557,31 @@ if (changing_status == "changed") { # function to transform M-values back to Beta values: expit2 <- function(x) 2^x/(1+2^x) - Betas_clean_quantile_bmiq_combated <- expit2(corrected_data) - save(Betas_clean_quantile_bmiq_combated, - file = paste0(user_choices$project_name, "/processed_data/Betas_clean_quantile_bmiq_combated.Rdata")) + Betas_clean_unfiltered_quantile_bmiq_combated <- expit2(corrected_data) + save(Betas_clean_unfiltered_quantile_bmiq_combated, + file = paste0(user_choices$project_name, "/processed_data/Betas_clean_unfiltered_quantile_bmiq_combated.Rdata")) # save an ExpressionSet - if (all.equal(colnames(Betas_clean_quantile_bmiq_combated),rownames(PhenoData_clean))){ + if (all.equal(colnames(Betas_clean_unfiltered_quantile_bmiq_combated),rownames(PhenoData_clean))){ annotated_PhenoData_clean <- new("AnnotatedDataFrame", data = as.data.frame(PhenoData_clean)) # required for ExpressionSet - Betas_clean_quantile_bmiq_combated_ExprSet = new("ExpressionSet", - exprs = as.matrix(Betas_clean_quantile_bmiq_combated), + Betas_clean_unfiltered_quantile_bmiq_combated_ExprSet = new("ExpressionSet", + exprs = as.matrix(Betas_clean_unfiltered_quantile_bmiq_combated), phenoData = annotated_PhenoData_clean) - save(Betas_clean_quantile_bmiq_combated_ExprSet, - file = paste0(user_choices$project_name, "/processed_data/Betas_clean_quantile_bmiq_combated_ExprSet.Rdata")) - RGSet_clean_final <- RGSet_clean[ ,(colnames(RGSet_clean) %in% sampleNames(Betas_clean_quantile_bmiq_combated_ExprSet))] - save(RGSet_clean_final, file = "/processed_data/finalData/RGSet_clean_final.Rdata") + save(Betas_clean_unfiltered_quantile_bmiq_combated_ExprSet, + file = paste0(user_choices$project_name, "/processed_data/Betas_clean_unfiltered_quantile_bmiq_combated_ExprSet.Rdata")) + RGSet_clean_unfiltered_final <- RGSet_clean[ ,(colnames(RGSet_clean) %in% sampleNames(Betas_clean_unfiltered_quantile_bmiq_combated_ExprSet))] + save(RGSet_clean_unfiltered_final, file = "/processed_data/finalData/RGSet_clean_unfiltered_final.Rdata") } } else { - if (all.equal(colnames(Betas_clean_quantile_bmiq),rownames(PhenoData_clean))){ + if (all.equal(colnames(Betas_clean_unfiltered_quantile_bmiq),rownames(PhenoData_clean))){ annotated_PhenoData_clean <- new("AnnotatedDataFrame", data = as.data.frame(PhenoData_clean)) # required for ExpressionSet - Betas_clean_quantile_bmiq_ExprSet = new("ExpressionSet", - exprs = as.matrix(Betas_clean_quantile_bmiq), + Betas_clean_unfiltered_quantile_bmiq_ExprSet = new("ExpressionSet", + exprs = as.matrix(Betas_clean_unfiltered_quantile_bmiq), phenoData = annotated_PhenoData_clean) - save(Betas_clean_quantile_bmiq_ExprSet, - file = paste0(user_choices$project_name, "/processed_data/Betas_clean_quantile_bmiq_ExprSet.Rdata")) - RGSet_clean_final <- RGSet_clean[ ,(colnames(RGSet_clean) %in% sampleNames(Betas_clean_quantile_bmiq_ExprSet))] - save(RGSet_clean_final, file = "/processed_data/finalData/RGSet_clean_final.Rdata") + save(Betas_clean_unfiltered_quantile_bmiq_ExprSet, + file = paste0(user_choices$project_name, "/processed_data/Betas_clean_unfiltered_quantile_bmiq_ExprSet.Rdata")) + RGSet_clean_unfiltered_final <- RGSet_clean[ ,(colnames(RGSet_clean) %in% sampleNames(Betas_clean_unfiltered_quantile_bmiq_ExprSet))] + save(RGSet_clean_unfiltered_final, file = "/processed_data/finalData/RGSet_clean_unfiltered_final.Rdata") } } ``` @@ -590,8 +590,8 @@ if (changing_status == "changed") { ```{r, beta values distribution report after batch effect correction} if (changing_status == "changed") { -densityPlot(Betas_clean_quantile_bmiq_combated, sampGroups = PhenoData_clean$Slide, - legend = FALSE, pal = "darkblue", main = "Post Batch Effect Correction - Normalized Unfiltered Beta", xlab = "Beta Value") +densityPlot(Betas_clean_unfiltered_quantile_bmiq_combated, sampGroups = PhenoData_clean$Slide, + legend = FALSE, pal = "darkblue", main = "Post Batch Effect Correction - Unfiltered and Normalized Beta", xlab = "Beta Value") } else { print("No data corrections were performed, therefore no new report available") } @@ -599,21 +599,21 @@ densityPlot(Betas_clean_quantile_bmiq_combated, sampGroups = PhenoData_clean$Sli ```{r, save beta values distribution report after batch effect correction, include=FALSE} if (changing_status == "changed") { -pdf(paste0(user_choices$project_name, "/reports/beta_densities_quantile_normalized_combated.pdf")) -densityPlot(Betas_clean_quantile_bmiq_combated, sampGroups = PhenoData_clean$Slide, - legend = FALSE, pal = "darkblue", main = "Post Batch Effect Correction - Normalized Unfilterd Beta", xlab = "Beta Value") +pdf(paste0(user_choices$project_name, "/reports/beta_densities_unfiltered_quantile_normalized_combated.pdf")) +densityPlot(Betas_clean_unfiltered_quantile_bmiq_combated, sampGroups = PhenoData_clean$Slide, + legend = FALSE, pal = "darkblue", main = "Post Batch Effect Correction - Unfiltered and Normalized Beta", xlab = "Beta Value") dev.off() } ``` ```{r, update summary table file for preprocessing steps, include=FALSE} if (changing_status == "changed") { - dim_Betas_combated <- dim(Betas_clean_quantile_bmiq_combated) + dim_Betas_unfiltered_combated <- dim(Betas_clean_unfiltered_quantile_bmiq_combated) step_number <- "12" step <- "Batch Correction - Unfiltered" data_class <- "Batch Correction - Unfiltered" - samples <- dim_Betas_combated[2] - probes <- dim_Betas_combated[1] + samples <- dim_Betas_unfiltered_combated[2] + probes <- dim_Betas_unfiltered_combated[1] table_preprocessing_adding <- data.frame(step_number, step, data_class, samples, probes) summary_table_preprocessing <- bind_rows(summary_table_preprocessing, table_preprocessing_adding) @@ -628,11 +628,11 @@ if (changing_status == "changed") { cat("PCA results and M-values in the different correction steps and final batch-corrected dataset were saved to the processed_data folder", sep = " threshold set by user) @@ -53,15 +53,15 @@ We will remove probes that: ## Removal of failed probes in one or more samples ```{r filter out failed probes, include=FALSE} -# ensure probes are in the same order in cleaned-normalized beta and cleaned detP objects, then subset -detP_clean_filtered <- detP_clean[match(rownames(Betas_clean_quantile_bmiq),rownames(detP_clean)),] -keep_betas <- rowSums(detP_clean_filtered < 0.01) == ncol(Betas_clean_quantile_bmiq) -Betas_clean_quantile_bmiq_filtered <- Betas_clean_quantile_bmiq[keep_betas,] - -# ensure probes are in the same order in cleaned-normalized RGChannelSet and detP objects, then subset -detP_clean_filtered2 <- detP_clean[match(featureNames(RGSet_clean_quantile),rownames(detP_clean)),] -keep_rgset <- rowSums(detP_clean_filtered2 < 0.01) == ncol(RGSet_clean_quantile) -RGSet_clean_quantile_filtered <- RGSet_clean_quantile[keep_rgset,] +# ensure probes are in the same order in cleaned beta and cleaned detP objects, then subset +detP_clean_filtered <- detP_clean[match(rownames(Betas_clean),rownames(detP_clean)),] +keep_betas <- rowSums(detP_clean_filtered < 0.01) == ncol(Betas_clean) +Betas_clean_filtered <- Betas_clean[keep_betas,] + +# ensure probes are in the same order in cleaned RGChannelSet and detP objects, then subset +detP_clean_filtered2 <- detP_clean[match(featureNames(RGSet_clean),rownames(detP_clean)),] +keep_rgset <- rowSums(detP_clean_filtered2 < 0.01) == ncol(RGSet_clean) +RGSet_clean_filtered <- RGSet_clean[keep_rgset,] ``` ```{r, info output for user failed probes, results='asis'} @@ -69,8 +69,8 @@ keep_betas_df <- as.data.frame(keep_betas) cat(paste0(keep_betas_df %>% filter(keep_betas == FALSE) %>% nrow(), " probes failed in one or more samples"), sep = "
\n") -dim_RGSet_filtered <- dim(RGSet_clean_quantile_filtered) -dim_Betas_filtered <- dim(Betas_clean_quantile_bmiq_filtered) +dim_RGSet_filtered <- dim(RGSet_clean_filtered) +dim_Betas_filtered <- dim(Betas_clean_filtered) step_number <- c("4", "4") step <- c("Filter technical", "Filter technical") data_class <- c("RGSet", "Betas") @@ -86,11 +86,11 @@ summary_table_preprocessing <- bind_rows(summary_table_preprocessing, table_prep ### X chromosome ```{r filter out probes on x chromosomes, include=FALSE} -keep_betas <- !(rownames(Betas_clean_quantile_bmiq_filtered) %in% annotations_clean$Name[annotations_clean$chr == "chrX"]) -Betas_clean_quantile_bmiq_filtered <- Betas_clean_quantile_bmiq_filtered[keep_betas,] +keep_betas <- !(rownames(Betas_clean_filtered) %in% annotations_clean$Name[annotations_clean$chr == "chrX"]) +Betas_clean_filtered <- Betas_clean_filtered[keep_betas,] -keep_rgset <- !(featureNames(RGSet_clean_quantile_filtered) %in% annotations_clean$Name[annotations_clean$chr == "chrX"]) -RGSet_clean_quantile_filtered <- RGSet_clean_quantile_filtered[keep_rgset,] +keep_rgset <- !(featureNames(RGSet_clean_filtered) %in% annotations_clean$Name[annotations_clean$chr == "chrX"]) +RGSet_clean_filtered <- RGSet_clean_filtered[keep_rgset,] ``` ```{r, info output for user x chromosomes, results='asis'} @@ -98,8 +98,8 @@ keep_betas_df <- as.data.frame(keep_betas) cat(paste0(keep_betas_df %>% filter(keep_betas == FALSE) %>% nrow(), " probes were on the X chromosome"), sep = "
\n") -dim_RGSet_filtered <- dim(RGSet_clean_quantile_filtered) -dim_Betas_filtered <- dim(Betas_clean_quantile_bmiq_filtered) +dim_RGSet_filtered <- dim(RGSet_clean_filtered) +dim_Betas_filtered <- dim(Betas_clean_filtered) step_number <- c("5", "5") step <- c("Filter X Chromosome", "Filter X Chromosome") data_class <- c("RGSet", "Betas") @@ -113,11 +113,11 @@ summary_table_preprocessing <- bind_rows(summary_table_preprocessing, table_prep ### Y chromosome ```{r filter out probes on y chromosomes, include=FALSE} -keep_betas <- !(rownames(Betas_clean_quantile_bmiq_filtered) %in% annotations_clean$Name[annotations_clean$chr == "chrY"]) -Betas_clean_quantile_bmiq_filtered <- Betas_clean_quantile_bmiq_filtered[keep_betas,] +keep_betas <- !(rownames(Betas_clean_filtered) %in% annotations_clean$Name[annotations_clean$chr == "chrY"]) +Betas_clean_filtered <- Betas_clean_filtered[keep_betas,] -keep_rgset <- !(featureNames(RGSet_clean_quantile_filtered) %in% annotations_clean$Name[annotations_clean$chr == "chrY"]) -RGSet_clean_quantile_filtered <- RGSet_clean_quantile_filtered[keep_rgset,] +keep_rgset <- !(featureNames(RGSet_clean_filtered) %in% annotations_clean$Name[annotations_clean$chr == "chrY"]) +RGSet_clean_filtered <- RGSet_clean_filtered[keep_rgset,] ``` ```{r, info output for user y chromosomes, results='asis'} @@ -125,8 +125,8 @@ keep_betas_df <- as.data.frame(keep_betas) cat(paste0(keep_betas_df %>% filter(keep_betas == FALSE) %>% nrow(), " probes were on the Y chromosome"), sep = "
\n") -dim_RGSet_filtered <- dim(RGSet_clean_quantile_filtered) -dim_Betas_filtered <- dim(Betas_clean_quantile_bmiq_filtered) +dim_RGSet_filtered <- dim(RGSet_clean_filtered) +dim_Betas_filtered <- dim(Betas_clean_filtered) step_number <- c("6", "6") step <- c("Filter Y Chromosome", "Filter Y Chromosome") data_class <- c("RGSet", "Betas") @@ -141,16 +141,16 @@ summary_table_preprocessing <- bind_rows(summary_table_preprocessing, table_prep ```{r filter out probes affected by snps, include=FALSE} if (user_choices$array_type == "v1") { - RGSet_clean_quantile_filtered <- dropLociWithSnps(RGSet_clean_quantile_filtered) # function requires a GenomicRatioSet + RGSet_clean_filtered <- dropLociWithSnps(RGSet_clean_filtered) # function requires a GenomicRatioSet - Betas_clean_quantile_bmiq_filtered <- Betas_clean_quantile_bmiq_filtered[rownames(Betas_clean_quantile_bmiq_filtered) %in% featureNames(RGSet_clean_quantile_filtered),] + Betas_clean_filtered <- Betas_clean_filtered[rownames(Betas_clean_filtered) %in% featureNames(RGSet_clean_filtered),] } ``` ```{r, info output for user probes affected by snps, results='asis'} if (user_choices$array_type == "v1") { - dim_RGSet_filtered <- dim(RGSet_clean_quantile_filtered) - dim_Betas_filtered <- dim(Betas_clean_quantile_bmiq_filtered) + dim_RGSet_filtered <- dim(RGSet_clean_filtered) + dim_Betas_filtered <- dim(Betas_clean_filtered) step_number <- c("7", "7") step <- c("Filter SNP Affected", "Filter SNP Affected") data_class <- c("RGSet", "Betas") @@ -183,11 +183,11 @@ if (user_choices$array_type == "v1") { index <- which(annot2$sex == "Exclude" | annot2$CRSH == "Exclude" | annot2$EUR == "Exclude") exclude_Chen <- annot2[index,] - keep_betas <- !(rownames(Betas_clean_quantile_bmiq_filtered) %in% exclude_Chen$Name) - Betas_clean_quantile_bmiq_filtered <- Betas_clean_quantile_bmiq_filtered[keep_betas,] + keep_betas <- !(rownames(Betas_clean_filtered) %in% exclude_Chen$Name) + Betas_clean_filtered <- Betas_clean_filtered[keep_betas,] - keep_rgset <- !(featureNames(RGSet_clean_quantile_filtered) %in% exclude_Chen$Name) - RGSet_clean_quantile_filtered <- RGSet_clean_quantile_filtered[keep_rgset,] + keep_rgset <- !(featureNames(RGSet_clean_filtered) %in% exclude_Chen$Name) + RGSet_clean_filtered <- RGSet_clean_filtered[keep_rgset,] } ``` @@ -197,8 +197,8 @@ if (user_choices$array_type == "v1") { cat(paste0(keep_betas_df %>% filter(keep_betas == FALSE) %>% nrow(), " probes were cross-hybridizing according to Chen et al."), sep = "
\n") - dim_RGSet_filtered <- dim(RGSet_clean_quantile_filtered) - dim_Betas_filtered <- dim(Betas_clean_quantile_bmiq_filtered) + dim_RGSet_filtered <- dim(RGSet_clean_filtered) + dim_Betas_filtered <- dim(Betas_clean_filtered) step_number <- c("8", "9") step <- c("Filter Cross-Hybridizing Chen", "Filter Cross-Hybridizing Chen") data_class <- c("RGSet", "Betas") @@ -216,11 +216,11 @@ if (user_choices$array_type == "v1") { if (user_choices$array_type == "v1") { exclude_McCartney <- read.table("data/probes_crosshybridizing_EPIC.txt", sep = "", header = F) - keep_betas <- !(rownames(Betas_clean_quantile_bmiq_filtered) %in% exclude_McCartney$V1) - Betas_clean_quantile_bmiq_filtered <- Betas_clean_quantile_bmiq_filtered[keep_betas,] + keep_betas <- !(rownames(Betas_clean_filtered) %in% exclude_McCartney$V1) + Betas_clean_filtered <- Betas_clean_filtered[keep_betas,] - keep_rgset <- !(featureNames(RGSet_clean_quantile_filtered) %in% exclude_McCartney$V1) - RGSet_clean_quantile_filtered <- RGSet_clean_quantile_filtered[keep_rgset,] + keep_rgset <- !(featureNames(RGSet_clean_filtered) %in% exclude_McCartney$V1) + RGSet_clean_filtered <- RGSet_clean_filtered[keep_rgset,] } ``` @@ -230,8 +230,8 @@ if (user_choices$array_type == "v1") { cat(paste0(keep_betas_df %>% filter(keep_betas == FALSE) %>% nrow(), " probes may be cross-hybridizing according to McCartney et al."), sep = "
\n") - dim_RGSet_filtered <- dim(RGSet_clean_quantile_filtered) - dim_Betas_filtered <- dim(Betas_clean_quantile_bmiq_filtered) + dim_RGSet_filtered <- dim(RGSet_clean_filtered) + dim_Betas_filtered <- dim(Betas_clean_filtered) step_number <- c("9", "9") step <- c("Filter Cross-Hybridizing McCartney", "Filter Cross-Hybridizing McCartney") data_class <- c("RGSet", "Betas") @@ -251,11 +251,11 @@ if (user_choices$array_type == "v1") { index <- which(exclude_polymorphic$EUR_AF > 0.05) exclude_polymorphic <- exclude_polymorphic[index,] - keep_betas <- !(rownames(Betas_clean_quantile_bmiq_filtered) %in% exclude_polymorphic$IlmnID) - Betas_clean_quantile_bmiq_filtered <- Betas_clean_quantile_bmiq_filtered[keep_betas,] + keep_betas <- !(rownames(Betas_clean_filtered) %in% exclude_polymorphic$IlmnID) + Betas_clean_filtered <- Betas_clean_filtered[keep_betas,] - keep_rgset <- !(featureNames(RGSet_clean_quantile_filtered) %in% exclude_polymorphic$IlmnID) - RGSet_clean_quantile_filtered <- RGSet_clean_quantile_filtered[keep_rgset,] + keep_rgset <- !(featureNames(RGSet_clean_filtered) %in% exclude_polymorphic$IlmnID) + RGSet_clean_filtered <- RGSet_clean_filtered[keep_rgset,] } ``` @@ -265,8 +265,8 @@ if (user_choices$array_type == "v1") { cat(paste0(keep_betas_df %>% filter(keep_betas == FALSE) %>% nrow(), " probes were polymorphic"), sep = "
\n") - dim_RGSet_filtered <- dim(RGSet_clean_quantile_filtered) - dim_Betas_filtered <- dim(Betas_clean_quantile_bmiq_filtered) + dim_RGSet_filtered <- dim(RGSet_clean_filtered) + dim_Betas_filtered <- dim(Betas_clean_filtered) step_number <- c("10", "10") step <- c("Filter Polymorphic", "Filter Polymorphic") data_class <- c("RGSet", "Betas") @@ -280,27 +280,27 @@ saveRDS(summary_table_preprocessing, paste0(user_choices$project_name, "/reports ``` ```{r, save data, include=FALSE} -save(Betas_clean_quantile_bmiq_filtered, file = paste0(user_choices$project_name, "/processed_data/Betas_clean_quantile_bmiq_filtered.Rdata")) -save(RGSet_clean_quantile_filtered, file = paste0(user_choices$project_name, "/processed_data/RGSet_clean_quantile_filtered.Rdata")) +save(Betas_clean_filtered, file = paste0(user_choices$project_name, "/processed_data/Betas_clean_filtered.Rdata")) +save(RGSet_clean_filtered, file = paste0(user_choices$project_name, "/processed_data/RGSet_clean_filtered.Rdata")) -Betas_clean_quantile_filtered <- getBeta(RGSet_clean_quantile_filtered) # beta values -save(Betas_clean_quantile_filtered, file = paste0(user_choices$project_name, "/processed_data/Betas_clean_quantile_filtered.Rdata")) +Betas_clean_filtered <- getBeta(RGSet_clean_filtered) # beta values +save(Betas_clean_filtered, file = paste0(user_choices$project_name, "/processed_data/Betas_clean_filtered.Rdata")) -Ms_clean_quantile_filtered <- getM(RGSet_clean_quantile_filtered) -save(Ms_clean_quantile_filtered, file = paste0(user_choices$project_name, "/processed_data/Ms_clean_quantile_filtered.Rdata")) +Ms_clean_filtered <- getM(RGSet_clean_filtered) +save(Ms_clean_filtered, file = paste0(user_choices$project_name, "/processed_data/Ms_clean_filtered.Rdata")) ``` ## Beta values distribution report ```{r, beta values distribution report after excluding the poorly-detected probes} -densityPlot(Betas_clean_quantile_bmiq_filtered, sampGroups = PhenoData_clean$Slide, - legend = FALSE, pal = "darkblue", main = "Post Filtering - Normalized Beta", xlab = "Beta Value") +densityPlot(Betas_clean_filtered, sampGroups = PhenoData_clean$Slide, + legend = FALSE, pal = "darkblue", main = "Post Filtering Beta", xlab = "Beta Value") ``` ```{r, save beta values distribution report after excluding the poorly-detected probes, include=FALSE} -pdf(paste0(user_choices$project_name, "/reports/beta_densities_quantile_normalized_filtered.pdf")) -densityPlot(Betas_clean_quantile_bmiq_filtered, sampGroups = PhenoData_clean$Slide, - legend = FALSE, pal = "darkblue", main = "Post Filtering - Normalized Beta", xlab = "Beta Value") +pdf(paste0(user_choices$project_name, "/reports/beta_densities_filtered.pdf")) +densityPlot(Betas_clean_filtered, sampGroups = PhenoData_clean$Slide, + legend = FALSE, pal = "darkblue", main = "Post Filtering Beta", xlab = "Beta Value") dev.off() ``` @@ -320,7 +320,7 @@ summary_table_preprocessing %>% paged_table() connection <- file(paste0(user_choices$personal_path, "/", user_choices$project_name, "/session_info_all.txt"), "a+") writeLines("######################################################################", connection) -writeLines("############################ Script 7: #############################", connection) +writeLines("############################ Script 6: #############################", connection) writeLines("######################################################################", connection) sessioninfo <- sessionInfo() writeLines("\nR Version:", connection) diff --git a/scripts/7_normalization.Rmd b/scripts/7_normalization.Rmd index f10737e..ef4718a 100644 --- a/scripts/7_normalization.Rmd +++ b/scripts/7_normalization.Rmd @@ -1,9 +1,11 @@ --- -title: "Script 6: Normalization" +title: "Script 7: Normalization" date: '`r format(Sys.time(), "%d %B, %Y")`' output: html_document --- +->>>>>>> To do: Do normalizetion also for _unfiltered data!!!! + ```{r runtime start, include=FALSE} start_time <- Sys.time() ``` @@ -40,8 +42,8 @@ if (user_choices$array_type == "v2") { ``` ```{r load data, include=FALSE} -load(paste0(user_choices$project_name, "/processed_data/RGSet_clean.Rdata")) -load(paste0(user_choices$project_name, "/processed_data/Betas_clean.Rdata")) +load(paste0(user_choices$project_name, "/processed_data/RGSet_clean_filtered.Rdata")) +load(paste0(user_choices$project_name, "/processed_data/Betas_clean_filtered.Rdata")) load(paste0(user_choices$project_name, "/processed_data/PhenoData_clean.Rdata")) ``` @@ -62,30 +64,30 @@ data for further investigation. These files are saved to reports folder. The user can further alter the normalization methods if needed (e.g. include noob etc.) ```{r get annotation of probes, include=FALSE} -annotations_clean = getAnnotation(RGSet_clean) -save(annotations_clean, file = paste0(user_choices$project_name, "/reports/annotations_clean.Rdata")) +annotations_clean_filtered = getAnnotation(RGSet_clean_filtered) +save(annotations_clean_filtered, file = paste0(user_choices$project_name, "/reports/annotations_clean_filtered.Rdata")) ``` ## Normalization of data performed ```{r normalize probes, include=FALSE} -RGSet_clean_quantile <- preprocessQuantile(RGSet_clean) -save(RGSet_clean_quantile, file = paste0(user_choices$project_name, - "/processed_data/RGSet_clean_quantile.Rdata")) # output: GenomicRatioSet - -Betas_clean_quantile <- getBeta(RGSet_clean_quantile) -save(Betas_clean_quantile, file = paste0(user_choices$project_name, - "/processed_data/Betas_clean_quantile.Rdata")) -Mset_clean_qunatile <- getM(RGSet_clean_quantile) -save(Mset_clean_qunatile, file = paste0(user_choices$project_name, - "/processed_data/Mset_clean_qunatile.Rdata")) +RGSet_clean_filtered_quantile <- preprocessQuantile(RGSet_clean_filtered) +save(RGSet_clean_filtered_quantile, file = paste0(user_choices$project_name, + "/processed_data/RGSet_clean_filtered_quantile.Rdata")) # output: GenomicRatioSet + +Betas_clean_filtered_quantile <- getBeta(RGSet_clean_filtered_quantile) +save(Betas_clean_filtered_quantile, file = paste0(user_choices$project_name, + "/processed_data/Betas_clean_filtered_quantile.Rdata")) +Mset_clean_filtered_qunatile <- getM(RGSet_clean_filtered_quantile) +save(Mset_clean_filtered_qunatile, file = paste0(user_choices$project_name, + "/processed_data/Mset_clean_filtered_qunatile.Rdata")) # further normalization with BMIQ: -probeType <- as.data.frame(annotations_clean[rownames(Betas_clean_quantile),c("Name","Type")]) +probeType <- as.data.frame(annotations_clean_filtered[rownames(Betas_clean_filtered_quantile),c("Name","Type")]) probeType$probeType = ifelse(probeType$Type %in% "I", 1, 2) -Betas_clean_quantile_bmiq <- apply(Betas_clean_quantile[,1:length(colnames(Betas_clean_quantile))],2, +Betas_clean_filtered_quantile_bmiq <- apply(Betas_clean_filtered_quantile[,1:length(colnames(Betas_clean_filtered_quantile))],2, function(a) BMIQ(a,probeType$probeType,plots=FALSE)$nbeta) # sourced from script "BMIQ_1.6_Teschendorff.R" -save(Betas_clean_quantile_bmiq, file = paste0(user_choices$project_name, - "/processed_data/Betas_clean_quantile_bmiq.Rdata")) +save(Betas_clean_filtered_quantile_bmiq, file = paste0(user_choices$project_name, + "/processed_data/Betas_clean_filtered_quantile_bmiq.Rdata")) ``` ## Beta values distribution report @@ -93,21 +95,25 @@ save(Betas_clean_quantile_bmiq, file = paste0(user_choices$project_name, ```{r beta values distribution} densityPlot(Betas_clean, sampGroups = PhenoData_clean$Slide, legend = FALSE, pal = "darkblue", main = "Raw Betas", xlab = "Beta Value") -densityPlot(Betas_clean_quantile, sampGroups = PhenoData_clean$Slide, legend = FALSE, - pal = "darkblue", main = "Quantile Adjusted Betas", xlab = "Beta Value") -densityPlot(Betas_clean_quantile_bmiq, sampGroups = PhenoData_clean$Slide, legend = FALSE, - pal = "darkblue", main = "Quantile-BMIQ Adjusted Betas", xlab = "Beta Value") +densityPlot(Betas_clean_filtered, sampGroups = PhenoData_clean$Slide, legend = FALSE, + pal = "darkblue", main = "Filtered Betas", xlab = "Beta Value") +densityPlot(Betas_clean_filtered_quantile, sampGroups = PhenoData_clean$Slide, legend = FALSE, + pal = "darkblue", main = "Filtered and Quantile Adjusted Betas", xlab = "Beta Value") +densityPlot(Betas_clean_filtered_quantile_bmiq, sampGroups = PhenoData_clean$Slide, legend = FALSE, + pal = "darkblue", main = "Filtered and Quantile-BMIQ Adjusted Betas", xlab = "Beta Value") ``` ```{r, save reports as pdf, include=FALSE} pdf(file = paste0(user_choices$project_name, - "/reports/beta_distributions_raw_normalized.pdf")) + "/reports/beta_distributions_raw_filtered_normalized.pdf")) densityPlot(Betas_clean, sampGroups = PhenoData_clean$Slide, legend = FALSE, pal = "darkblue", main = "Raw Betas", xlab = "Beta Value") -densityPlot(Betas_clean_quantile, sampGroups = PhenoData_clean$Slide, legend = FALSE, - pal = "darkblue", main = "Quantile Adjusted Betas", xlab = "Beta Value") -densityPlot(Betas_clean_quantile_bmiq, sampGroups = PhenoData_clean$Slide, legend = FALSE, - pal = "darkblue", main = "Quantile-BMIQ Adjusted Betas", xlab = "Beta Value") +densityPlot(Betas_clean_filtered, sampGroups = PhenoData_clean$Slide, legend = FALSE, + pal = "darkblue", main = "Filtered Betas", xlab = "Beta Value") +densityPlot(Betas_clean_filtered_quantile, sampGroups = PhenoData_clean$Slide, legend = FALSE, + pal = "darkblue", main = "Filtered and Quantile Adjusted Betas", xlab = "Beta Value") +densityPlot(Betas_clean_filtered_quantile_bmiq, sampGroups = PhenoData_clean$Slide, legend = FALSE, + pal = "darkblue", main = "Filtered and Quantile-BMIQ Adjusted Betas", xlab = "Beta Value") dev.off() # additional plots to look for outliers in distribution: @@ -119,12 +125,12 @@ for (i in 1:ncol(Betas_clean)) { } dev.off() -names_quantileBetas <- colnames(Betas_clean_quantile) -pdf(paste0(user_choices$project_name, "/reports/beta_densities_quantile_normalized.pdf")) -for (i in 1:ncol(Betas_clean_quantile)) { - i_mat <- as.matrix(Betas_clean_quantile[ ,i]) - name <- colnames(Betas_clean_quantile[,i]) - densityPlot(i_mat, pal = "darkblue", main = names_quantileBetas[i]) +names_filtered_quantileBetas <- colnames(Betas_clean_filtered_quantile) +pdf(paste0(user_choices$project_name, "/reports/beta_densities_filtered_quantile_normalized.pdf")) +for (i in 1:ncol(Betas_clean_filtered_quantile)) { + i_mat <- as.matrix(Betas_clean_filtered_quantile[ ,i]) + name <- colnames(Betas_clean_filtered_quantile[,i]) + densityPlot(i_mat, pal = "darkblue", main = names_filtered_quantileBetas[i]) } dev.off() ``` @@ -134,7 +140,7 @@ dev.off() ```{r, info output for user 1, results='asis'} cat("Data underwent stratified quantile normalization followed by BMIQ and were saved to the processed data folder. ", sep = "
\n") -if (length(which(is.nan(Betas_clean_quantile_bmiq))) == 0) { +if (length(which(is.nan(Betas_clean_filtered_quantile_bmiq))) == 0) { cat("There were no NAs after normalization", sep = "
\n") } else { cat("NAs were produced during normalization, please re-check!", sep = "
\n") @@ -149,7 +155,7 @@ cat("Session info text file was updated", sep = "
\n") connection <- file(paste0(user_choices$personal_path, "/", user_choices$project_name, "/session_info_all.txt"), "a+") writeLines("######################################################################", connection) -writeLines("############################ Script 6: #############################", connection) +writeLines("############################ Script 7: #############################", connection) writeLines("######################################################################", connection) sessioninfo <- sessionInfo() writeLines("\nR Version:", connection) diff --git a/scripts/8_removal_of_outliers_detect_batch_effects.Rmd b/scripts/8_removal_of_outliers_detect_batch_effects.Rmd index 6b325a7..e7719c0 100644 --- a/scripts/8_removal_of_outliers_detect_batch_effects.Rmd +++ b/scripts/8_removal_of_outliers_detect_batch_effects.Rmd @@ -69,9 +69,9 @@ knitr::opts_chunk$set(echo = FALSE) ```{r load data, include=FALSE} summary_table_preprocessing <- readRDS(paste0(user_choices$project_name, "/reports/summary_table_preprocessing.rds")) -load(paste0(user_choices$project_name, "/processed_data/Betas_clean_quantile.Rdata")) -load(paste0(user_choices$project_name, "/processed_data/Betas_clean_quantile_bmiq.Rdata")) -load(paste0(user_choices$project_name, "/processed_data/Betas_clean_quantile_bmiq_filtered.Rdata")) +load(paste0(user_choices$project_name, "/processed_data/Betas_clean_filtered_quantile.Rdata")) +load(paste0(user_choices$project_name, "/processed_data/Betas_clean_unfiltered_quantile_bmiq.Rdata")) +load(paste0(user_choices$project_name, "/processed_data/Betas_clean_filtered_quantile_bmiq.Rdata")) load(paste0(user_choices$project_name, "/processed_data/PhenoData_clean.Rdata")) phenotype_data <- readRDS(paste0(user_choices$project_name, "/processed_data/phenotype_data.rds")) ``` @@ -91,7 +91,7 @@ rowVars <- function(x, na.rm=FALSE, dims=1, unbiased=TRUE, SumSquares=FALSE, two (rowSums(x^2, na.rm, dims) - rowSums(x, na.rm, dims)^2/N) / Nm1 } -m_values <- apply(Betas_clean_quantile_bmiq_filtered, 2, function(x) log2((x)/(1-x))) # get M-values +m_values <- apply(Betas_clean_filtered_quantile_bmiq, 2, function(x) log2((x)/(1-x))) # get M-values variance_m_values = as.matrix(rowVars(m_values)) # calculate variances per row if (identical(which(variance_m_values == 0), integer(0))){ @@ -102,7 +102,7 @@ if (identical(which(variance_m_values == 0), integer(0))){ variance_m_values[variance_m_values == 0] <- NA variance_m_values <- na.omit(variance_m_values) intersect <- intersect(rownames(variance_m_values), rownames(m_values)) - Betas_clean_quantile_bmiq_filtered_combat <- Betas_clean_quantile_bmiq_filtered[intersect,] # original name _batch + Betas_clean_filtered_quantile_bmiq_combat <- Betas_clean_filtered_quantile_bmiq[intersect,] # original name _batch m_values <- m_values[intersect,] } ``` @@ -227,11 +227,11 @@ dev.off() if(!identical(extreme_outliers_filtered %>% pull(personid),character(0))){ print("PCA is repeated after exclusion of extreme outliers") # remove outliers - Betas_clean_quantile_bmiq_filtered <- Betas_clean_quantile_bmiq_filtered[, !(colnames(Betas_clean_quantile_bmiq_filtered) %in% extreme_outliers_filtered)] + Betas_clean_filtered_quantile_bmiq <- Betas_clean_filtered_quantile_bmiq[, !(colnames(Betas_clean_filtered_quantile_bmiq) %in% extreme_outliers_filtered)] PhenoData_clean <- PhenoData_clean[, !(colnames(PhenoData_clean) %in% extreme_outliers_filtered)] # repeat all code above after outlier exclusion - m_values <- apply(Betas_clean_quantile_bmiq_filtered, 2, function(x) log2((x)/(1-x))) # get M-values + m_values <- apply(Betas_clean_filtered_quantile_bmiq, 2, function(x) log2((x)/(1-x))) # get M-values variance_m_values = as.matrix(rowVars(m_values)) # calculate variances per row if (identical(which(variance_m_values == 0), integer(0))){ @@ -242,7 +242,7 @@ if(!identical(extreme_outliers_filtered %>% pull(personid),character(0))){ variance_m_values[variance_m_values == 0] <- NA variance_m_values = na.omit(variance_m_values) intersect <- intersect(rownames(variance_m_values), rownames(m_values)) - Betas_clean_quantile_bmiq_filtered_combat <- Betas_clean_quantile_bmiq_filtered[intersect,] # original name _batch + Betas_clean_filtered_quantile_bmiq_combat <- Betas_clean_filtered_quantile_bmiq[intersect,] # original name _batch m_values = m_values[intersect,] } PCA_object <- prcomp(t(m_values), retx = T, center = T, scale. = T) @@ -437,7 +437,7 @@ if (sum(is.na(anova_lm_pvalue_df$P_Value)) == nrow(anova_lm_pvalue_df)){ ### Removal of probes with zero variance across samples ```{r check probe variance and exclude if 0 - unfilterd, results='asis'} -m_values <- apply(Betas_clean_quantile_bmiq, 2, function(x) log2((x)/(1-x))) # get M-values +m_values <- apply(Betas_clean_unfiltered_quantile_bmiq, 2, function(x) log2((x)/(1-x))) # get M-values variance_m_values = as.matrix(rowVars(m_values)) # calculate variances per row if (identical(which(variance_m_values == 0), integer(0))){ @@ -448,7 +448,7 @@ if (identical(which(variance_m_values == 0), integer(0))){ variance_m_values[variance_m_values == 0] <- NA variance_m_values <- na.omit(variance_m_values) intersect <- intersect(rownames(variance_m_values), rownames(m_values)) - Betas_clean_quantile_bmiq_combat <- Betas_clean_quantile_bmiq[intersect,] # original name _batch + Betas_clean_unfiltered_quantile_bmiq_combat <- Betas_clean_unfiltered_quantile_bmiq[intersect,] # original name _batch m_values <- m_values[intersect,] } ``` @@ -573,12 +573,12 @@ dev.off() if(!identical(extreme_outliers_unfiltered %>% pull(personid),character(0))){ print("PCA is repeated after exclusion of extreme outliers") # remove outliers - Betas_clean_quantile_bmiq <- Betas_clean_quantile_bmiq[, !(colnames(Betas_clean_quantile_bmiq) %in% extreme_outliers_unfiltered)] - Betas_clean_quantile <- Betas_clean_quantile[, !(colnames(Betas_clean_quantile) %in% extreme_outliers_unfiltered)] + Betas_clean_unfiltered_quantile_bmiq <- Betas_clean_unfiltered_quantile_bmiq[, !(colnames(Betas_clean_unfiltered_quantile_bmiq) %in% extreme_outliers_unfiltered)] + Betas_clean_unfiltered_quantile <- Betas_clean_unfiltered_quantile[, !(colnames(Betas_clean_unfiltered_quantile) %in% extreme_outliers_unfiltered)] PhenoData_clean <- PhenoData_clean[, !(colnames(PhenoData_clean) %in% extreme_outliers_unfiltered)] # repeat all code above after outlier exclusion - m_values <- apply(Betas_clean_quantile_bmiq, 2, function(x) log2((x)/(1-x))) # get M-values + m_values <- apply(Betas_clean_unfiltered_quantile_bmiq, 2, function(x) log2((x)/(1-x))) # get M-values variance_m_values = as.matrix(rowVars(m_values)) # calculate variances per row if (identical(which(variance_m_values == 0), integer(0))){ @@ -589,7 +589,7 @@ if(!identical(extreme_outliers_unfiltered %>% pull(personid),character(0))){ variance_m_values[variance_m_values == 0] <- NA variance_m_values = na.omit(variance_m_values) intersect <- intersect(rownames(variance_m_values), rownames(m_values)) - Betas_clean_quantile_bmiq_combat <- Betas_clean_quantile_bmiq[intersect,] # original name _batch + Betas_clean_unfiltered_quantile_bmiq_combat <- Betas_clean_unfiltered_quantile_bmiq[intersect,] # original name _batch m_values = m_values[intersect,] } PCA_object <- prcomp(t(m_values), retx = T, center = T, scale. = T) @@ -780,13 +780,13 @@ if (sum(is.na(anova_lm_pvalue_df$P_Value)) == nrow(anova_lm_pvalue_df)){ ``` ```{r, update summary table file for preprocessing steps, include=FALSE} -dim_Betas_filtered_no_outliers <- dim(Betas_clean_quantile_bmiq_filtered) -dim_Betas_no_outliers <- dim(Betas_clean_quantile_bmiq) +dim_Betas_filtered_no_outliers <- dim(Betas_clean_filtered_quantile_bmiq) +dim_Betas_unfiltered_no_outliers <- dim(Betas_clean_unfiltered_quantile_bmiq) step_number <- c("10", "10") step <- c("Outlier Removal", "Outlier Removal") data_class <- c("Betas filtered", "Betas unfiltered") -samples <- c(dim_Betas_filtered_no_outliers[2], dim_Betas_no_outliers[2]) -probes <- c(dim_Betas_filtered_no_outliers[1], dim_Betas_no_outliers[1]) +samples <- c(dim_Betas_filtered_no_outliers[2], dim_Betas_unfiltered_no_outliers[2]) +probes <- c(dim_Betas_filtered_no_outliers[1], dim_Betas_unfiltered_no_outliers[1]) table_preprocessing_adding <- data.frame(step_number, step, data_class, samples, probes) summary_table_preprocessing <- bind_rows(summary_table_preprocessing, table_preprocessing_adding) diff --git a/scripts/9_correction_batch_effects_filtered.Rmd b/scripts/9_correction_batch_effects_filtered.Rmd index 2081997..6f3ba58 100644 --- a/scripts/9_correction_batch_effects_filtered.Rmd +++ b/scripts/9_correction_batch_effects_filtered.Rmd @@ -62,7 +62,7 @@ knitr::opts_chunk$set(echo = FALSE) ```{r load data, include=FALSE} summary_table_preprocessing <- readRDS(paste0(user_choices$project_name, "/reports/summary_table_preprocessing.rds")) -load(paste0(user_choices$project_name, "/processed_data/Betas_clean_quantile_bmiq_filtered.Rdata")) +load(paste0(user_choices$project_name, "/processed_data/Betas_clean_filtered_quantile_bmiq.Rdata")) load(paste0(user_choices$project_name, "/processed_data/PhenoData_clean.Rdata")) phenotype_data <- readRDS(paste0(user_choices$project_name, "/processed_data/phenotype_data.rds")) ``` @@ -70,7 +70,7 @@ phenotype_data <- readRDS(paste0(user_choices$project_name, "/processed_data/phe ## Correction of technical batch effects with ComBat ```{r, correction of first batch effect and PCA} -m_values <- apply(Betas_clean_quantile_bmiq_filtered, 2, function(x) log2((x)/(1-x))) # get M-values +m_values <- apply(Betas_clean_filtered_quantile_bmiq, 2, function(x) log2((x)/(1-x))) # get M-values if (correction_variable_1 == "PLEASE FILL IN"){ print("No first batch correction variable was specified by user, data remains unchanged") @@ -566,27 +566,27 @@ if (changing_status == "changed") { # function to transform M-values back to Beta values: expit2 <- function(x) 2^x/(1+2^x) - Betas_clean_quantile_bmiq_filtered_combated <- expit2(corrected_data) - save(Betas_clean_quantile_bmiq_filtered_combated, - file = paste0(user_choices$project_name, "/processed_data/Betas_clean_quantile_bmiq_filtered_combated.Rdata")) + Betas_clean_filtered_quantile_bmiq_combated <- expit2(corrected_data) + save(Betas_clean_filtered_quantile_bmiq_combated, + file = paste0(user_choices$project_name, "/processed_data/Betas_clean_filtered_quantile_bmiq_combated.Rdata")) # save an ExpressionSet - if (all.equal(colnames(Betas_clean_quantile_bmiq_filtered_combated),rownames(PhenoData_clean))){ + if (all.equal(colnames(Betas_clean_filtered_quantile_bmiq_combated),rownames(PhenoData_clean))){ annotated_PhenoData_clean <- new("AnnotatedDataFrame", data = as.data.frame(PhenoData_clean)) # required for ExpressionSet - Betas_clean_quantile_bmiq_filtered_combated_ExprSet = new("ExpressionSet", - exprs = as.matrix(Betas_clean_quantile_bmiq_filtered_combated), + Betas_clean_filtered_quantile_bmiq_combated_ExprSet = new("ExpressionSet", + exprs = as.matrix(Betas_clean_filtered_quantile_bmiq_combated), phenoData = annotated_PhenoData_clean) - save(Betas_clean_quantile_bmiq_filtered_combated_ExprSet, - file = paste0(user_choices$project_name, "/processed_data/Betas_clean_quantile_bmiq_filtered_combated_ExprSet.Rdata")) + save(Betas_clean_filtered_quantile_bmiq_combated_ExprSet, + file = paste0(user_choices$project_name, "/processed_data/Betas_clean_filtered_quantile_bmiq_combated_ExprSet.Rdata")) } } else { - if (all.equal(colnames(Betas_clean_quantile_bmiq_filtered),rownames(PhenoData_clean))){ + if (all.equal(colnames(Betas_clean_filtered_quantile_bmiq),rownames(PhenoData_clean))){ annotated_PhenoData_clean <- new("AnnotatedDataFrame", data = as.data.frame(PhenoData_clean)) # required for ExpressionSet - Betas_clean_quantile_bmiq_filtered_ExprSet = new("ExpressionSet", - exprs = as.matrix(Betas_clean_quantile_bmiq_filtered), + Betas_clean_filtered_quantile_bmiq_ExprSet = new("ExpressionSet", + exprs = as.matrix(Betas_clean_filtered_quantile_bmiq), phenoData = annotated_PhenoData_clean) - save(Betas_clean_quantile_bmiq_filtered_ExprSet, - file = paste0(user_choices$project_name, "/processed_data/Betas_clean_quantile_bmiq_filtered_ExprSet.Rdata")) + save(Betas_clean_filtered_quantile_bmiq_ExprSet, + file = paste0(user_choices$project_name, "/processed_data/Betas_clean_filtered_quantile_bmiq_ExprSet.Rdata")) } } ``` @@ -595,8 +595,8 @@ if (changing_status == "changed") { ```{r, beta values distribution report after batch effect correction} if (changing_status == "changed") { -densityPlot(Betas_clean_quantile_bmiq_filtered_combated, sampGroups = PhenoData_clean$Slide, - legend = FALSE, pal = "darkblue", main = "Post Batch Effect Correction - Normalized Filtered Beta", xlab = "Beta Value") +densityPlot(Betas_clean_filtered_quantile_bmiq_combated, sampGroups = PhenoData_clean$Slide, + legend = FALSE, pal = "darkblue", main = "Post Batch Effect Correction - Filtered and Normalized Beta", xlab = "Beta Value") } else { print("No batch corrections were performed, therefore no new reports are available") } @@ -604,16 +604,16 @@ densityPlot(Betas_clean_quantile_bmiq_filtered_combated, sampGroups = PhenoData_ ```{r, save beta values distribution report after batch effect correction, include=FALSE} if (changing_status == "changed") { -pdf(paste0(user_choices$project_name, "/reports/beta_densities_quantile_normalized_filtered_combated.pdf")) -densityPlot(Betas_clean_quantile_bmiq_filtered_combated, sampGroups = PhenoData_clean$Slide, - legend = FALSE, pal = "darkblue", main = "Post Batch Effect Correction - Normalized Filtered Beta", xlab = "Beta Value") +pdf(paste0(user_choices$project_name, "/reports/beta_densities_filtered_quantile_normalized_combated.pdf")) +densityPlot(Betas_clean_filtered_quantile_bmiq_combated, sampGroups = PhenoData_clean$Slide, + legend = FALSE, pal = "darkblue", main = "Post Batch Effect Correction - Filtered and Normalized Beta", xlab = "Beta Value") dev.off() } ``` ```{r, update summary table file for preprocessing steps, include=FALSE} if (changing_status == "changed") { - dim_Betas_filtered_combated <- dim(Betas_clean_quantile_bmiq_filtered_combated) + dim_Betas_filtered_combated <- dim(Betas_clean_filtered_quantile_bmiq_combated) step_number <- "11" step <- "Batch Correction - Filtered" data_class <- "Batch Correction - Filtered" @@ -633,11 +633,11 @@ if (changing_status == "changed") { cat("PCA results and M-values in the different correction steps and final batch-corrected dataset were saved to the processed_data folder", sep = "