Skip to content

Commit

Permalink
naming adapted to scripts sequence change
Browse files Browse the repository at this point in the history
  • Loading branch information
@NatanYusupov
NatanYusupov committed Mar 15, 2024
1 parent d1c9c4b commit f63ed96
Show file tree
Hide file tree
Showing 5 changed files with 176 additions and 170 deletions.
62 changes: 31 additions & 31 deletions scripts/10_correction_batch_effects_unfiltered.Rmd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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"))
Expand All @@ -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")
Expand Down Expand Up @@ -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")
}
}
```
Expand All @@ -590,30 +590,30 @@ 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")
}
```

```{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)
Expand All @@ -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 = "<br<\n")
if (all.equal(colnames(Betas_clean_quantile_bmiq_filtered_combated),rownames(PhenoData_clean))){
cat("The sample names between the corrected-filtered-normalized Beta dataset and the phenotype data were identical",
"ExpressionSet for filtered data was created and saved to the processed data folder", sep = "<br>\n")
if (all.equal(colnames(Betas_clean_unfiltered_quantile_bmiq_filtered_combated),rownames(PhenoData_clean))){
cat("The sample names between the unfiltered-normalized-corrected Beta dataset and the phenotype data were identical",
"ExpressionSet for unfiltered data was created and saved to the processed data folder", sep = "<br>\n")
} else {
cat("The sample names between the corrected-filtered-normalized Beta dataset and the phenotype data were not identical, please consult an expert", sep = "<br>\n")
cat("The sample names between the unfiltered-normalized-corrected Beta dataset and the phenotype data were not identical, please consult an expert", sep = "<br>\n")
}
cat("Summary table of preprocessing steps was updated", sep = "<br>\n")
Expand Down
Loading

0 comments on commit f63ed96

Please sign in to comment.