CPTACのデータとRを使ってがんや腫瘍のタンパク質発現量と患者の生存率の関係を解析する（CPTAC, v.4.4, 2024） 

# DEG
library(baySeq)
library(edgeR)
library(DESeq2)
library(limma)

# enrichment analysis
library(clusterProfiler)
library(fgsea)
library(ssGSEA2) # devtools::install_github("nicolerg/ssGSEA2") # devtools::install_github(broadinstitute/ssGSEA2.0) and devtools::install_github(broadinstitute/ssGSEA2) did not work.

library(escape)
library(GSVA)
library(pathview)

# immune cell deconvolution
library(ConsensusTME)
library(estimate)
library(immunedeconv)
library(MCPcounter)
library(msigdbr)
library(quantiseqr)
library(TIMER) #devtools::install_github('hanfeisun/TIMER'); it will not work!! use TIMER of immunodeconv if you need.
library(xCell)

# stats
library(NSM3) # sudo apt install libgmp-dev libgmp10 libgmp3-dev
library(pwr)
library(missForest)

# basic
library(beeswarm)
library(circlize)
library(org.Hs.eg.db)
library(org.Mm.eg.db)
library(doParallel)
library(foreach)
library(parallel) # Need it for use of baySeq
library(pheatmap)
library(rtracklayer)
library(rvest) # this is for extraction of html table; results of leading edge analysis.
library(xml2) # this is for extraction of html table; results of leading edge analysis.
library(openxlsx) # for read xlsx file.
library(ComplexHeatmap)
library(tidyverse)
library(data.table)
library(survminer)
library(survival)
library(gridExtra)

library(MASS)
library(coin)

library(AnnotationDbi)
library(org.Hs.eg.db)

library(RhpcBLASctl)

library(psych) # for geometric mean

# # shiny
# library(shiny)
# library(bslib)

# Check global environment
Sys.getenv()

# DEG
library(baySeq)
library(edgeR)
library(DESeq2)
library(limma)

# enrichment analysis
library(clusterProfiler)
library(fgsea)
library(ssGSEA2) # devtools::install_github("nicolerg/ssGSEA2") # devtools::install_github(broadinstitute/ssGSEA2.0) and devtools::install_github(broadinstitute/ssGSEA2) did not work.

library(escape)
library(GSVA)
library(pathview)

# immune cell deconvolution
library(ConsensusTME)
library(estimate)
library(immunedeconv)
library(MCPcounter)
library(msigdbr)
library(quantiseqr)
library(TIMER) #devtools::install_github('hanfeisun/TIMER'); it will not work!! use TIMER of immunodeconv if you need.
library(xCell)

# stats
library(NSM3) # sudo apt install libgmp-dev libgmp10 libgmp3-dev
library(pwr)
library(missForest)

# basic
library(beeswarm)
library(circlize)
library(org.Hs.eg.db)
library(org.Mm.eg.db)
library(doParallel)
library(foreach)
library(parallel) # Need it for use of baySeq
library(pheatmap)
library(rtracklayer)
library(rvest) # this is for extraction of html table; results of leading edge analysis.
library(xml2) # this is for extraction of html table; results of leading edge analysis.
library(openxlsx) # for read xlsx file.
library(ComplexHeatmap)
library(tidyverse)
library(data.table)
library(survminer)
library(survival)
library(gridExtra)

library(MASS)
library(coin)

library(AnnotationDbi)
library(org.Hs.eg.db)

library(RhpcBLASctl)

library(psych) # for geometric mean

# # shiny
# library(shiny)
# library(bslib)

# Check global environment
Sys.getenv()

# code 1

CPTAC <- list()
CPTAC <- list.files(path = "/mnt/seqdata/public_data/Blog/CPTAC/count", all.files = TRUE, full.names = TRUE, recursive = TRUE)

# code 1

CPTAC <- list()
CPTAC <- list.files(path = "/mnt/seqdata/public_data/Blog/CPTAC/count", all.files = TRUE, full.names = TRUE, recursive = TRUE)

# code 2

CPTAC_df <- unlist(CPTAC) %>% data.frame()
colnames(CPTAC_df) <- "path"

all_count_file <- CPTAC_df[str_detect(CPTAC_df$path, "tmt.*\\.tsv$")|str_detect(CPTAC_df$path, "itraq.tsv$"),] %>% data.frame()
colnames(all_count_file) <- "path"

proteome_and_peptide_count_file <- all_count_file[str_detect(all_count_file$path,"^(?!.*phospho).*"),] %>% data.frame()
colnames(proteome_and_peptide_count_file) <- "path"

proteome_count_file <- proteome_and_peptide_count_file[str_detect(proteome_and_peptide_count_file$path,"^(?!.*peptide).*"),] %>% data.frame()
colnames(proteome_count_file) <- "path"

# just check the file extention of iTRAQ data.
# itraq <- IRAKM_CPTAC_df[str_detect(IRAKM_CPTAC_df$path, "itraq.tsv$"),]



foreach(i = 1:nrow(proteome_count_file)) %do% {print(proteome_count_file[i,])}

# code 2

CPTAC_df <- unlist(CPTAC) %>% data.frame()
colnames(CPTAC_df) <- "path"

all_count_file <- CPTAC_df[str_detect(CPTAC_df$path, "tmt.*\\.tsv$")|str_detect(CPTAC_df$path, "itraq.tsv$"),] %>% data.frame()
colnames(all_count_file) <- "path"

proteome_and_peptide_count_file <- all_count_file[str_detect(all_count_file$path,"^(?!.*phospho).*"),] %>% data.frame()
colnames(proteome_and_peptide_count_file) <- "path"

proteome_count_file <- proteome_and_peptide_count_file[str_detect(proteome_and_peptide_count_file$path,"^(?!.*peptide).*"),] %>% data.frame()
colnames(proteome_count_file) <- "path"

# just check the file extention of iTRAQ data.
# itraq <- IRAKM_CPTAC_df[str_detect(IRAKM_CPTAC_df$path, "itraq.tsv$"),]



foreach(i = 1:nrow(proteome_count_file)) %do% {print(proteome_count_file[i,])}

# code 3

# Important
# Need to connect E drive.
# Read the dataset of PDC mass spec data for all patients, which are associated with "IRAK3".

proteome_collection <- list()

foreach(i = 1:nrow(proteome_count_file), .packages = 'readr', .combine = 'list') %do% {
  proteome_collection[[i]] <- read_tsv(proteome_count_file[i,])
}

# Need to check how many files are stored in proteome_collection.
# The data frame "file_number" is for just checking file number. It is not recqured for other analysis.
view(proteome_collection[[1]])
inner_join_for_check_file_number <- function(x, y) {inner_join(x, y, by = "Gene")}
file_number <- purrr::reduce(proteome_collection, inner_join_for_check_file_number)
ncol(file_number)

# code 3

# Important
# Need to connect E drive.
# Read the dataset of PDC mass spec data for all patients, which are associated with "IRAK3".

proteome_collection <- list()

foreach(i = 1:nrow(proteome_count_file), .packages = 'readr', .combine = 'list') %do% {
  proteome_collection[[i]] <- read_tsv(proteome_count_file[i,])
}

# Need to check how many files are stored in proteome_collection.
# The data frame "file_number" is for just checking file number. It is not recqured for other analysis.
view(proteome_collection[[1]])
inner_join_for_check_file_number <- function(x, y) {inner_join(x, y, by = "Gene")}
file_number <- purrr::reduce(proteome_collection, inner_join_for_check_file_number)
ncol(file_number)

# code 4

# Remove unnecessary columns from the dataset.

proteome_collection_wo_unnecessary_column <- list()

foreach(i = 1:(length(proteome_collection)), .packages = c('dplyr', 'purrr'), .combine = 'list') %do% {
  proteome_collection_wo_unnecessary_column[[i]] <- proteome_collection[[i]] %>% dplyr::select(-one_of("NCBIGeneID", "Authority", "Description", "Organism", "Chromosome", "Locus"))
}


length(proteome_collection_wo_unnecessary_column) # 53

PDC_study_ID <- str_extract(proteome_count_file$path, "PDC\\d+(?=/)")

length(PDC_study_ID) # 53
names(proteome_collection_wo_unnecessary_column) <- PDC_study_ID

for(i in 1:length(proteome_collection_wo_unnecessary_column)) {
temp_df <- colnames(proteome_collection_wo_unnecessary_column[[i]])[2:ncol(proteome_collection_wo_unnecessary_column[[i]])] %>% data.frame()
colnames(temp_df) <- "file_name"
temp_df$PDC_study_ID <- PDC_study_ID[i]
temp_df <- temp_df %>% unite("file_name_W_PDC_study_ID", c("PDC_study_ID", "file_name"), sep = "_", remove = FALSE)
colnames(proteome_collection_wo_unnecessary_column[[i]])[2:ncol(proteome_collection_wo_unnecessary_column[[i]])] <- temp_df$file_name_W_PDC_study_ID}

# code 4

# Remove unnecessary columns from the dataset.

proteome_collection_wo_unnecessary_column <- list()

foreach(i = 1:(length(proteome_collection)), .packages = c('dplyr', 'purrr'), .combine = 'list') %do% {
  proteome_collection_wo_unnecessary_column[[i]] <- proteome_collection[[i]] %>% dplyr::select(-one_of("NCBIGeneID", "Authority", "Description", "Organism", "Chromosome", "Locus"))
}


length(proteome_collection_wo_unnecessary_column) # 53

PDC_study_ID <- str_extract(proteome_count_file$path, "PDC\\d+(?=/)")

length(PDC_study_ID) # 53
names(proteome_collection_wo_unnecessary_column) <- PDC_study_ID

for(i in 1:length(proteome_collection_wo_unnecessary_column)) {
temp_df <- colnames(proteome_collection_wo_unnecessary_column[[i]])[2:ncol(proteome_collection_wo_unnecessary_column[[i]])] %>% data.frame()
colnames(temp_df) <- "file_name"
temp_df$PDC_study_ID <- PDC_study_ID[i]
temp_df <- temp_df %>% unite("file_name_W_PDC_study_ID", c("PDC_study_ID", "file_name"), sep = "_", remove = FALSE)
colnames(proteome_collection_wo_unnecessary_column[[i]])[2:ncol(proteome_collection_wo_unnecessary_column[[i]])] <- temp_df$file_name_W_PDC_study_ID}

# code 5

# Clean up the data frame Column name contain a lot of unnecessary characters, and therefore replace them with "_".

foreach(i = 1:(length(proteome_collection)), .packages = c('dplyr', 'purrr'), .combine = 'list') %do% {
colnames(proteome_collection_wo_unnecessary_column[[i]]) <- gsub("[- ]", "_", colnames(proteome_collection_wo_unnecessary_column[[i]]))
colnames(proteome_collection_wo_unnecessary_column[[i]]) <- gsub("[/]", "_", colnames(proteome_collection_wo_unnecessary_column[[i]]))
colnames(proteome_collection_wo_unnecessary_column[[i]]) <- gsub("[:]", "_", colnames(proteome_collection_wo_unnecessary_column[[i]]))
colnames(proteome_collection_wo_unnecessary_column[[i]]) <- gsub("[()]", "_", colnames(proteome_collection_wo_unnecessary_column[[i]]))
colnames(proteome_collection_wo_unnecessary_column[[i]]) <- gsub("[+]", "_", colnames(proteome_collection_wo_unnecessary_column[[i]]))
colnames(proteome_collection_wo_unnecessary_column[[i]]) <- gsub("__+", "_", colnames(proteome_collection_wo_unnecessary_column[[i]]))
colnames(proteome_collection_wo_unnecessary_column[[i]]) <- gsub("[.]", "", colnames(proteome_collection_wo_unnecessary_column[[i]]))
}

# code 5

# Clean up the data frame Column name contain a lot of unnecessary characters, and therefore replace them with "_".

foreach(i = 1:(length(proteome_collection)), .packages = c('dplyr', 'purrr'), .combine = 'list') %do% {
colnames(proteome_collection_wo_unnecessary_column[[i]]) <- gsub("[- ]", "_", colnames(proteome_collection_wo_unnecessary_column[[i]]))
colnames(proteome_collection_wo_unnecessary_column[[i]]) <- gsub("[/]", "_", colnames(proteome_collection_wo_unnecessary_column[[i]]))
colnames(proteome_collection_wo_unnecessary_column[[i]]) <- gsub("[:]", "_", colnames(proteome_collection_wo_unnecessary_column[[i]]))
colnames(proteome_collection_wo_unnecessary_column[[i]]) <- gsub("[()]", "_", colnames(proteome_collection_wo_unnecessary_column[[i]]))
colnames(proteome_collection_wo_unnecessary_column[[i]]) <- gsub("[+]", "_", colnames(proteome_collection_wo_unnecessary_column[[i]]))
colnames(proteome_collection_wo_unnecessary_column[[i]]) <- gsub("__+", "_", colnames(proteome_collection_wo_unnecessary_column[[i]]))
colnames(proteome_collection_wo_unnecessary_column[[i]]) <- gsub("[.]", "", colnames(proteome_collection_wo_unnecessary_column[[i]]))
}

# code 6

# Merge All dataset by "Gene".

full_join_by_gene <- function(x, y) {full_join(x, y, by = "Gene")}

merged <- reduce(proteome_collection_wo_unnecessary_column, full_join_by_gene)

ncol(merged)

check_colnames <- colnames(merged)[2:length(merged)] %>% data.frame()
colnames(check_colnames) <- "name"

# code 6

# Merge All dataset by "Gene".

full_join_by_gene <- function(x, y) {full_join(x, y, by = "Gene")}

merged <- reduce(proteome_collection_wo_unnecessary_column, full_join_by_gene)

ncol(merged)

check_colnames <- colnames(merged)[2:length(merged)] %>% data.frame()
colnames(check_colnames) <- "name"

# code 7

# Extract statistics (mean, median and sd for each library calculated by PDC portal)
merged_mean_median_sd <- merged[1:3,]

# Remove statistics from the dataset. 
merged_normalized_log2ratio <- merged[4:nrow(merged),]

# Add row name into statistics by PDC portal.
merged_mean_median_sd <- merged_mean_median_sd %>% column_to_rownames(var = "Gene")

# Add row name into the dataset.
merged_normalized_log2ratio <- merged_normalized_log2ratio %>% column_to_rownames(var = "Gene")

# code 7

# Extract statistics (mean, median and sd for each library calculated by PDC portal)
merged_mean_median_sd <- merged[1:3,]

# Remove statistics from the dataset. 
merged_normalized_log2ratio <- merged[4:nrow(merged),]

# Add row name into statistics by PDC portal.
merged_mean_median_sd <- merged_mean_median_sd %>% column_to_rownames(var = "Gene")

# Add row name into the dataset.
merged_normalized_log2ratio <- merged_normalized_log2ratio %>% column_to_rownames(var = "Gene")

# code 8

## Discrepancies in summary statistics at this url (https://pdc.cancer.gov/pdc-docs/data-analysis-guides#data-processing-and-harmonization) described CPTAC is using median normalization, and the values in data frame "merged_mean_median_sd" are a value before median normalization. 

## The description is "The summary statistics (mean, median, standard deviation) in quantitation files are computed before median normalization. To match these values, add the median value back to each log2ratio column prior to performing statistical calculations."

## frequent question (https://proteomic.datacommons.cancer.gov/pdc/faq) discribed how to calculate back normalization. The discription is "Why don't the summary statistics (Mean, Median, StdDev) in the quantitation files match my calculations? The summary statistics in the quantitation files are computed before median normalization of the log2ratio values. To match these statistics, you need to add the median back to each column."

# these are the median and mean value by CPTAC
merged_mean_median_sd[2,1] # median = -1.024735
merged_mean_median_sd[1,1] # mean = -1.030376

# Each value without NA. 
normalized_log2ratio <- merged_normalized_log2ratio[,1]
before_normalized_log2ratio <- normalized_log2ratio + merged_mean_median_sd[2,1]
before_normalized_log2ratio_wo_na <- before_normalized_log2ratio[!is.na(before_normalized_log2ratio)]

# Check distribution before substituting NA with 0
normalized_log2ratio_w_NA <- normalized_log2ratio
hist(normalized_log2ratio, breaks = 100, ylim = c(0, 14000))

merged_mean_median_sd[2,1] # -1.024735
median(before_normalized_log2ratio_wo_na) # -1.029058
mean(before_normalized_log2ratio_wo_na) # -1.039549

# NA replaced with 0
normalized_log2ratio[is.na(normalized_log2ratio)] <- 0
median(normalized_log2ratio) # 0
mean(normalized_log2ratio) # -0.005616025
before_normalized_log2ratio_NA_is_0 <- normalized_log2ratio + merged_mean_median_sd[2,1]
2^median(before_normalized_log2ratio_NA_is_0) # 0.4915004

# distribution if NA were sustituted with 0. The distribition was totally weird. 
hist(normalized_log2ratio, breaks = 100, ylim = c(0, 14000))

merged_mean_median_sd[2,1] # -1.024735
median(before_normalized_log2ratio_NA_is_0) # -1.024735 # It is consistent with CPTAC calculation.
mean(before_normalized_log2ratio_NA_is_0) # -1.030351
# exp(mean(log(before_normalized_log2ratio_NA_is_0))) # NaN
geometric.mean(before_normalized_log2ratio_NA_is_0) # 0.2752213

## When 0 was substituted into NA in median normalized intensity by CPTAC or submitters, the median value was same as CPTAC's median. Following order of calculation was able to reproduce the median value by CPTAC or submitters.

# 1. Replace NA with 0 (Because the value should be centered with median value.)
# 2. Plus median (turned the value before normalization.)
# 3. Then median will be same as median by CPTAC. The result is log2ratio before normalization.
# 4. calculate power of 2. The result is ratio before normalization.
# 5. the value (intensity) before normalization  

## If 0 was substituted into NA in median-normalized intensity, the data set provided by CPTAC, it means that the value was centered with median value. 


# So, What meaning of NA in original dataset is very important. Why is NA generated? In mass spec data, NA means "Not detected" or "the value is under detection limit". The real value is unknown. In addition, each column has different median. what should be.

# code 8

## Discrepancies in summary statistics at this url (https://pdc.cancer.gov/pdc-docs/data-analysis-guides#data-processing-and-harmonization) described CPTAC is using median normalization, and the values in data frame "merged_mean_median_sd" are a value before median normalization. 

## The description is "The summary statistics (mean, median, standard deviation) in quantitation files are computed before median normalization. To match these values, add the median value back to each log2ratio column prior to performing statistical calculations."

## frequent question (https://proteomic.datacommons.cancer.gov/pdc/faq) discribed how to calculate back normalization. The discription is "Why don't the summary statistics (Mean, Median, StdDev) in the quantitation files match my calculations? The summary statistics in the quantitation files are computed before median normalization of the log2ratio values. To match these statistics, you need to add the median back to each column."

# these are the median and mean value by CPTAC
merged_mean_median_sd[2,1] # median = -1.024735
merged_mean_median_sd[1,1] # mean = -1.030376

# Each value without NA. 
normalized_log2ratio <- merged_normalized_log2ratio[,1]
before_normalized_log2ratio <- normalized_log2ratio + merged_mean_median_sd[2,1]
before_normalized_log2ratio_wo_na <- before_normalized_log2ratio[!is.na(before_normalized_log2ratio)]

# Check distribution before substituting NA with 0
normalized_log2ratio_w_NA <- normalized_log2ratio
hist(normalized_log2ratio, breaks = 100, ylim = c(0, 14000))

merged_mean_median_sd[2,1] # -1.024735
median(before_normalized_log2ratio_wo_na) # -1.029058
mean(before_normalized_log2ratio_wo_na) # -1.039549

# NA replaced with 0
normalized_log2ratio[is.na(normalized_log2ratio)] <- 0
median(normalized_log2ratio) # 0
mean(normalized_log2ratio) # -0.005616025
before_normalized_log2ratio_NA_is_0 <- normalized_log2ratio + merged_mean_median_sd[2,1]
2^median(before_normalized_log2ratio_NA_is_0) # 0.4915004

# distribution if NA were sustituted with 0. The distribition was totally weird. 
hist(normalized_log2ratio, breaks = 100, ylim = c(0, 14000))

merged_mean_median_sd[2,1] # -1.024735
median(before_normalized_log2ratio_NA_is_0) # -1.024735 # It is consistent with CPTAC calculation.
mean(before_normalized_log2ratio_NA_is_0) # -1.030351
# exp(mean(log(before_normalized_log2ratio_NA_is_0))) # NaN
geometric.mean(before_normalized_log2ratio_NA_is_0) # 0.2752213

## When 0 was substituted into NA in median normalized intensity by CPTAC or submitters, the median value was same as CPTAC's median. Following order of calculation was able to reproduce the median value by CPTAC or submitters.

# 1. Replace NA with 0 (Because the value should be centered with median value.)
# 2. Plus median (turned the value before normalization.)
# 3. Then median will be same as median by CPTAC. The result is log2ratio before normalization.
# 4. calculate power of 2. The result is ratio before normalization.
# 5. the value (intensity) before normalization  

## If 0 was substituted into NA in median-normalized intensity, the data set provided by CPTAC, it means that the value was centered with median value. 


# So, What meaning of NA in original dataset is very important. Why is NA generated? In mass spec data, NA means "Not detected" or "the value is under detection limit". The real value is unknown. In addition, each column has different median. what should be.

# code 9

# Following is stats of median value used for normalization at CPTAC or submitters.
merged_mean_median_sd_shared <- merged_mean_median_sd[,str_detect(colnames(merged_mean_median_sd), "Unshared_") == FALSE]
merged_mean_median_sd_shared_vector <- merged_mean_median_sd_shared[2,] %>% unlist()

hist(merged_mean_median_sd_shared_vector, 
     breaks = 1000,
     main = "median value used for normalization at CPTAC or submitters",
     xlab = "median",
     xlim = c(-3.0, 1.0))

quantile(merged_mean_median_sd_shared_vector, probs = c(0.25, 0.5, 0.75))

# code 9

# Following is stats of median value used for normalization at CPTAC or submitters.
merged_mean_median_sd_shared <- merged_mean_median_sd[,str_detect(colnames(merged_mean_median_sd), "Unshared_") == FALSE]
merged_mean_median_sd_shared_vector <- merged_mean_median_sd_shared[2,] %>% unlist()

hist(merged_mean_median_sd_shared_vector, 
     breaks = 1000,
     main = "median value used for normalization at CPTAC or submitters",
     xlab = "median",
     xlim = c(-3.0, 1.0))

quantile(merged_mean_median_sd_shared_vector, probs = c(0.25, 0.5, 0.75))

# code 10

# Calculate percent of NA in a row, meaning decide which genes were how many NAs are conteined. 

# MARGIN = ; a vector giving the subscripts which the function will be applied over. E.g., for a matrix 1 indicates rows, 2 indicates columns, c(1, 2) indicates rows and columns. Where X has named dimnames, it can be a character vector selecting dimension names.

percent_na <- apply(merged_normalized_log2ratio, 1, function(x) {sum(is.na(x))/ncol(merged_normalized_log2ratio)})

merged_normalized_log2ratio_w_percent_na <- cbind(merged_normalized_log2ratio, percent_na)
merged_normalized_log2ratio_w_percent_na$percent_na

nrow(merged_normalized_log2ratio_w_percent_na) # 20414 genes
ncol(merged_normalized_log2ratio_w_percent_na) # 11435 = 11434 cases + percent_na 

# code 10

# Calculate percent of NA in a row, meaning decide which genes were how many NAs are conteined. 

# MARGIN = ; a vector giving the subscripts which the function will be applied over. E.g., for a matrix 1 indicates rows, 2 indicates columns, c(1, 2) indicates rows and columns. Where X has named dimnames, it can be a character vector selecting dimension names.

percent_na <- apply(merged_normalized_log2ratio, 1, function(x) {sum(is.na(x))/ncol(merged_normalized_log2ratio)})

merged_normalized_log2ratio_w_percent_na <- cbind(merged_normalized_log2ratio, percent_na)
merged_normalized_log2ratio_w_percent_na$percent_na

nrow(merged_normalized_log2ratio_w_percent_na) # 20414 genes
ncol(merged_normalized_log2ratio_w_percent_na) # 11435 = 11434 cases + percent_na 

# code 11

# Distribution of percent of NAs. 
hist(merged_normalized_log2ratio_w_percent_na$percent_na,
     breaks = 10000,
     main = "Histogram of NA ratio",
     xlab = "NA ratio",
     xlim = c(0, 1.00))


# Distribution of percent of NAs. 
hist(merged_normalized_log2ratio_w_percent_na$percent_na,
     breaks = 10000,
     main = "Histogram of NA ratio",
     xlab = "NA ratio",
     xlim = c(0.85, 1.00))

# Distribution of percent of NAs. 
hist(merged_normalized_log2ratio_w_percent_na$percent_na,
     breaks = 10000,
     main = "Histogram of NA ratio",
     xlab = "NA ratio",
     xlim = c(0.0, 0.2))

# code 11

# Distribution of percent of NAs. 
hist(merged_normalized_log2ratio_w_percent_na$percent_na,
     breaks = 10000,
     main = "Histogram of NA ratio",
     xlab = "NA ratio",
     xlim = c(0, 1.00))


# Distribution of percent of NAs. 
hist(merged_normalized_log2ratio_w_percent_na$percent_na,
     breaks = 10000,
     main = "Histogram of NA ratio",
     xlab = "NA ratio",
     xlim = c(0.85, 1.00))

# Distribution of percent of NAs. 
hist(merged_normalized_log2ratio_w_percent_na$percent_na,
     breaks = 10000,
     main = "Histogram of NA ratio",
     xlab = "NA ratio",
     xlim = c(0.0, 0.2))

# code 12

# If a gene was detected in more than 90 % of patients, a gene was including in analysis (alternatively, a gene was detected in only 10 % of patients, a gene was omitted from analysis.)

merged_normalized_log2ratio_2 <- merged_normalized_log2ratio_w_percent_na %>% filter(percent_na < 0.1)

nrow(merged_normalized_log2ratio_2) # 5208 genes were remained.
ncol(merged_normalized_log2ratio_2) # 11435 = 11434 cases + percent_na 

# Add "Gene" into column name.
merged_normalized_log2ratio_2 <- merged_normalized_log2ratio_2 %>% rownames_to_column(var = "Gene")

# Pick up genes with NA ratio > 0.1. 

gene_to_be_omitted <- merged_normalized_log2ratio_w_percent_na[merged_normalized_log2ratio_w_percent_na$percent_na >= 0.1,]$Gene

length(gene_to_be_omitted) # 15206 genes

# code 12

# If a gene was detected in more than 90 % of patients, a gene was including in analysis (alternatively, a gene was detected in only 10 % of patients, a gene was omitted from analysis.)

merged_normalized_log2ratio_2 <- merged_normalized_log2ratio_w_percent_na %>% filter(percent_na < 0.1)

nrow(merged_normalized_log2ratio_2) # 5208 genes were remained.
ncol(merged_normalized_log2ratio_2) # 11435 = 11434 cases + percent_na 

# Add "Gene" into column name.
merged_normalized_log2ratio_2 <- merged_normalized_log2ratio_2 %>% rownames_to_column(var = "Gene")

# Pick up genes with NA ratio > 0.1. 

gene_to_be_omitted <- merged_normalized_log2ratio_w_percent_na[merged_normalized_log2ratio_w_percent_na$percent_na >= 0.1,]$Gene

length(gene_to_be_omitted) # 15206 genes

# code 13

# Transform the dataset to tibble.

merged_normalized_log2ratio_3 <- merged_normalized_log2ratio_2 %>% as_tibble()

# Remove unnecessary columns.

merged_normalized_log2ratio_3 <- merged_normalized_log2ratio_3 %>% 
  select(-one_of("percent_na"))

# Set column name as "Gene".
merged_normalized_log2ratio_3 <- merged_normalized_log2ratio_3 %>% column_to_rownames(var = "Gene")

ncol(merged_normalized_log2ratio_3) # 11434 cases
nrow(merged_normalized_log2ratio_3) # 5208 genes

# code 13

# Transform the dataset to tibble.

merged_normalized_log2ratio_3 <- merged_normalized_log2ratio_2 %>% as_tibble()

# Remove unnecessary columns.

merged_normalized_log2ratio_3 <- merged_normalized_log2ratio_3 %>% 
  select(-one_of("percent_na"))

# Set column name as "Gene".
merged_normalized_log2ratio_3 <- merged_normalized_log2ratio_3 %>% column_to_rownames(var = "Gene")

ncol(merged_normalized_log2ratio_3) # 11434 cases
nrow(merged_normalized_log2ratio_3) # 5208 genes

# code 14

# This collection of the files contain QC files, which were probably used for calculating ratio of proteins for each sample at other laboratory. Therefore they will be omitted from the dataset.

merged_normalized_log2ratio_4 <-merged_normalized_log2ratio_3 %>% select(
  colnames(merged_normalized_log2ratio_3)[!grepl(pattern = "^(QC|Withdrawn)", colnames(merged_normalized_log2ratio_3), ignore.case = FALSE)]
)

merged_normalized_log2ratio_5 <-merged_normalized_log2ratio_4 %>% select(
  colnames(merged_normalized_log2ratio_4)[!grepl(pattern = "POOL|CONTROL", colnames(merged_normalized_log2ratio_4), ignore.case = TRUE)]
)

# merged_normalized_log2ratio_6 <-merged_normalized_log2ratio_5 %>% select(
#   colnames(merged_normalized_log2ratio_5)[!grepl(pattern = "normal", colnames(merged_normalized_log2ratio_5), ignore.case = TRUE)]
# )

merged_normalized_log2ratio_6 <-merged_normalized_log2ratio_5 %>% select(
  colnames(merged_normalized_log2ratio_5)[!grepl(pattern = "Unshared", colnames(merged_normalized_log2ratio_5), ignore.case = TRUE)]
)

ncol(merged_normalized_log2ratio_6) # 5144

file.remove("/mnt/team4tb/Dropbox/Blog/2026 01 03 blog CPTAC/colnames_merged_normalized_log2ratio_6.txt")
cat(colnames(merged_normalized_log2ratio_6), sep = "\n", file = "/mnt/team4tb/Dropbox/Blog/2026 01 03 blog CPTAC/colnames_merged_normalized_log2ratio_6.txt")

# 5144 cases were remained, and they will be used for further analysis.

# code 14

# This collection of the files contain QC files, which were probably used for calculating ratio of proteins for each sample at other laboratory. Therefore they will be omitted from the dataset.

merged_normalized_log2ratio_4 <-merged_normalized_log2ratio_3 %>% select(
  colnames(merged_normalized_log2ratio_3)[!grepl(pattern = "^(QC|Withdrawn)", colnames(merged_normalized_log2ratio_3), ignore.case = FALSE)]
)

merged_normalized_log2ratio_5 <-merged_normalized_log2ratio_4 %>% select(
  colnames(merged_normalized_log2ratio_4)[!grepl(pattern = "POOL|CONTROL", colnames(merged_normalized_log2ratio_4), ignore.case = TRUE)]
)

# merged_normalized_log2ratio_6 <-merged_normalized_log2ratio_5 %>% select(
#   colnames(merged_normalized_log2ratio_5)[!grepl(pattern = "normal", colnames(merged_normalized_log2ratio_5), ignore.case = TRUE)]
# )

merged_normalized_log2ratio_6 <-merged_normalized_log2ratio_5 %>% select(
  colnames(merged_normalized_log2ratio_5)[!grepl(pattern = "Unshared", colnames(merged_normalized_log2ratio_5), ignore.case = TRUE)]
)

ncol(merged_normalized_log2ratio_6) # 5144

file.remove("/mnt/team4tb/Dropbox/Blog/2026 01 03 blog CPTAC/colnames_merged_normalized_log2ratio_6.txt")
cat(colnames(merged_normalized_log2ratio_6), sep = "\n", file = "/mnt/team4tb/Dropbox/Blog/2026 01 03 blog CPTAC/colnames_merged_normalized_log2ratio_6.txt")

# 5144 cases were remained, and they will be used for further analysis.

# code 15

## Transform data frame "merged_normalized_log2ratio_6" into matrix "merged_normalized_log2ratio_6_mat"
merged_normalized_log2ratio_6_mat <- as.matrix(merged_normalized_log2ratio_6) 
view(merged_normalized_log2ratio_6_mat[1:100, 1:100]) # check the data.

## Retrieve the cases in data frame "merged_mean_median_sd" that are same as the cases in the matrix "merged_normalized_log2ratio_6_mat".
merged_mean_median_sd_2 <-merged_mean_median_sd[,colnames(merged_normalized_log2ratio_6_mat)]

## Transform the data frame "merged_mean_median_sd_2" into matrix "merged_mean_median_sd_2_mat"
merged_mean_median_sd_2_mat <- as.matrix(merged_mean_median_sd_2)
view(merged_mean_median_sd_2_mat[, 1:100]) # check the data.

## the matrix "merged_normalized_log2ratio_6_mat" was already centered with median value. To get back to the value before normalization, the median value was added into each column. The matrix "merged_log2ratio_6_mat" was the data before normalization, but it is still log2 transformed value.
merged_log2ratio_6_mat <- sweep(x = merged_normalized_log2ratio_6_mat,
                                            MARGIN = 2,
                                            STATS = merged_mean_median_sd_2_mat[2,],
                                            FUN = "+") # The matrix "merged_log2ratio_6_mat" is normalized_log2ratio + median = log2ratio

## the matrix "merged_analysis_before_normalization_ratio" was argument of the data.
merged_analysis_before_normalization_ratio <- 2^merged_log2ratio_6_mat # the matrix "merged_analysis_before_normalization_ratio" is raw ratio before normalization. If the value gets back to normalized ratio, following equation should be used.

# log2(2^merged_log2ratio_6_mat) - median = normalized_log2ratio

view(merged_analysis_before_normalization_ratio[1:100, 1:100]) # check the data.

## make sure the value of the matrix "merged_analysis_before_normalization_ratio" was numeric.
merged_analysis_before_normalization_ratio <- apply(merged_analysis_before_normalization_ratio, 2, as.numeric)

## make sure the rownames of the matrix "merged_analysis_before_normalization_ratio" is same as that of the matrix "merged_log2ratio_6_mat", which is original dataset. 
rownames(merged_analysis_before_normalization_ratio) <- rownames(merged_log2ratio_6_mat)

view(merged_analysis_before_normalization_ratio[1:100, 1:100]) # check the data.

## Get the matrix "merged_analysis_before_normalization_ratio" back to data frame "merged_analysis_before_normalization_ratio_df"
merged_analysis_before_normalization_ratio_df <- data.frame(merged_analysis_before_normalization_ratio)
colnames(merged_analysis_before_normalization_ratio_df) <- colnames(merged_normalized_log2ratio_6)
# 
# merged_analysis_before_normalization_ratio_df[is.na(merged_analysis_before_normalization_ratio_df)] <- 0
# 
# options(scipen = 100)

# code 15

## Transform data frame "merged_normalized_log2ratio_6" into matrix "merged_normalized_log2ratio_6_mat"
merged_normalized_log2ratio_6_mat <- as.matrix(merged_normalized_log2ratio_6) 
view(merged_normalized_log2ratio_6_mat[1:100, 1:100]) # check the data.

## Retrieve the cases in data frame "merged_mean_median_sd" that are same as the cases in the matrix "merged_normalized_log2ratio_6_mat".
merged_mean_median_sd_2 <-merged_mean_median_sd[,colnames(merged_normalized_log2ratio_6_mat)]

## Transform the data frame "merged_mean_median_sd_2" into matrix "merged_mean_median_sd_2_mat"
merged_mean_median_sd_2_mat <- as.matrix(merged_mean_median_sd_2)
view(merged_mean_median_sd_2_mat[, 1:100]) # check the data.

## the matrix "merged_normalized_log2ratio_6_mat" was already centered with median value. To get back to the value before normalization, the median value was added into each column. The matrix "merged_log2ratio_6_mat" was the data before normalization, but it is still log2 transformed value.
merged_log2ratio_6_mat <- sweep(x = merged_normalized_log2ratio_6_mat,
                                            MARGIN = 2,
                                            STATS = merged_mean_median_sd_2_mat[2,],
                                            FUN = "+") # The matrix "merged_log2ratio_6_mat" is normalized_log2ratio + median = log2ratio

## the matrix "merged_analysis_before_normalization_ratio" was argument of the data.
merged_analysis_before_normalization_ratio <- 2^merged_log2ratio_6_mat # the matrix "merged_analysis_before_normalization_ratio" is raw ratio before normalization. If the value gets back to normalized ratio, following equation should be used.

# log2(2^merged_log2ratio_6_mat) - median = normalized_log2ratio

view(merged_analysis_before_normalization_ratio[1:100, 1:100]) # check the data.

## make sure the value of the matrix "merged_analysis_before_normalization_ratio" was numeric.
merged_analysis_before_normalization_ratio <- apply(merged_analysis_before_normalization_ratio, 2, as.numeric)

## make sure the rownames of the matrix "merged_analysis_before_normalization_ratio" is same as that of the matrix "merged_log2ratio_6_mat", which is original dataset. 
rownames(merged_analysis_before_normalization_ratio) <- rownames(merged_log2ratio_6_mat)

view(merged_analysis_before_normalization_ratio[1:100, 1:100]) # check the data.

## Get the matrix "merged_analysis_before_normalization_ratio" back to data frame "merged_analysis_before_normalization_ratio_df"
merged_analysis_before_normalization_ratio_df <- data.frame(merged_analysis_before_normalization_ratio)
colnames(merged_analysis_before_normalization_ratio_df) <- colnames(merged_normalized_log2ratio_6)
# 
# merged_analysis_before_normalization_ratio_df[is.na(merged_analysis_before_normalization_ratio_df)] <- 0
# 
# options(scipen = 100)

# code 16

# "merged_analysis_before_normalization_ratio_df" is the dataset that was put the normalized dataset back to the (possiblly) original dataset. 

# "samples_to_be_analyzed" at following line was column names, which are file names of tmt10, in dataset "merged_analysis_before_normalization_ratio_df". 

# Files in "samples_to_be_analyzed" will be used for further analysis, such as calculation of immune suppression score. 
view(colnames(merged_analysis_before_normalization_ratio_df))

samples_to_be_analyzed <- colnames(merged_analysis_before_normalization_ratio_df) 

# In case of file will be saved on C drive.

file.remove("/mnt/team4tb/Dropbox/Blog/2026 01 03 blog CPTAC/samples_to_be_analyzed.txt")
cat(samples_to_be_analyzed, sep = "\n", file = "/mnt/team4tb/Dropbox/Blog/2026 01 03 blog CPTAC/samples_to_be_analyzed.txt")

# code 16

# "merged_analysis_before_normalization_ratio_df" is the dataset that was put the normalized dataset back to the (possiblly) original dataset. 

# "samples_to_be_analyzed" at following line was column names, which are file names of tmt10, in dataset "merged_analysis_before_normalization_ratio_df". 

# Files in "samples_to_be_analyzed" will be used for further analysis, such as calculation of immune suppression score. 
view(colnames(merged_analysis_before_normalization_ratio_df))

samples_to_be_analyzed <- colnames(merged_analysis_before_normalization_ratio_df) 

# In case of file will be saved on C drive.

file.remove("/mnt/team4tb/Dropbox/Blog/2026 01 03 blog CPTAC/samples_to_be_analyzed.txt")
cat(samples_to_be_analyzed, sep = "\n", file = "/mnt/team4tb/Dropbox/Blog/2026 01 03 blog CPTAC/samples_to_be_analyzed.txt")

# code 17

## the data frame "merged_normalized_log2ratio_6" was the dataframe without unnecessary columns. This data will be used for further analysis. Here, column names (which are cases) and column number (which is number of cases) will be checked by colnames(), ncol(), and length(). 

colnames(merged_normalized_log2ratio_6)
ncol(merged_normalized_log2ratio_6) # 5144
samples_to_be_analyzed
length(samples_to_be_analyzed) # 5144

## Check if the case number between dataframe "merged_normalized_log2ratio_6" and "samples_to_be_analyzed" is same or not.

table(!colnames(merged_normalized_log2ratio_6) %in% samples_to_be_analyzed) # FALSE 5144 means that all column if the data frame "merged_normalized_log2ratio_6" is same as the vectior "samples_to_be_analyzed".

check_match <- colnames(merged_normalized_log2ratio_6)[colnames(merged_normalized_log2ratio_6) %in% samples_to_be_analyzed] 

check_unmatch <- colnames(merged_normalized_log2ratio_6)[!colnames(merged_normalized_log2ratio_6) %in% samples_to_be_analyzed] 

length(check_match) # 5144
length(check_unmatch) # 0
length(check_match) + length(check_unmatch) # 5144
view(colnames(merged_normalized_log2ratio_6))
view(samples_to_be_analyzed)
view(check_unmatch)

## the case number between dataframe "merged_normalized_log2ratio_6" and "samples_to_be_analyzed" is same.

# code 17

## the data frame "merged_normalized_log2ratio_6" was the dataframe without unnecessary columns. This data will be used for further analysis. Here, column names (which are cases) and column number (which is number of cases) will be checked by colnames(), ncol(), and length(). 

colnames(merged_normalized_log2ratio_6)
ncol(merged_normalized_log2ratio_6) # 5144
samples_to_be_analyzed
length(samples_to_be_analyzed) # 5144

## Check if the case number between dataframe "merged_normalized_log2ratio_6" and "samples_to_be_analyzed" is same or not.

table(!colnames(merged_normalized_log2ratio_6) %in% samples_to_be_analyzed) # FALSE 5144 means that all column if the data frame "merged_normalized_log2ratio_6" is same as the vectior "samples_to_be_analyzed".

check_match <- colnames(merged_normalized_log2ratio_6)[colnames(merged_normalized_log2ratio_6) %in% samples_to_be_analyzed] 

check_unmatch <- colnames(merged_normalized_log2ratio_6)[!colnames(merged_normalized_log2ratio_6) %in% samples_to_be_analyzed] 

length(check_match) # 5144
length(check_unmatch) # 0
length(check_match) + length(check_unmatch) # 5144
view(colnames(merged_normalized_log2ratio_6))
view(samples_to_be_analyzed)
view(check_unmatch)

## the case number between dataframe "merged_normalized_log2ratio_6" and "samples_to_be_analyzed" is same.

# code 18

## prepare the dataset that is argument before normalization.
cases_pre_argu_normalization <- colnames(merged_analysis_before_normalization_ratio_df)
proteins_pre_argu_normalization <- rownames(merged_analysis_before_normalization_ratio_df)

length(cases_pre_argu_normalization) # 5144
length(proteins_pre_argu_normalization) # 5208

merged_analysis_before_normalization_ratio_mat <- as.matrix(merged_analysis_before_normalization_ratio_df)

merged_analysis_before_normalization_ratio_mat_t <- t(merged_analysis_before_normalization_ratio_mat)
view(merged_analysis_before_normalization_ratio_mat_t[1:1000, 1:10]) # check the data.

## prepare the dataset before normalization, it is still logarithm. 

cases_pre_normalization_log <- colnames(merged_log2ratio_6_mat)
proteins_pre_normalization_log <- rownames(merged_log2ratio_6_mat)

length(cases_pre_normalization_log) # 5144
length(proteins_pre_normalization_log) # 5208

merged_log2ratio_6_mat_t <- t(merged_log2ratio_6_mat)
view(merged_log2ratio_6_mat_t[1:1000, 1:10]) # check the data.

## prepare the dataset that is already normalized by CPTAC or submitters.

cases_post_normalization <- colnames(merged_normalized_log2ratio_6_mat)
proteins_post_normalization <- rownames(merged_normalized_log2ratio_6_mat)

length(cases_post_normalization) # 5144
length(proteins_post_normalization) # 5208

merged_normalized_log2ratio_6_mat_t <- t(merged_normalized_log2ratio_6_mat)
view(merged_normalized_log2ratio_6_mat_t[1:1000, 1:10]) # check the data.

## check cases and proteins

table(cases_pre_argu_normalization %in% cases_pre_normalization_log) # TRUE 5144
table(proteins_pre_argu_normalization %in% proteins_pre_normalization_log) # TRUE 5208

table(cases_pre_argu_normalization %in% cases_post_normalization) # TRUE 5144
table(proteins_pre_argu_normalization %in% proteins_post_normalization) # TRUE 5208

table(cases_pre_normalization_log %in% cases_post_normalization) # TRUE 5144
table(proteins_pre_normalization_log %in% proteins_post_normalization) # TRUE 5208

# code 18

## prepare the dataset that is argument before normalization.
cases_pre_argu_normalization <- colnames(merged_analysis_before_normalization_ratio_df)
proteins_pre_argu_normalization <- rownames(merged_analysis_before_normalization_ratio_df)

length(cases_pre_argu_normalization) # 5144
length(proteins_pre_argu_normalization) # 5208

merged_analysis_before_normalization_ratio_mat <- as.matrix(merged_analysis_before_normalization_ratio_df)

merged_analysis_before_normalization_ratio_mat_t <- t(merged_analysis_before_normalization_ratio_mat)
view(merged_analysis_before_normalization_ratio_mat_t[1:1000, 1:10]) # check the data.

## prepare the dataset before normalization, it is still logarithm. 

cases_pre_normalization_log <- colnames(merged_log2ratio_6_mat)
proteins_pre_normalization_log <- rownames(merged_log2ratio_6_mat)

length(cases_pre_normalization_log) # 5144
length(proteins_pre_normalization_log) # 5208

merged_log2ratio_6_mat_t <- t(merged_log2ratio_6_mat)
view(merged_log2ratio_6_mat_t[1:1000, 1:10]) # check the data.

## prepare the dataset that is already normalized by CPTAC or submitters.

cases_post_normalization <- colnames(merged_normalized_log2ratio_6_mat)
proteins_post_normalization <- rownames(merged_normalized_log2ratio_6_mat)

length(cases_post_normalization) # 5144
length(proteins_post_normalization) # 5208

merged_normalized_log2ratio_6_mat_t <- t(merged_normalized_log2ratio_6_mat)
view(merged_normalized_log2ratio_6_mat_t[1:1000, 1:10]) # check the data.

## check cases and proteins

table(cases_pre_argu_normalization %in% cases_pre_normalization_log) # TRUE 5144
table(proteins_pre_argu_normalization %in% proteins_pre_normalization_log) # TRUE 5208

table(cases_pre_argu_normalization %in% cases_post_normalization) # TRUE 5144
table(proteins_pre_argu_normalization %in% proteins_post_normalization) # TRUE 5208

table(cases_pre_normalization_log %in% cases_post_normalization) # TRUE 5144
table(proteins_pre_normalization_log %in% proteins_post_normalization) # TRUE 5208

# code 19

# save.image("/mnt/team4tb/Dropbox/Blog/2026 01 03 blog CPTAC/2026 01 10 CPTAC expression matrix_1.RData")

# code 19

# save.image("/mnt/team4tb/Dropbox/Blog/2026 01 03 blog CPTAC/2026 01 10 CPTAC expression matrix_1.RData")

# code 20

# Impute missing values using missForest for dataset BEFORE normalization
cl <- makeCluster(detectCores())
registerDoParallel(cl)
start_time_pre_argu <- Sys.time()
set.seed(20260121) # for reproducibility
result_missForest_pre_argu <- missForest(merged_analysis_before_normalization_ratio_mat_t, parallelize = "variables", verbose = TRUE)
end_time_pre_argu <- Sys.time()
stopCluster(cl)

gc(T,T,T)
gc(T,T,T)

spendintg_time_pre_argu <- end_time_pre_argu - start_time_pre_argu
spendintg_time_pre_argu # Time difference of 20.40168 hours

# 15.24111 hours will be spent for missForest with NH-U12A.

# code 20

# Impute missing values using missForest for dataset BEFORE normalization
cl <- makeCluster(detectCores())
registerDoParallel(cl)
start_time_pre_argu <- Sys.time()
set.seed(20260121) # for reproducibility
result_missForest_pre_argu <- missForest(merged_analysis_before_normalization_ratio_mat_t, parallelize = "variables", verbose = TRUE)
end_time_pre_argu <- Sys.time()
stopCluster(cl)

gc(T,T,T)
gc(T,T,T)

spendintg_time_pre_argu <- end_time_pre_argu - start_time_pre_argu
spendintg_time_pre_argu # Time difference of 20.40168 hours

# 15.24111 hours will be spent for missForest with NH-U12A.

# code 21
# Impute missing values using missForest for dataset AFTER normalization
# The matrix "merged_normalized_log2ratio_6_mat_t" is from the matrix "merged_normalized_log2ratio_6_mat".
# The matrix "merged_normalized_log2ratio_6_mat" is normalized_log2ratio (median-normalized log2ratio) + median = log2ratio. log2ratio is the value before log2 transform.
cl <- makeCluster(detectCores())
registerDoParallel(cl)
start_time_pre <- Sys.time()
set.seed(20260121) # for reproducibility
result_missForest_pre <- missForest(merged_log2ratio_6_mat_t, parallelize = "variables", verbose = TRUE)
end_time_pre <- Sys.time()
stopCluster(cl)

gc(T,T,T)
gc(T,T,T)

spendintg_time_pre <- end_time_pre - start_time_pre # Time difference of 1.196548 days
spendintg_time_pre # Time difference of 1.217632 days

# code 21
# Impute missing values using missForest for dataset AFTER normalization
# The matrix "merged_normalized_log2ratio_6_mat_t" is from the matrix "merged_normalized_log2ratio_6_mat".
# The matrix "merged_normalized_log2ratio_6_mat" is normalized_log2ratio (median-normalized log2ratio) + median = log2ratio. log2ratio is the value before log2 transform.
cl <- makeCluster(detectCores())
registerDoParallel(cl)
start_time_pre <- Sys.time()
set.seed(20260121) # for reproducibility
result_missForest_pre <- missForest(merged_log2ratio_6_mat_t, parallelize = "variables", verbose = TRUE)
end_time_pre <- Sys.time()
stopCluster(cl)

gc(T,T,T)
gc(T,T,T)

spendintg_time_pre <- end_time_pre - start_time_pre # Time difference of 1.196548 days
spendintg_time_pre # Time difference of 1.217632 days

# code 22
# Impute missing values using missForest for dataset AFTER normalization
# The matrix "merged_normalized_log2ratio_6_mat_t" is from the matrix "merged_normalized_log2ratio_6_mat".
# The matrix "merged_normalized_log2ratio_6_mat" is normalized_log2ratio (median-normalized log2ratio) + median = log2ratio. log2ratio is the value before log2 transform.
cl <- makeCluster(detectCores())
registerDoParallel(cl)
start_time_post <- Sys.time()
set.seed(20260121) # for reproducibility
result_missForest_post <- missForest(merged_normalized_log2ratio_6_mat_t, parallelize = "variables", verbose = TRUE)
end_time_post <- Sys.time()
stopCluster(cl)

spendintg_time_post <- end_time_post - start_time_post
spendintg_time_post # Time difference of 1.129523 days

gc(T,T,T)
gc(T,T,T)

# code 22
# Impute missing values using missForest for dataset AFTER normalization
# The matrix "merged_normalized_log2ratio_6_mat_t" is from the matrix "merged_normalized_log2ratio_6_mat".
# The matrix "merged_normalized_log2ratio_6_mat" is normalized_log2ratio (median-normalized log2ratio) + median = log2ratio. log2ratio is the value before log2 transform.
cl <- makeCluster(detectCores())
registerDoParallel(cl)
start_time_post <- Sys.time()
set.seed(20260121) # for reproducibility
result_missForest_post <- missForest(merged_normalized_log2ratio_6_mat_t, parallelize = "variables", verbose = TRUE)
end_time_post <- Sys.time()
stopCluster(cl)

spendintg_time_post <- end_time_post - start_time_post
spendintg_time_post # Time difference of 1.129523 days

gc(T,T,T)
gc(T,T,T)

# code 23
# save.image("/mnt/team4tb/Dropbox/Blog/2026 01 03 blog CPTAC/2026 01 10 CPTAC expression matrix_2.RData")

# code 23
# save.image("/mnt/team4tb/Dropbox/Blog/2026 01 03 blog CPTAC/2026 01 10 CPTAC expression matrix_2.RData")

# code 24

# check histogram

hist(merged_analysis_before_normalization_ratio_mat_t, breaks = 1000) # , ylim = c(0, 14000)
hist(result_missForest_pre_argu$ximp[1:100, 1:100], breaks = 1000)

hist(merged_analysis_before_normalization_ratio_mat_t[1:100, 1:100], breaks = 1000, xlim = c(0, 4))
hist(result_missForest_pre_argu$ximp[1:100, 1:100], breaks = 1000, xlim = c(0, 4))

hist(merged_log2ratio_6_mat_t[1:100, 1:100], breaks = 1000)
hist(result_missForest_pre$ximp[1:100, 1:100], breaks = 1000)

hist(merged_normalized_log2ratio_6_mat_t[1:100, 1:100], breaks = 1000)
hist(result_missForest_post$ximp[1:100, 1:100], breaks = 1000)

# code 24

# check histogram

hist(merged_analysis_before_normalization_ratio_mat_t, breaks = 1000) # , ylim = c(0, 14000)
hist(result_missForest_pre_argu$ximp[1:100, 1:100], breaks = 1000)

hist(merged_analysis_before_normalization_ratio_mat_t[1:100, 1:100], breaks = 1000, xlim = c(0, 4))
hist(result_missForest_pre_argu$ximp[1:100, 1:100], breaks = 1000, xlim = c(0, 4))

hist(merged_log2ratio_6_mat_t[1:100, 1:100], breaks = 1000)
hist(result_missForest_pre$ximp[1:100, 1:100], breaks = 1000)

hist(merged_normalized_log2ratio_6_mat_t[1:100, 1:100], breaks = 1000)
hist(result_missForest_post$ximp[1:100, 1:100], breaks = 1000)

# code 25

# merged_analysis_before_normalization_ratio_df # data frame, for "result_missForest_pre_argu" object
# merged_log2ratio_6_mat # matrix, for "result_missForest_pre" object
# merged_normalized_log2ratio_6 # data frame, for "result_missForest_post" object

## transform matrix "merged_log2ratio_6_mat" to data frame
merged_log2ratio_6_df <- data.frame(merged_log2ratio_6_mat)

## Check the dimension of each data frame. They had same dimension.
dim(merged_log2ratio_6_df) # 5208 5202
dim(merged_normalized_log2ratio_6) # 5208 5202
dim(merged_analysis_before_normalization_ratio_df) # 5208 5202

## Retrieve location of NA value in the dataframe "merged_normalized_log2ratio_6" (which is normalized intensity at CPTAC or submitter) and "merged_analysis_before_normalization_ratio_df"(which is the intensity before normalization). 

NA_location_pre_argu <- which(is.na(merged_analysis_before_normalization_ratio_df), arr.ind = TRUE)  %>% data.frame()
NA_location_pre <- which(is.na(merged_log2ratio_6_df), arr.ind = TRUE)  %>% data.frame()
NA_location_post <- which(is.na(merged_normalized_log2ratio_6), arr.ind = TRUE) %>% data.frame()


# check if the location of NA of dataframe "merged_normalized_log2ratio_6" and the dataframe "merged_analysis_before_normalization_ratio_df" is same or not. The results showed all location for both dataframe was same.

#
table(NA_location_pre$col %in% NA_location_pre_argu$col) # TRUE 748682, All true
table(NA_location_pre$row %in% NA_location_pre_argu$row) # TRUE 748682, All true
#
table(NA_location_pre_argu$col %in% NA_location_post$col) # TRUE 748682, All true
table(NA_location_pre_argu$row %in% NA_location_post$row) # TRUE 748682, All true
#
table(NA_location_pre$col %in% NA_location_post$col) # TRUE 748682, All true
table(NA_location_pre$row %in% NA_location_post$row) # TRUE 748682, All true

# code 25

# merged_analysis_before_normalization_ratio_df # data frame, for "result_missForest_pre_argu" object
# merged_log2ratio_6_mat # matrix, for "result_missForest_pre" object
# merged_normalized_log2ratio_6 # data frame, for "result_missForest_post" object

## transform matrix "merged_log2ratio_6_mat" to data frame
merged_log2ratio_6_df <- data.frame(merged_log2ratio_6_mat)

## Check the dimension of each data frame. They had same dimension.
dim(merged_log2ratio_6_df) # 5208 5202
dim(merged_normalized_log2ratio_6) # 5208 5202
dim(merged_analysis_before_normalization_ratio_df) # 5208 5202

## Retrieve location of NA value in the dataframe "merged_normalized_log2ratio_6" (which is normalized intensity at CPTAC or submitter) and "merged_analysis_before_normalization_ratio_df"(which is the intensity before normalization). 

NA_location_pre_argu <- which(is.na(merged_analysis_before_normalization_ratio_df), arr.ind = TRUE)  %>% data.frame()
NA_location_pre <- which(is.na(merged_log2ratio_6_df), arr.ind = TRUE)  %>% data.frame()
NA_location_post <- which(is.na(merged_normalized_log2ratio_6), arr.ind = TRUE) %>% data.frame()


# check if the location of NA of dataframe "merged_normalized_log2ratio_6" and the dataframe "merged_analysis_before_normalization_ratio_df" is same or not. The results showed all location for both dataframe was same.

#
table(NA_location_pre$col %in% NA_location_pre_argu$col) # TRUE 748682, All true
table(NA_location_pre$row %in% NA_location_pre_argu$row) # TRUE 748682, All true
#
table(NA_location_pre_argu$col %in% NA_location_post$col) # TRUE 748682, All true
table(NA_location_pre_argu$row %in% NA_location_post$row) # TRUE 748682, All true
#
table(NA_location_pre$col %in% NA_location_post$col) # TRUE 748682, All true
table(NA_location_pre$row %in% NA_location_post$row) # TRUE 748682, All true

# code 26

# Check what the data looks like.

view(merged_analysis_before_normalization_ratio_mat_t[1:118, 1:118]) # I do not know why I showed 1:118, but the area contained NA values.
view(result_missForest_pre_argu$ximp[1:100, 1:100])# missForest results of dataset before normalization and back to arguments

view(merged_log2ratio_6_mat_t[1:100, 1:100])
view(result_missForest_pre$ximp[1:100, 1:100])  # missForest results of dataset before normalizatiom.

view(merged_normalized_log2ratio_6_mat_t[1:100, 1:100])
view(result_missForest_post$ximp[1:100, 1:100])  # missForest results of normalized dataset 

# code 26

# Check what the data looks like.

view(merged_analysis_before_normalization_ratio_mat_t[1:118, 1:118]) # I do not know why I showed 1:118, but the area contained NA values.
view(result_missForest_pre_argu$ximp[1:100, 1:100])# missForest results of dataset before normalization and back to arguments

view(merged_log2ratio_6_mat_t[1:100, 1:100])
view(result_missForest_pre$ximp[1:100, 1:100])  # missForest results of dataset before normalizatiom.

view(merged_normalized_log2ratio_6_mat_t[1:100, 1:100])
view(result_missForest_post$ximp[1:100, 1:100])  # missForest results of normalized dataset 

# code 27

# Results of missForest were matrix. Transform the matrix into dataframe. Following is for the intensity before normalization ("result_missForest_pre"). The column and the row should respectively be cases and peptids and therefore it transposed by t().

# # cases_pre_argu_normalization
# # proteins_pre_argu_normalization
# # 
# # cases_pre_normalization_log
# # proteins_pre_normalization_log
# # 
# # cases_post_normalization
# # proteins_post_normalization

result_missForest_pre_argu_df <- t(result_missForest_pre_argu$ximp) %>% data.frame()
colnames(result_missForest_pre_argu_df) <- cases_pre_argu_normalization

# code 27

# Results of missForest were matrix. Transform the matrix into dataframe. Following is for the intensity before normalization ("result_missForest_pre"). The column and the row should respectively be cases and peptids and therefore it transposed by t().

# # cases_pre_argu_normalization
# # proteins_pre_argu_normalization
# # 
# # cases_pre_normalization_log
# # proteins_pre_normalization_log
# # 
# # cases_post_normalization
# # proteins_post_normalization

result_missForest_pre_argu_df <- t(result_missForest_pre_argu$ximp) %>% data.frame()
colnames(result_missForest_pre_argu_df) <- cases_pre_argu_normalization

# code 28

# Check if the median value of "result_missForest_pre_df" is closer to the median value that is calculated by CPTAC or submitters. 

## At first, check data.

view(result_missForest_pre_argu_df[1:200, 1:200]) # check data. This is data before normalization that was calculated with median. 
view(merged_mean_median_sd_2[1:3, 1:200]) # check data. This is mean, median and sd value that was used for median normalization.
view(merged_normalized_log2ratio_6_mat[1:200, 1:200]) # check data. This is normalized intensity by CPTAC or submitters.

## Match the cases in data frame "merged_mean_median_sd_2" with the cases in data frame "result_missForest_pre_df".

merged_mean_median_sd_3 <- merged_mean_median_sd_2 %>% dplyr::select(all_of(colnames(result_missForest_pre_argu_df)))
merged_mean_median_sd_3_mat <- as.matrix(merged_mean_median_sd_3)

view(merged_mean_median_sd_3_mat[1:3, 1:100]) # check data
view(t(result_missForest_pre_argu$ximp)[1:100, 1:100]) # check data

## log2 transform of missForest result before normalization.
result_missForest_pre_argu_log2 <- log2(t(result_missForest_pre_argu$ximp))
view(result_missForest_pre_argu_log2)[1:10, 1:10]# check data

## Normalization again. Later, check if the estimated NA value is close to median value by CPTAC or submitter. 
## Remind; log2(merged_normalized_log2ratio_6_mat) - median = normalized_log2ratio

# result_missForest_pre_argu_df_median_normalized <- sweep(x = result_missForest_pre_argu_log2,
#                                             MARGIN = 2,
#                                             STATS = merged_mean_median_sd_3_mat[2,],
#                                             FUN = "-")

# view(result_missForest_pre_argu_df_median_normalized[1:100, 1:100]) # check data

## Make sure the value of data frame "result_missForest_pre_df_median_normalization" is numeric. 
# result_missForest_pre_argu_df_median_normalized <- apply(result_missForest_pre_argu_df_median_normalized, 2, as.numeric)
result_missForest_pre_argu_log2 <- apply(result_missForest_pre_argu_log2, 2, as.numeric)

## Then the matrix "result_missForest_pre_df_median_normalization" transformed to data frame "result_missForest_pre_df_median_normalization_df"
# result_missForest_pre_argu_df_median_normalized_df <- data.frame(result_missForest_pre_argu_df_median_normalized)
result_missForest_pre_argu_log2_df <- data.frame(result_missForest_pre_argu_log2)

## Retrieve column name (case) and row name(protein).
colnames(result_missForest_pre_argu_log2_df) <- cases_pre_argu_normalization
rownames(result_missForest_pre_argu_log2_df) <- proteins_pre_argu_normalization

# # cases_pre_argu_normalization
# # proteins_pre_argu_normalization
# # 
# # cases_pre_normalization_log
# # proteins_pre_normalization_log
# # 
# # cases_post_normalization
# # proteins_post_normalization

table(cases_pre_argu_normalization %in% cases_pre_normalization_log) # TRUE 5202
table(proteins_pre_argu_normalization %in% proteins_pre_normalization_log) # TRUE 5208

table(cases_pre_argu_normalization %in% cases_post_normalization) # TRUE 5202
table(proteins_pre_argu_normalization %in% proteins_post_normalization) # TRUE 5208

table(cases_pre_normalization_log %in% cases_post_normalization) # TRUE 5202
table(proteins_pre_normalization_log %in% proteins_post_normalization) # TRUE 5208

# code 28

# Check if the median value of "result_missForest_pre_df" is closer to the median value that is calculated by CPTAC or submitters. 

## At first, check data.

view(result_missForest_pre_argu_df[1:200, 1:200]) # check data. This is data before normalization that was calculated with median. 
view(merged_mean_median_sd_2[1:3, 1:200]) # check data. This is mean, median and sd value that was used for median normalization.
view(merged_normalized_log2ratio_6_mat[1:200, 1:200]) # check data. This is normalized intensity by CPTAC or submitters.

## Match the cases in data frame "merged_mean_median_sd_2" with the cases in data frame "result_missForest_pre_df".

merged_mean_median_sd_3 <- merged_mean_median_sd_2 %>% dplyr::select(all_of(colnames(result_missForest_pre_argu_df)))
merged_mean_median_sd_3_mat <- as.matrix(merged_mean_median_sd_3)

view(merged_mean_median_sd_3_mat[1:3, 1:100]) # check data
view(t(result_missForest_pre_argu$ximp)[1:100, 1:100]) # check data

## log2 transform of missForest result before normalization.
result_missForest_pre_argu_log2 <- log2(t(result_missForest_pre_argu$ximp))
view(result_missForest_pre_argu_log2)[1:10, 1:10]# check data

## Normalization again. Later, check if the estimated NA value is close to median value by CPTAC or submitter. 
## Remind; log2(merged_normalized_log2ratio_6_mat) - median = normalized_log2ratio

# result_missForest_pre_argu_df_median_normalized <- sweep(x = result_missForest_pre_argu_log2,
#                                             MARGIN = 2,
#                                             STATS = merged_mean_median_sd_3_mat[2,],
#                                             FUN = "-")

# view(result_missForest_pre_argu_df_median_normalized[1:100, 1:100]) # check data

## Make sure the value of data frame "result_missForest_pre_df_median_normalization" is numeric. 
# result_missForest_pre_argu_df_median_normalized <- apply(result_missForest_pre_argu_df_median_normalized, 2, as.numeric)
result_missForest_pre_argu_log2 <- apply(result_missForest_pre_argu_log2, 2, as.numeric)

## Then the matrix "result_missForest_pre_df_median_normalization" transformed to data frame "result_missForest_pre_df_median_normalization_df"
# result_missForest_pre_argu_df_median_normalized_df <- data.frame(result_missForest_pre_argu_df_median_normalized)
result_missForest_pre_argu_log2_df <- data.frame(result_missForest_pre_argu_log2)

## Retrieve column name (case) and row name(protein).
colnames(result_missForest_pre_argu_log2_df) <- cases_pre_argu_normalization
rownames(result_missForest_pre_argu_log2_df) <- proteins_pre_argu_normalization

# # cases_pre_argu_normalization
# # proteins_pre_argu_normalization
# # 
# # cases_pre_normalization_log
# # proteins_pre_normalization_log
# # 
# # cases_post_normalization
# # proteins_post_normalization

table(cases_pre_argu_normalization %in% cases_pre_normalization_log) # TRUE 5202
table(proteins_pre_argu_normalization %in% proteins_pre_normalization_log) # TRUE 5208

table(cases_pre_argu_normalization %in% cases_post_normalization) # TRUE 5202
table(proteins_pre_argu_normalization %in% proteins_post_normalization) # TRUE 5208

table(cases_pre_normalization_log %in% cases_post_normalization) # TRUE 5202
table(proteins_pre_normalization_log %in% proteins_post_normalization) # TRUE 5208

# code 29

# merged_log2ratio_6_mat_t 

result_missForest_pre_df <- t(result_missForest_pre$ximp) %>% data.frame()
colnames(result_missForest_pre_df) <- cases_pre_normalization_log

# merged_mean_median_sd_4 <- merged_mean_median_sd_2 %>% dplyr::select(all_of(colnames(result_missForest_pre_df)))
# merged_mean_median_sd_4_mat <- as.matrix(merged_mean_median_sd_4)

result_missForest_pre_t <- t(result_missForest_pre$ximp)
view(result_missForest_pre_t[1:10, 1:10])# check data

## Normalization again. Later, check if the estimated NA value is close to median value by CPTAC or submitter. 
## Remind; log2(merged_normalized_log2ratio_6_mat) - median = normalized_log2ratio

# result_missForest_pre_df_median_normalized <- sweep(x = result_missForest_pre_t,
#                                             MARGIN = 2,
#                                             STATS = merged_mean_median_sd_4_mat[2,],
#                                             FUN = "-")

result_missForest_pre_t_df <- data.frame(result_missForest_pre_t)
colnames(result_missForest_pre_t_df) <- cases_pre_normalization_log

view(result_missForest_pre_t_df[1:100, 1:100])
view(cases_pre_normalization_log[1:100])
table(colnames(result_missForest_pre_t_df) %in% cases_pre_normalization_log) # TRUE 5144
table(colnames(result_missForest_pre_t_df) %in% cases_pre_normalization_log) # TRUE 5144
table(colnames(result_missForest_pre_df) %in% cases_pre_normalization_log) # TRUE 5144 

# colnames(result_missForest_pre_df_median_normalized_df) <- cases_pre_normalization_log
# rownames(result_missForest_pre_df_median_normalized_df) <- proteins_pre_normalization_log

# # cases_pre_argu_normalization
# # proteins_pre_argu_normalization
# # 
# # cases_pre_normalization_log
# # proteins_pre_normalization_log
# # 
# # cases_post_normalization
# # proteins_post_normalization

# code 29

# merged_log2ratio_6_mat_t 

result_missForest_pre_df <- t(result_missForest_pre$ximp) %>% data.frame()
colnames(result_missForest_pre_df) <- cases_pre_normalization_log

# merged_mean_median_sd_4 <- merged_mean_median_sd_2 %>% dplyr::select(all_of(colnames(result_missForest_pre_df)))
# merged_mean_median_sd_4_mat <- as.matrix(merged_mean_median_sd_4)

result_missForest_pre_t <- t(result_missForest_pre$ximp)
view(result_missForest_pre_t[1:10, 1:10])# check data

## Normalization again. Later, check if the estimated NA value is close to median value by CPTAC or submitter. 
## Remind; log2(merged_normalized_log2ratio_6_mat) - median = normalized_log2ratio

# result_missForest_pre_df_median_normalized <- sweep(x = result_missForest_pre_t,
#                                             MARGIN = 2,
#                                             STATS = merged_mean_median_sd_4_mat[2,],
#                                             FUN = "-")

result_missForest_pre_t_df <- data.frame(result_missForest_pre_t)
colnames(result_missForest_pre_t_df) <- cases_pre_normalization_log

view(result_missForest_pre_t_df[1:100, 1:100])
view(cases_pre_normalization_log[1:100])
table(colnames(result_missForest_pre_t_df) %in% cases_pre_normalization_log) # TRUE 5144
table(colnames(result_missForest_pre_t_df) %in% cases_pre_normalization_log) # TRUE 5144
table(colnames(result_missForest_pre_df) %in% cases_pre_normalization_log) # TRUE 5144 

# colnames(result_missForest_pre_df_median_normalized_df) <- cases_pre_normalization_log
# rownames(result_missForest_pre_df_median_normalized_df) <- proteins_pre_normalization_log

# # cases_pre_argu_normalization
# # proteins_pre_argu_normalization
# # 
# # cases_pre_normalization_log
# # proteins_pre_normalization_log
# # 
# # cases_post_normalization
# # proteins_post_normalization

# code 30

## matched the dimension (I mean, they are columns and rows) of "result_missForest_post_df" with data frame "result_missForest_pre_df_median_normalization", and then transformed it to data frame "result_missForest_post_df"

view(t(result_missForest_post$ximp)[1:100, 1:100]) # check data

merged_mean_median_sd_4 <- merged_mean_median_sd_2 %>% dplyr::select(all_of(colnames(t(result_missForest_post$ximp))))
merged_mean_median_sd_4_mat <- as.matrix(merged_mean_median_sd_4)

table(colnames(t(result_missForest_post$ximp)) %in% colnames(merged_mean_median_sd_4)) # TRUE 5144 

result_missForest_post_df <- sweep(x = t(result_missForest_post$ximp),
                                            MARGIN = 2,
                                            STATS = merged_mean_median_sd_4_mat[2,],
                                            FUN = "+")
## Retrieve column name for futher analysis.
colnames(result_missForest_post_df) <- cases_post_normalization

# code 30

## matched the dimension (I mean, they are columns and rows) of "result_missForest_post_df" with data frame "result_missForest_pre_df_median_normalization", and then transformed it to data frame "result_missForest_post_df"

view(t(result_missForest_post$ximp)[1:100, 1:100]) # check data

merged_mean_median_sd_4 <- merged_mean_median_sd_2 %>% dplyr::select(all_of(colnames(t(result_missForest_post$ximp))))
merged_mean_median_sd_4_mat <- as.matrix(merged_mean_median_sd_4)

table(colnames(t(result_missForest_post$ximp)) %in% colnames(merged_mean_median_sd_4)) # TRUE 5144 

result_missForest_post_df <- sweep(x = t(result_missForest_post$ximp),
                                            MARGIN = 2,
                                            STATS = merged_mean_median_sd_4_mat[2,],
                                            FUN = "+")
## Retrieve column name for futher analysis.
colnames(result_missForest_post_df) <- cases_post_normalization

# code 31

## check dimension of data frame.

dim(result_missForest_pre_argu_log2_df) # [1] 5208 5144 # This is centered with median value.
dim(result_missForest_pre_t_df) # [1] 5208 5144 # This is centered with median value.
dim(result_missForest_post_df) # [1] 5208 5144 # This is centered with median value at CPTAC or submitters.

## check if columns (cases) and rows (protein) between data frame "result_missForest_pre_argu_df_median_normalized_df" (intensity normalized by median value after imputation with missForest) and "result_missForest_post_df" (intensity of normalized by CPTAC or submitters and then missForest imputation)
table(colnames(result_missForest_pre_argu_log2_df) %in% colnames(result_missForest_pre_t_df)) # TRUE 5144 
table(colnames(result_missForest_pre_argu_log2_df) %in% colnames(result_missForest_post_df)) # TRUE 5144 
table(colnames(result_missForest_pre_t_df) %in% colnames(result_missForest_post_df)) # TRUE 5144 

table(rownames(result_missForest_pre_argu_log2_df) %in% rownames(result_missForest_pre_t_df)) # TRUE 5208 
table(rownames(result_missForest_pre_argu_log2_df) %in% rownames(result_missForest_post_df)) # TRUE 5208 
table(rownames(result_missForest_pre_t_df) %in% rownames(result_missForest_post_df)) # TRUE 5208 
## Prepare empty data frame for comparison of missForest results.
### column "pre" is for median value of data frame "result_missForest_pre_argu_df_median_normalized_df".
### column "post" is for median value of data frame "result_missForest_post_df".

check_after_missForest <- NA_location_post %>% mutate(
  pre_argu = NA,
  pre = NA,
  post = NA
)

# Put median value for each NA location in data frame "result_missForest_pre_argu_df_median_normalized_df" and "result_missForest_post_df" into column "pre" and "post" in data frame "check_after_missForest", respectively.

start_check_after_missForest <- Sys.time()
for(i in 1:nrow(NA_location_post)){
## for(i in 1:1000){ ## this line is just a test with small number of rows.
  print(paste0("Processing row ", i, " of ", nrow(NA_location_post)))
  check_after_missForest$pre_argu[i] <- result_missForest_pre_argu_log2_df[NA_location_post$row[i], NA_location_post$col[i]]
  check_after_missForest$pre[i] <- result_missForest_pre_t_df[NA_location_post$row[i], NA_location_post$col[i]]
  check_after_missForest$post[i] <- result_missForest_post_df[NA_location_post$row[i], NA_location_post$col[i]]
}
end_check_after_missForest <- Sys.time()

time_check_after_missForest <- end_check_after_missForest - start_check_after_missForest
time_check_after_missForest # Time difference of 18.79237 mins

# code 31

## check dimension of data frame.

dim(result_missForest_pre_argu_log2_df) # [1] 5208 5144 # This is centered with median value.
dim(result_missForest_pre_t_df) # [1] 5208 5144 # This is centered with median value.
dim(result_missForest_post_df) # [1] 5208 5144 # This is centered with median value at CPTAC or submitters.

## check if columns (cases) and rows (protein) between data frame "result_missForest_pre_argu_df_median_normalized_df" (intensity normalized by median value after imputation with missForest) and "result_missForest_post_df" (intensity of normalized by CPTAC or submitters and then missForest imputation)
table(colnames(result_missForest_pre_argu_log2_df) %in% colnames(result_missForest_pre_t_df)) # TRUE 5144 
table(colnames(result_missForest_pre_argu_log2_df) %in% colnames(result_missForest_post_df)) # TRUE 5144 
table(colnames(result_missForest_pre_t_df) %in% colnames(result_missForest_post_df)) # TRUE 5144 

table(rownames(result_missForest_pre_argu_log2_df) %in% rownames(result_missForest_pre_t_df)) # TRUE 5208 
table(rownames(result_missForest_pre_argu_log2_df) %in% rownames(result_missForest_post_df)) # TRUE 5208 
table(rownames(result_missForest_pre_t_df) %in% rownames(result_missForest_post_df)) # TRUE 5208 
## Prepare empty data frame for comparison of missForest results.
### column "pre" is for median value of data frame "result_missForest_pre_argu_df_median_normalized_df".
### column "post" is for median value of data frame "result_missForest_post_df".

check_after_missForest <- NA_location_post %>% mutate(
  pre_argu = NA,
  pre = NA,
  post = NA
)

# Put median value for each NA location in data frame "result_missForest_pre_argu_df_median_normalized_df" and "result_missForest_post_df" into column "pre" and "post" in data frame "check_after_missForest", respectively.

start_check_after_missForest <- Sys.time()
for(i in 1:nrow(NA_location_post)){
## for(i in 1:1000){ ## this line is just a test with small number of rows.
  print(paste0("Processing row ", i, " of ", nrow(NA_location_post)))
  check_after_missForest$pre_argu[i] <- result_missForest_pre_argu_log2_df[NA_location_post$row[i], NA_location_post$col[i]]
  check_after_missForest$pre[i] <- result_missForest_pre_t_df[NA_location_post$row[i], NA_location_post$col[i]]
  check_after_missForest$post[i] <- result_missForest_post_df[NA_location_post$row[i], NA_location_post$col[i]]
}
end_check_after_missForest <- Sys.time()

time_check_after_missForest <- end_check_after_missForest - start_check_after_missForest
time_check_after_missForest # Time difference of 18.79237 mins

# code 32

## median and modal value of imputed NA value. modal value was calculated by modal() in raster package.
### median value of imputed NA value
median(check_after_missForest$pre_argu) # -0.124715
median(check_after_missForest$pre) # -0.2962415
median(check_after_missForest$post) # -0.2952342

### mean value of imputed NA value
mean(check_after_missForest$pre_argu) # -0.1505831
mean(check_after_missForest$pre) # -0.3874793
mean(check_after_missForest$post) # -0.3849859


## range of the imputed NA value
quantile(check_after_missForest$pre_argu, probs = c(0.25, 0.5, 0.75)) # -0.5652116 -0.1247150  0.2215660
quantile(check_after_missForest$pre, probs = c(0.25, 0.5, 0.75)) # -0.76488152 -0.29624152  0.05727549
quantile(check_after_missForest$post, probs = c(0.25, 0.5, 0.75)) # -0.76853020 -0.29523424  0.05929829 

# IQR
IQR(check_after_missForest$pre_argu) # 0.7867776
IQR(check_after_missForest$pre) # 0.822157
IQR(check_after_missForest$post) # 0.8278285

## SD
sd(check_after_missForest$pre_argu) # 0.9131712
sd(check_after_missForest$pre) # 0.7586585
sd(check_after_missForest$post) # 0.8080771

## %CV
sd(check_after_missForest$pre_argu)/mean(check_after_missForest$pre_argu) # -6.064233 # mean value is below 0.
sd(check_after_missForest$pre)/mean(check_after_missForest$pre) # -1.957933 # mean value is below 0.
sd(check_after_missForest$post)/mean(check_after_missForest$post) # -2.098979 # mean value is below 0.

## range
range(check_after_missForest$pre_argu) # -6.474371 25.291331
max(check_after_missForest$pre_argu) - min(check_after_missForest$pre_argu) # 31.7657

range(check_after_missForest$pre) # -11.205002   8.976766
max(check_after_missForest$pre) - min(check_after_missForest$pre) # 20.18177

range(check_after_missForest$post) # -8.107361  9.021220
max(check_after_missForest$post) - min(check_after_missForest$post) # 17.12858

## histogram
hist_check_after_missForest_pre_argu <- hist(check_after_missForest$pre_argu, breaks = 10000)
hist_check_after_missForest_pre <- hist(check_after_missForest$pre, breaks = 10000)
hist_check_after_missForest_post <- hist(check_after_missForest$post, breaks = 10000)

hist_check_after_missForest_pre_argu_df <- data.frame(mids = hist_check_after_missForest_pre_argu$mids, counts = hist_check_after_missForest_pre_argu$counts)
hist_check_after_missForest_pre_df <- data.frame(mids = hist_check_after_missForest_pre$mids, counts = hist_check_after_missForest_pre$counts)
hist_check_after_missForest_post_df <- data.frame(mids = hist_check_after_missForest_post$mids, counts = hist_check_after_missForest_post$counts)

hist_check_after_missForest_pre_argu_df[hist_check_after_missForest_pre_argu_df$counts == max(hist_check_after_missForest_pre_argu_df$counts), ]$counts # 3551
hist_check_after_missForest_pre_argu_df[hist_check_after_missForest_pre_argu_df$counts == max(hist_check_after_missForest_pre_argu_df$counts), ]$mids # -0.0075

hist_check_after_missForest_pre_df[hist_check_after_missForest_pre_df$counts == max(hist_check_after_missForest_pre_df$counts), ]$counts # 1425
hist_check_after_missForest_pre_df[hist_check_after_missForest_pre_df$counts == max(hist_check_after_missForest_pre_df$counts), ]$mids # -0.005

hist_check_after_missForest_post_df[hist_check_after_missForest_post_df$counts == max(hist_check_after_missForest_post_df$counts), ]$counts # 1463
hist_check_after_missForest_post_df[hist_check_after_missForest_post_df$counts == max(hist_check_after_missForest_post_df$counts), ]$mids # 0.009

# code 32

## median and modal value of imputed NA value. modal value was calculated by modal() in raster package.
### median value of imputed NA value
median(check_after_missForest$pre_argu) # -0.124715
median(check_after_missForest$pre) # -0.2962415
median(check_after_missForest$post) # -0.2952342

### mean value of imputed NA value
mean(check_after_missForest$pre_argu) # -0.1505831
mean(check_after_missForest$pre) # -0.3874793
mean(check_after_missForest$post) # -0.3849859


## range of the imputed NA value
quantile(check_after_missForest$pre_argu, probs = c(0.25, 0.5, 0.75)) # -0.5652116 -0.1247150  0.2215660
quantile(check_after_missForest$pre, probs = c(0.25, 0.5, 0.75)) # -0.76488152 -0.29624152  0.05727549
quantile(check_after_missForest$post, probs = c(0.25, 0.5, 0.75)) # -0.76853020 -0.29523424  0.05929829 

# IQR
IQR(check_after_missForest$pre_argu) # 0.7867776
IQR(check_after_missForest$pre) # 0.822157
IQR(check_after_missForest$post) # 0.8278285

## SD
sd(check_after_missForest$pre_argu) # 0.9131712
sd(check_after_missForest$pre) # 0.7586585
sd(check_after_missForest$post) # 0.8080771

## %CV
sd(check_after_missForest$pre_argu)/mean(check_after_missForest$pre_argu) # -6.064233 # mean value is below 0.
sd(check_after_missForest$pre)/mean(check_after_missForest$pre) # -1.957933 # mean value is below 0.
sd(check_after_missForest$post)/mean(check_after_missForest$post) # -2.098979 # mean value is below 0.

## range
range(check_after_missForest$pre_argu) # -6.474371 25.291331
max(check_after_missForest$pre_argu) - min(check_after_missForest$pre_argu) # 31.7657

range(check_after_missForest$pre) # -11.205002   8.976766
max(check_after_missForest$pre) - min(check_after_missForest$pre) # 20.18177

range(check_after_missForest$post) # -8.107361  9.021220
max(check_after_missForest$post) - min(check_after_missForest$post) # 17.12858

## histogram
hist_check_after_missForest_pre_argu <- hist(check_after_missForest$pre_argu, breaks = 10000)
hist_check_after_missForest_pre <- hist(check_after_missForest$pre, breaks = 10000)
hist_check_after_missForest_post <- hist(check_after_missForest$post, breaks = 10000)

hist_check_after_missForest_pre_argu_df <- data.frame(mids = hist_check_after_missForest_pre_argu$mids, counts = hist_check_after_missForest_pre_argu$counts)
hist_check_after_missForest_pre_df <- data.frame(mids = hist_check_after_missForest_pre$mids, counts = hist_check_after_missForest_pre$counts)
hist_check_after_missForest_post_df <- data.frame(mids = hist_check_after_missForest_post$mids, counts = hist_check_after_missForest_post$counts)

hist_check_after_missForest_pre_argu_df[hist_check_after_missForest_pre_argu_df$counts == max(hist_check_after_missForest_pre_argu_df$counts), ]$counts # 3551
hist_check_after_missForest_pre_argu_df[hist_check_after_missForest_pre_argu_df$counts == max(hist_check_after_missForest_pre_argu_df$counts), ]$mids # -0.0075

hist_check_after_missForest_pre_df[hist_check_after_missForest_pre_df$counts == max(hist_check_after_missForest_pre_df$counts), ]$counts # 1425
hist_check_after_missForest_pre_df[hist_check_after_missForest_pre_df$counts == max(hist_check_after_missForest_pre_df$counts), ]$mids # -0.005

hist_check_after_missForest_post_df[hist_check_after_missForest_post_df$counts == max(hist_check_after_missForest_post_df$counts), ]$counts # 1463
hist_check_after_missForest_post_df[hist_check_after_missForest_post_df$counts == max(hist_check_after_missForest_post_df$counts), ]$mids # 0.009

# code 33

## Check distribution of imputed values.
### histogram of imputed value. the data was stored at list "hist_imputed_pre" and "hist_imputed_post". 
hist_imputed_pre_argu <- hist(check_after_missForest$pre_argu, 
     breaks = 10000,
     main = "Distribution of imputed values by missForest (pre normalization and back to arguments)",
     xlab = "Expression level")

hist_imputed_pre <- hist(check_after_missForest$pre, 
     breaks = 10000,
     main = "Distribution of imputed values by missForest (pre normalization)",
     xlab = "Expression level")

hist_imputed_post <- hist(check_after_missForest$post, 
     breaks = 10000,
     main = "Distribution of imputed values by missForest (post normalization)",
     xlab = "Expression level")

### Distribution of the imputed NA values were stored at data frame "hist_imputed_pre_df" and "hist_imputed_post_df".
hist_imputed_pre_argu_df <- data.frame(mids = hist_imputed_pre_argu$mids, counts = hist_imputed_pre_argu$counts)
hist_imputed_pre_df <- data.frame(mids = hist_imputed_pre$mids, counts = hist_imputed_pre$counts)
hist_imputed_post_df <- data.frame(mids = hist_imputed_post$mids, counts = hist_imputed_post$counts)

### Take modal value for missForest results using the intensity before normalization.
hist_imputed_pre_argu_df[hist_imputed_pre_argu_df$counts == max(hist_imputed_pre_argu_df$counts), ]$counts # 3551
hist_imputed_pre_argu_df[hist_imputed_pre_argu_df$counts == max(hist_imputed_pre_argu_df$counts), ]$mids # -0.0075 # modal values when the break point is 10000

### Take modal value for missForest results using the intensity before normalization.

hist_imputed_pre_df[hist_imputed_pre_df$counts == max(hist_imputed_pre_df$counts), ]$counts # 1425
hist_imputed_pre_df[hist_imputed_pre_df$counts == max(hist_imputed_pre_df$counts), ]$mids # -0.005 # modal values when the break point is 10000

plot(hist_imputed_pre) # NA values for 9014 cases were imputed around 0 (modal value was 0.0175)

### Take modal value for missForest results using the normalized intensity.

hist_imputed_post_df[hist_imputed_post_df$counts == max(hist_imputed_post_df$counts), ]$counts # 1463
hist_imputed_post_df[hist_imputed_post_df$counts == max(hist_imputed_post_df$counts), ]$mids # 0.009 # modal values when the break point is 10000

plot(hist_imputed_post_df) # NA values for 1427 cases were imputed around 0 (modal value was 0.009)

## After prepare the data frame, clear unused memory as possible as it can be.
gc(T,T,T)
gc(T,T,T)

## Proteins contains the parcentage of NA around 0.1 were selected. 
NA_more_than_0.1_gene <- merged_normalized_log2ratio_w_percent_na %>% filter(percent_na <= 0.1 & percent_na >= 0.099) %>% dplyr::select(percent_na)

## WDR59 was selected as a representative proteins for checking which dataset (missForest using before or after normalization) os better. 
# result_missForest_pre_df_median_normalization_df["WDR59", ] %>% view()

## Checked with histogram. However, it is difficult to decide which dataset is better for missForest immputation.

### result_missForest_pre_argu_df_median_normalized_df
### result_missForest_pre_df_median_normalized_df
### result_missForest_post_df

hist(t(result_missForest_pre_argu_log2_df["WDR59", ]), 
     breaks = 1000,
     main = "Distribution after missForest imputation then normalizaion again (pre normalization)",
     xlab = "Expression level")

hist(t(result_missForest_pre_t_df["WDR59", ]), 
     breaks = 1000,
     main = "Distribution after missForest imputation then normalizaion again (pre normalization)",
     xlab = "Expression level")

hist(t(result_missForest_post_df["WDR59", ]), 
     breaks = 1000,
     main = "Distribution after missForest imputation (post normalization)",
     xlab = "Expression level")



## Other way to decide which dataset is better for missFortest imputation is comparison of difference between median values by CPTAC, missForest result using the dataset before normalization (data frame "result_missForest_pre_df_median_normalization_df"; the dataset of which is normalized again with median value from missForest result) and after normalization (data frame "result_missForest_post_df"; the dataset is median normalization by CPTAC or submitters).

### Check difference in the median value between data frame "result_missForest_pre_df_median_normalization_df" and "merged_mean_median_sd_3" (median value by CPTAC or submitter). The first column was used as the representative. but It was also difficult to decide which is better.

median(result_missForest_pre_argu_log2_df[,1]) # -1.055684
median(result_missForest_pre_t_df[,1]) # -1.075251 # median value of median-normalized data should be close to approximately 0. It looks OK. 
median(result_missForest_post_df[,1]) # -1.073708 # I think it is weird. Why the value is close to median value by CPTAC or submitters? 
median(merged_mean_median_sd_3[,1]) # -1.024735 # the value is median value that was calculated by CPTAC or submitters.

## The median value of first column of data frame "result_missForest_post_df" (-1.075444) is close to the median value by CPTAC or submitters (-1.024735)... I think it is weird. Why median value of median normalized intensity was close to median by CPTAC or submitters. The data frame "result_missForest_post_df" is  MEDIAN NORMALIZED DATA and therefore the median value should be closed to 0 when the NA value means that the peptid was difficult to be detected and it resulted in the super low intensity. 

## median_after_missForest <- data.frame(matrix(nrow = 1, ncol = ncol(result_missForest_post_df)))

## Next, calculate median value over the cases in both data frame "result_missForest_pre_df_median_normalization_df" and "result_missForest_post_df".

### required data;
#### matrix "result_missForest_pre_log2"
#### data frame "result_missForest_post_df"; the data frame is already normalized. it should be gotten back to the value before normalization, and then calculate median value.
#### data frame "merged_mean_median_sd_3"; this is mean, median, and sd for case of interest calculated by CPTAC or submitters.

# code 33

## Check distribution of imputed values.
### histogram of imputed value. the data was stored at list "hist_imputed_pre" and "hist_imputed_post". 
hist_imputed_pre_argu <- hist(check_after_missForest$pre_argu, 
     breaks = 10000,
     main = "Distribution of imputed values by missForest (pre normalization and back to arguments)",
     xlab = "Expression level")

hist_imputed_pre <- hist(check_after_missForest$pre, 
     breaks = 10000,
     main = "Distribution of imputed values by missForest (pre normalization)",
     xlab = "Expression level")

hist_imputed_post <- hist(check_after_missForest$post, 
     breaks = 10000,
     main = "Distribution of imputed values by missForest (post normalization)",
     xlab = "Expression level")

### Distribution of the imputed NA values were stored at data frame "hist_imputed_pre_df" and "hist_imputed_post_df".
hist_imputed_pre_argu_df <- data.frame(mids = hist_imputed_pre_argu$mids, counts = hist_imputed_pre_argu$counts)
hist_imputed_pre_df <- data.frame(mids = hist_imputed_pre$mids, counts = hist_imputed_pre$counts)
hist_imputed_post_df <- data.frame(mids = hist_imputed_post$mids, counts = hist_imputed_post$counts)

### Take modal value for missForest results using the intensity before normalization.
hist_imputed_pre_argu_df[hist_imputed_pre_argu_df$counts == max(hist_imputed_pre_argu_df$counts), ]$counts # 3551
hist_imputed_pre_argu_df[hist_imputed_pre_argu_df$counts == max(hist_imputed_pre_argu_df$counts), ]$mids # -0.0075 # modal values when the break point is 10000

### Take modal value for missForest results using the intensity before normalization.

hist_imputed_pre_df[hist_imputed_pre_df$counts == max(hist_imputed_pre_df$counts), ]$counts # 1425
hist_imputed_pre_df[hist_imputed_pre_df$counts == max(hist_imputed_pre_df$counts), ]$mids # -0.005 # modal values when the break point is 10000

plot(hist_imputed_pre) # NA values for 9014 cases were imputed around 0 (modal value was 0.0175)

### Take modal value for missForest results using the normalized intensity.

hist_imputed_post_df[hist_imputed_post_df$counts == max(hist_imputed_post_df$counts), ]$counts # 1463
hist_imputed_post_df[hist_imputed_post_df$counts == max(hist_imputed_post_df$counts), ]$mids # 0.009 # modal values when the break point is 10000

plot(hist_imputed_post_df) # NA values for 1427 cases were imputed around 0 (modal value was 0.009)

## After prepare the data frame, clear unused memory as possible as it can be.
gc(T,T,T)
gc(T,T,T)

## Proteins contains the parcentage of NA around 0.1 were selected. 
NA_more_than_0.1_gene <- merged_normalized_log2ratio_w_percent_na %>% filter(percent_na <= 0.1 & percent_na >= 0.099) %>% dplyr::select(percent_na)

## WDR59 was selected as a representative proteins for checking which dataset (missForest using before or after normalization) os better. 
# result_missForest_pre_df_median_normalization_df["WDR59", ] %>% view()

## Checked with histogram. However, it is difficult to decide which dataset is better for missForest immputation.

### result_missForest_pre_argu_df_median_normalized_df
### result_missForest_pre_df_median_normalized_df
### result_missForest_post_df

hist(t(result_missForest_pre_argu_log2_df["WDR59", ]), 
     breaks = 1000,
     main = "Distribution after missForest imputation then normalizaion again (pre normalization)",
     xlab = "Expression level")

hist(t(result_missForest_pre_t_df["WDR59", ]), 
     breaks = 1000,
     main = "Distribution after missForest imputation then normalizaion again (pre normalization)",
     xlab = "Expression level")

hist(t(result_missForest_post_df["WDR59", ]), 
     breaks = 1000,
     main = "Distribution after missForest imputation (post normalization)",
     xlab = "Expression level")



## Other way to decide which dataset is better for missFortest imputation is comparison of difference between median values by CPTAC, missForest result using the dataset before normalization (data frame "result_missForest_pre_df_median_normalization_df"; the dataset of which is normalized again with median value from missForest result) and after normalization (data frame "result_missForest_post_df"; the dataset is median normalization by CPTAC or submitters).

### Check difference in the median value between data frame "result_missForest_pre_df_median_normalization_df" and "merged_mean_median_sd_3" (median value by CPTAC or submitter). The first column was used as the representative. but It was also difficult to decide which is better.

median(result_missForest_pre_argu_log2_df[,1]) # -1.055684
median(result_missForest_pre_t_df[,1]) # -1.075251 # median value of median-normalized data should be close to approximately 0. It looks OK. 
median(result_missForest_post_df[,1]) # -1.073708 # I think it is weird. Why the value is close to median value by CPTAC or submitters? 
median(merged_mean_median_sd_3[,1]) # -1.024735 # the value is median value that was calculated by CPTAC or submitters.

## The median value of first column of data frame "result_missForest_post_df" (-1.075444) is close to the median value by CPTAC or submitters (-1.024735)... I think it is weird. Why median value of median normalized intensity was close to median by CPTAC or submitters. The data frame "result_missForest_post_df" is  MEDIAN NORMALIZED DATA and therefore the median value should be closed to 0 when the NA value means that the peptid was difficult to be detected and it resulted in the super low intensity. 

## median_after_missForest <- data.frame(matrix(nrow = 1, ncol = ncol(result_missForest_post_df)))

## Next, calculate median value over the cases in both data frame "result_missForest_pre_df_median_normalization_df" and "result_missForest_post_df".

### required data;
#### matrix "result_missForest_pre_log2"
#### data frame "result_missForest_post_df"; the data frame is already normalized. it should be gotten back to the value before normalization, and then calculate median value.
#### data frame "merged_mean_median_sd_3"; this is mean, median, and sd for case of interest calculated by CPTAC or submitters.

# code 34

# temp <- list()
# temp[[1]] <- data.frame(t(apply(result_missForest_pre_df_median_normalization_df, 2, median))) # apply(MARGIN = 2) means apply the function (here, median) to each column. If MARGIN = 1, the function will be applied to each row.
# temp[[2]] <- data.frame(t(apply(result_missForest_post_df, 2, median))) # apply(MARGIN = 2) means apply the function (here, median) to each column. If MARGIN = 1, the function will be applied to each row.
# names(temp) <- c("median_pre", "median_post")

# result_missForest_post_mat <- as.matrix(result_missForest_post_df)
# view(result_missForest_post_mat[1:100, 1:100]) # Check data
# result_missForest_post_mat_2 <- sweep(x = result_missForest_post_mat, MARGIN = 2, STATS = merged_mean_median_sd_3_mat[2,], FUN = "+")
# view(result_missForest_post_mat_2[1:100, 1:100]) # Check data

temp <- list()
temp[[1]] <- data.frame(t(apply(result_missForest_pre_argu_log2_df, 2, median)))
temp[[2]] <- data.frame(t(apply(result_missForest_pre_t_df, 2, median)))
temp[[3]] <- data.frame(t(apply(result_missForest_post_df, 2, median)))
names(temp) <- c("median_pre_argu", "median_pre", "median_post")



## marge the columns. 
median_after_missForest <- bind_rows(temp)

## name the columns.
rownames(median_after_missForest) <- names(temp)
colnames(median_after_missForest) <- colnames(result_missForest_post_df)

## Add median value by CPTAC or submitters to data frame "median_after_missForest". 
median_after_missForest <- bind_rows(median_after_missForest, merged_mean_median_sd_3[2,]) %>% t() %>% data.frame()

## column "median_pre"; median value that is calculated with imputed NA value by missForest of intensity before normalization
## column "median_post"; median value that is calculated with imputed NA value by missForest of normalized intensity.
## column "Median"; median value that is calculated by CPTAC and submitters.

# 
median_after_missForest$difference_pre_argu <- abs(median_after_missForest$median_pre_argu - median_after_missForest$Median)
median_after_missForest$difference_pre <- abs(median_after_missForest$median_pre - median_after_missForest$Median)
median_after_missForest$difference_post <- abs(median_after_missForest$median_post - median_after_missForest$Median)

median_after_missForest_difference <- colnames(median_after_missForest)[5:7]
median_after_missForest$smaller_difference_is <- apply(median_after_missForest[,5:7], 1, function(x) median_after_missForest_difference[which.min(x)])

table(median_after_missForest$smaller_difference_is)

## Following table is the benchmark of missForest imputation. 
## missForest using arguments (inversed log2 value) before normalization got the median value that was more the closer to median value of CPTAC or submitters.
##     difference_post      difference_pre difference_pre_argu 
##               1186                1612                2404 
                                 
                                 
# median_after_missForest <- median_after_missForest %>% mutate(
#   smaller_differece_is = case_when(
#     difference_pre < difference_post ~ "pre",
#     difference_pre > difference_post ~ "post",
#     difference_pre == difference_post ~ "even",
#     TRUE ~ NA)
#   )

# table(median_after_missForest$smaller_differece_is)

# even post  pre 
#    2  218 4924 

# code 34

# temp <- list()
# temp[[1]] <- data.frame(t(apply(result_missForest_pre_df_median_normalization_df, 2, median))) # apply(MARGIN = 2) means apply the function (here, median) to each column. If MARGIN = 1, the function will be applied to each row.
# temp[[2]] <- data.frame(t(apply(result_missForest_post_df, 2, median))) # apply(MARGIN = 2) means apply the function (here, median) to each column. If MARGIN = 1, the function will be applied to each row.
# names(temp) <- c("median_pre", "median_post")

# result_missForest_post_mat <- as.matrix(result_missForest_post_df)
# view(result_missForest_post_mat[1:100, 1:100]) # Check data
# result_missForest_post_mat_2 <- sweep(x = result_missForest_post_mat, MARGIN = 2, STATS = merged_mean_median_sd_3_mat[2,], FUN = "+")
# view(result_missForest_post_mat_2[1:100, 1:100]) # Check data

temp <- list()
temp[[1]] <- data.frame(t(apply(result_missForest_pre_argu_log2_df, 2, median)))
temp[[2]] <- data.frame(t(apply(result_missForest_pre_t_df, 2, median)))
temp[[3]] <- data.frame(t(apply(result_missForest_post_df, 2, median)))
names(temp) <- c("median_pre_argu", "median_pre", "median_post")



## marge the columns. 
median_after_missForest <- bind_rows(temp)

## name the columns.
rownames(median_after_missForest) <- names(temp)
colnames(median_after_missForest) <- colnames(result_missForest_post_df)

## Add median value by CPTAC or submitters to data frame "median_after_missForest". 
median_after_missForest <- bind_rows(median_after_missForest, merged_mean_median_sd_3[2,]) %>% t() %>% data.frame()

## column "median_pre"; median value that is calculated with imputed NA value by missForest of intensity before normalization
## column "median_post"; median value that is calculated with imputed NA value by missForest of normalized intensity.
## column "Median"; median value that is calculated by CPTAC and submitters.

# 
median_after_missForest$difference_pre_argu <- abs(median_after_missForest$median_pre_argu - median_after_missForest$Median)
median_after_missForest$difference_pre <- abs(median_after_missForest$median_pre - median_after_missForest$Median)
median_after_missForest$difference_post <- abs(median_after_missForest$median_post - median_after_missForest$Median)

median_after_missForest_difference <- colnames(median_after_missForest)[5:7]
median_after_missForest$smaller_difference_is <- apply(median_after_missForest[,5:7], 1, function(x) median_after_missForest_difference[which.min(x)])

table(median_after_missForest$smaller_difference_is)

## Following table is the benchmark of missForest imputation. 
## missForest using arguments (inversed log2 value) before normalization got the median value that was more the closer to median value of CPTAC or submitters.
##     difference_post      difference_pre difference_pre_argu 
##               1186                1612                2404 
                                 
                                 
# median_after_missForest <- median_after_missForest %>% mutate(
#   smaller_differece_is = case_when(
#     difference_pre < difference_post ~ "pre",
#     difference_pre > difference_post ~ "post",
#     difference_pre == difference_post ~ "even",
#     TRUE ~ NA)
#   )

# table(median_after_missForest$smaller_differece_is)

# even post  pre 
#    2  218 4924 

# code 35

hist(as.matrix(result_missForest_pre_argu_log2_df), breaks = 10000)
hist(as.matrix(result_missForest_pre_t_df), breaks = 10000)
hist(as.matrix(result_missForest_post_df), breaks = 10000)

## range of the imputed NA value
quantile(as.matrix(result_missForest_pre_argu_log2_df), probs = c(0.25, 0.5, 0.75)) # -0.5219597 -0.0705235  0.3350008
quantile(as.matrix(result_missForest_pre_t_df),, probs = c(0.25, 0.5, 0.75)) # -0.52889415 -0.07524554  0.33078758
quantile(as.matrix(result_missForest_post_df),, probs = c(0.25, 0.5, 0.75)) # -0.52895082 -0.07527012  0.33085971

## need wilcox.test and ks.test for comparison between before and after imputation.
ks.test(merged_analysis_before_normalization_ratio_mat_t, result_missForest_pre_argu$ximp) # p-value < 2.2e-16
wilcox.test(merged_analysis_before_normalization_ratio_mat_t, result_missForest_pre_argu$ximp) # p-value < 2.2e-16

ks.test(merged_log2ratio_6_mat_t, result_missForest_pre$ximp) # p-value < 2.2e-16
wilcox.test(merged_log2ratio_6_mat_t, result_missForest_pre$ximp) # p-value < 2.2e-16

ks.test(merged_normalized_log2ratio_6_mat_t, result_missForest_post$ximp) # p-value < 2.2e-16
wilcox.test(merged_normalized_log2ratio_6_mat_t, result_missForest_post$ximp) # p-value = 1.273e-07

#

hist(merged_analysis_before_normalization_ratio_mat_t, breaks = 100)
hist(result_missForest_pre_argu$ximp, breaks = 100)

hist(merged_analysis_before_normalization_ratio_mat_t[1:100, 1:100], breaks = 10000)
hist(result_missForest_pre_argu$ximp[1:100, 1:100], breaks = 10000)

max(merged_analysis_before_normalization_ratio_mat_t, na.rm = TRUE) # 327053101
max(result_missForest_pre_argu$ximp, na.rm = TRUE) # 327053101

#

hist(merged_log2ratio_6_mat_t, breaks = 100)
hist(result_missForest_pre$ximp, breaks = 100)

max(merged_log2ratio_6_mat_t, na.rm = TRUE) # 28.28495
max(result_missForest_pre$ximp, na.rm = TRUE) # 28.28495

#

hist(merged_normalized_log2ratio_6_mat_t, breaks = 100)
hist(result_missForest_post$ximp, breaks = 100)

max(merged_normalized_log2ratio_6_mat_t, na.rm = TRUE) # 28.95969
max(result_missForest_post$ximp, na.rm = TRUE) # 28.95969

# check maximum value.
max(merged_analysis_before_normalization_ratio_mat_t, na.rm = TRUE) # 327053101
range(merged_analysis_before_normalization_ratio_mat_t, na.rm = TRUE) # 3.510291e-09 3.270531e+08

max(merged_log2ratio_6_mat_t, na.rm = TRUE) # 28.28495
2^max(merged_log2ratio_6_mat_t, na.rm = TRUE) # 327053101

# code 35

hist(as.matrix(result_missForest_pre_argu_log2_df), breaks = 10000)
hist(as.matrix(result_missForest_pre_t_df), breaks = 10000)
hist(as.matrix(result_missForest_post_df), breaks = 10000)

## range of the imputed NA value
quantile(as.matrix(result_missForest_pre_argu_log2_df), probs = c(0.25, 0.5, 0.75)) # -0.5219597 -0.0705235  0.3350008
quantile(as.matrix(result_missForest_pre_t_df),, probs = c(0.25, 0.5, 0.75)) # -0.52889415 -0.07524554  0.33078758
quantile(as.matrix(result_missForest_post_df),, probs = c(0.25, 0.5, 0.75)) # -0.52895082 -0.07527012  0.33085971

## need wilcox.test and ks.test for comparison between before and after imputation.
ks.test(merged_analysis_before_normalization_ratio_mat_t, result_missForest_pre_argu$ximp) # p-value < 2.2e-16
wilcox.test(merged_analysis_before_normalization_ratio_mat_t, result_missForest_pre_argu$ximp) # p-value < 2.2e-16

ks.test(merged_log2ratio_6_mat_t, result_missForest_pre$ximp) # p-value < 2.2e-16
wilcox.test(merged_log2ratio_6_mat_t, result_missForest_pre$ximp) # p-value < 2.2e-16

ks.test(merged_normalized_log2ratio_6_mat_t, result_missForest_post$ximp) # p-value < 2.2e-16
wilcox.test(merged_normalized_log2ratio_6_mat_t, result_missForest_post$ximp) # p-value = 1.273e-07

#

hist(merged_analysis_before_normalization_ratio_mat_t, breaks = 100)
hist(result_missForest_pre_argu$ximp, breaks = 100)

hist(merged_analysis_before_normalization_ratio_mat_t[1:100, 1:100], breaks = 10000)
hist(result_missForest_pre_argu$ximp[1:100, 1:100], breaks = 10000)

max(merged_analysis_before_normalization_ratio_mat_t, na.rm = TRUE) # 327053101
max(result_missForest_pre_argu$ximp, na.rm = TRUE) # 327053101

#

hist(merged_log2ratio_6_mat_t, breaks = 100)
hist(result_missForest_pre$ximp, breaks = 100)

max(merged_log2ratio_6_mat_t, na.rm = TRUE) # 28.28495
max(result_missForest_pre$ximp, na.rm = TRUE) # 28.28495

#

hist(merged_normalized_log2ratio_6_mat_t, breaks = 100)
hist(result_missForest_post$ximp, breaks = 100)

max(merged_normalized_log2ratio_6_mat_t, na.rm = TRUE) # 28.95969
max(result_missForest_post$ximp, na.rm = TRUE) # 28.95969

# check maximum value.
max(merged_analysis_before_normalization_ratio_mat_t, na.rm = TRUE) # 327053101
range(merged_analysis_before_normalization_ratio_mat_t, na.rm = TRUE) # 3.510291e-09 3.270531e+08

max(merged_log2ratio_6_mat_t, na.rm = TRUE) # 28.28495
2^max(merged_log2ratio_6_mat_t, na.rm = TRUE) # 327053101

# code 36

## go to line 872 to check the matrix "result_missForest_pre_argu_log2".

## Normalization again. Later, check if the estimated NA value is close to median value by CPTAC or submitter. 
## Remind; log2(merged_normalized_log2ratio_6_mat) - median = normalized_log2ratio

result_missForest_pre_argu_df_median_normalized <- sweep(x = result_missForest_pre_argu_log2,
                                            MARGIN = 2,
                                            STATS = merged_mean_median_sd_3_mat[2,],
                                            FUN = "-")

# view(result_missForest_pre_argu_df_median_normalized[1:100, 1:100]) # check data

## Make sure the value of data frame "result_missForest_pre_df_median_normalization" is numeric. 
# result_missForest_pre_argu_df_median_normalized <- apply(result_missForest_pre_argu_df_median_normalized, 2, as.numeric)
result_missForest_pre_argu_df_median_normalized <- apply(result_missForest_pre_argu_df_median_normalized, 2, as.numeric)

## Then the matrix "result_missForest_pre_df_median_normalization" transformed to data frame "result_missForest_pre_df_median_normalization_df"
# result_missForest_pre_argu_df_median_normalized_df <- data.frame(result_missForest_pre_argu_df_median_normalized)
result_missForest_pre_argu_df_median_normalized_df <- data.frame(result_missForest_pre_argu_df_median_normalized)

median(result_missForest_pre_argu_df_median_normalized) # 0.008310289, looks good because this dataset was centered with median, as the above.

## put column name (case) and row name(protein).
colnames(result_missForest_pre_argu_df_median_normalized_df) <- cases_pre_argu_normalization
rownames(result_missForest_pre_argu_df_median_normalized_df) <- proteins_pre_argu_normalization

## colnames(result_missForest_pre_argu_log2_df) <- cases_pre_argu_normalization
## rownames(result_missForest_pre_argu_log2_df) <- proteins_pre_argu_normalization

write_tsv(result_missForest_pre_argu_df_median_normalized_df, "/mnt/team4tb/Dropbox/Blog/2026 01 03 blog CPTAC/result_missForest_pre_argu_df_median_normalized_df.tsv")

cat(cases_pre_argu_normalization, file = "/mnt/team4tb/Dropbox/Blog/2026 01 03 blog CPTAC/cases_pre_argu_normalization.txt", sep = "\n")
cat(proteins_pre_argu_normalization, file = "/mnt/team4tb/Dropbox/Blog/2026 01 03 blog CPTAC/proteins_pre_argu_normalization.txt", sep = "\n")

# code 36

## go to line 872 to check the matrix "result_missForest_pre_argu_log2".

## Normalization again. Later, check if the estimated NA value is close to median value by CPTAC or submitter. 
## Remind; log2(merged_normalized_log2ratio_6_mat) - median = normalized_log2ratio

result_missForest_pre_argu_df_median_normalized <- sweep(x = result_missForest_pre_argu_log2,
                                            MARGIN = 2,
                                            STATS = merged_mean_median_sd_3_mat[2,],
                                            FUN = "-")

# view(result_missForest_pre_argu_df_median_normalized[1:100, 1:100]) # check data

## Make sure the value of data frame "result_missForest_pre_df_median_normalization" is numeric. 
# result_missForest_pre_argu_df_median_normalized <- apply(result_missForest_pre_argu_df_median_normalized, 2, as.numeric)
result_missForest_pre_argu_df_median_normalized <- apply(result_missForest_pre_argu_df_median_normalized, 2, as.numeric)

## Then the matrix "result_missForest_pre_df_median_normalization" transformed to data frame "result_missForest_pre_df_median_normalization_df"
# result_missForest_pre_argu_df_median_normalized_df <- data.frame(result_missForest_pre_argu_df_median_normalized)
result_missForest_pre_argu_df_median_normalized_df <- data.frame(result_missForest_pre_argu_df_median_normalized)

median(result_missForest_pre_argu_df_median_normalized) # 0.008310289, looks good because this dataset was centered with median, as the above.

## put column name (case) and row name(protein).
colnames(result_missForest_pre_argu_df_median_normalized_df) <- cases_pre_argu_normalization
rownames(result_missForest_pre_argu_df_median_normalized_df) <- proteins_pre_argu_normalization

## colnames(result_missForest_pre_argu_log2_df) <- cases_pre_argu_normalization
## rownames(result_missForest_pre_argu_log2_df) <- proteins_pre_argu_normalization

write_tsv(result_missForest_pre_argu_df_median_normalized_df, "/mnt/team4tb/Dropbox/Blog/2026 01 03 blog CPTAC/result_missForest_pre_argu_df_median_normalized_df.tsv")

cat(cases_pre_argu_normalization, file = "/mnt/team4tb/Dropbox/Blog/2026 01 03 blog CPTAC/cases_pre_argu_normalization.txt", sep = "\n")
cat(proteins_pre_argu_normalization, file = "/mnt/team4tb/Dropbox/Blog/2026 01 03 blog CPTAC/proteins_pre_argu_normalization.txt", sep = "\n")

# code 37

# save.image("/mnt/team4tb/Dropbox/Blog/2026 01 03 blog CPTAC/2026 01 10 CPTAC expression matrix_3.RData")

# code 37

# save.image("/mnt/team4tb/Dropbox/Blog/2026 01 03 blog CPTAC/2026 01 10 CPTAC expression matrix_3.RData")

# code 38

# Read clinical information.
## biospeciemen
biospeciemen <- read_tsv("/mnt/seqdata/public_data/Blog/CPTAC/PDC_biospecimen_manifest_10202024_102342.tsv")
nrow(biospeciemen) # 5228 cases

## clinical
clinical <- read_tsv("/mnt/seqdata/public_data/Blog/CPTAC/PDC_clinical_manifest_10202024_102355.tsv")
nrow(clinical) # 3563 cases

## followup
followup <- read_tsv("/mnt/seqdata/public_data/Blog/CPTAC/PDC_clinical_followup_manifest_10202024_102355.tsv")
nrow(followup) # 279 cases

## treatment
treatment <- read_tsv("/mnt/seqdata/public_data/Blog/CPTAC/PDC_clinical_treatment_manifest_10202024_102355.tsv")
nrow(treatment) # 84 cases

## exposure
exposure <- read_tsv("/mnt/seqdata/public_data/Blog/CPTAC/PDC_clinical_exposure_manifest_10202024_102355.tsv")
nrow(exposure) # 1222 cases

## study
study <- read_tsv("/mnt/seqdata/public_data/Blog/CPTAC/PDC_study_manifest_10202024_102331.tsv")
nrow(study) # 55 studies

## manifest file for download protein expression data
manifest <- read_tsv("/mnt/seqdata/public_data/Blog/CPTAC/PDC_file_manifest_10202024_102427.tsv")
nrow(manifest) # 278 studies

# code 38

# Read clinical information.
## biospeciemen
biospeciemen <- read_tsv("/mnt/seqdata/public_data/Blog/CPTAC/PDC_biospecimen_manifest_10202024_102342.tsv")
nrow(biospeciemen) # 5228 cases

## clinical
clinical <- read_tsv("/mnt/seqdata/public_data/Blog/CPTAC/PDC_clinical_manifest_10202024_102355.tsv")
nrow(clinical) # 3563 cases

## followup
followup <- read_tsv("/mnt/seqdata/public_data/Blog/CPTAC/PDC_clinical_followup_manifest_10202024_102355.tsv")
nrow(followup) # 279 cases

## treatment
treatment <- read_tsv("/mnt/seqdata/public_data/Blog/CPTAC/PDC_clinical_treatment_manifest_10202024_102355.tsv")
nrow(treatment) # 84 cases

## exposure
exposure <- read_tsv("/mnt/seqdata/public_data/Blog/CPTAC/PDC_clinical_exposure_manifest_10202024_102355.tsv")
nrow(exposure) # 1222 cases

## study
study <- read_tsv("/mnt/seqdata/public_data/Blog/CPTAC/PDC_study_manifest_10202024_102331.tsv")
nrow(study) # 55 studies

## manifest file for download protein expression data
manifest <- read_tsv("/mnt/seqdata/public_data/Blog/CPTAC/PDC_file_manifest_10202024_102427.tsv")
nrow(manifest) # 278 studies

## code 39
### Remove unnecessary characters from the columns of clinical information. 

print("------------------------------biospeciemen----------------------------------------------")
colnames(biospeciemen)[1:10]
print("------------------------------treatment----------------------------------------------")
colnames(treatment)[1:10]
print("------------------------------exposure----------------------------------------------")
colnames(exposure)[1:10]
print("------------------------------followup----------------------------------------------")
colnames(followup)[1:10]
print("------------------------------clinical----------------------------------------------")
colnames(clinical)[1:10]

### biospeciemen
### "biospeciemen_renamed" will be used for matching the file name of protein expression later.

colname_biospeciemen <- colnames(biospeciemen)
colname_biospeciemen <- gsub(pattern = " ", "_", colname_biospeciemen)
colnames(biospeciemen) <- colname_biospeciemen

biospeciemen$Aliquot_Submitter_ID_2 <- gsub(pattern = "-", replacement = "_", x = biospeciemen$Aliquot_Submitter_ID)
biospeciemen <- biospeciemen %>% select(c(1,2, ncol(biospeciemen), 3:(ncol(biospeciemen)-1)))
ncol(biospeciemen)

### treatment
colname_treatment <- colnames(treatment)
colname_treatment <- gsub(pattern = " ", "_", colname_treatment)
colnames(treatment) <- colname_treatment

### exposure
colname_exposure <- colnames(exposure)
colname_exposure <- gsub(pattern = " ", "_", colname_exposure)
colnames(exposure) <- colname_exposure

### followup
colname_followup <- colnames(followup)
colname_followup <- gsub(pattern = " ", "_", colname_followup)
colnames(followup) <- colname_followup

### rename column "Case_Id", "Case_Submitter_Id", "Follow_Up_Submitter_Id", "Follow_Up_Id" for merge with other data frame.
followup <- followup %>% dplyr::rename("Case_ID" = "Case_Id") 
followup <- followup %>% dplyr::rename("Case_Submitter_ID" = "Case_Submitter_Id") 
followup <- followup %>% dplyr::rename("Follow_Up_Submitter_ID" = "Follow_Up_Submitter_Id") 
followup <- followup %>% dplyr::rename("Follow_Up_ID" = "Follow_Up_Id") 

### clinical
colname_clinical <- colnames(clinical)
colname_clinical <- gsub(pattern = " ", "_", colname_clinical)
colnames(clinical) <- colname_clinical

## study
colname_study <- colnames(study)
colname_study <- gsub(pattern = " ", "_", colname_study)
colnames(study) <- colname_study
study <- study %>% dplyr::rename("Case_number" = "Cases_#") # rename column "Case_Id" to "Case_ID" for merge with other data frame.

clinical_info_list <- list(treatment, exposure, followup, clinical)
names(clinical_info_list) <- c("treatment", "exposure", "followup", "clinical")

# Retrieve column name again. Becuase the "_" was removed from the columns at the above.
colname_biospeciemen_df <- data.frame(column = colnames(biospeciemen))
colname_clinical_df <- data.frame(column = colnames(clinical))
colname_treatment_df <- data.frame(column = colnames(treatment))
colname_exposure_df <- data.frame(column = colnames(exposure))
colname_followup_df <- data.frame(column = colnames(followup))
colname_study_df <- data.frame(column = colnames(study))

# this data frame will be used for checking the data. 
clinical_info_merged <- full_join(biospeciemen, clinical, by = "Case_ID", suffix = c("", "_clinical"))
nrow(clinical_info_merged) # 5306

clinical_info_merged <- full_join(clinical_info_merged, followup, by = "Case_ID", suffix = c("", "_followup"))
nrow(clinical_info_merged) # 5306

clinical_info_merged <- full_join(clinical_info_merged, treatment, by = "Case_ID", suffix = c("", "_treatment"))
nrow(clinical_info_merged) # 5339

clinical_info_merged <- full_join(clinical_info_merged, exposure, by = "Case_ID", suffix = c("", "_exposure"))
nrow(clinical_info_merged) # 5339

## code 39
### Remove unnecessary characters from the columns of clinical information. 

print("------------------------------biospeciemen----------------------------------------------")
colnames(biospeciemen)[1:10]
print("------------------------------treatment----------------------------------------------")
colnames(treatment)[1:10]
print("------------------------------exposure----------------------------------------------")
colnames(exposure)[1:10]
print("------------------------------followup----------------------------------------------")
colnames(followup)[1:10]
print("------------------------------clinical----------------------------------------------")
colnames(clinical)[1:10]

### biospeciemen
### "biospeciemen_renamed" will be used for matching the file name of protein expression later.

colname_biospeciemen <- colnames(biospeciemen)
colname_biospeciemen <- gsub(pattern = " ", "_", colname_biospeciemen)
colnames(biospeciemen) <- colname_biospeciemen

biospeciemen$Aliquot_Submitter_ID_2 <- gsub(pattern = "-", replacement = "_", x = biospeciemen$Aliquot_Submitter_ID)
biospeciemen <- biospeciemen %>% select(c(1,2, ncol(biospeciemen), 3:(ncol(biospeciemen)-1)))
ncol(biospeciemen)

### treatment
colname_treatment <- colnames(treatment)
colname_treatment <- gsub(pattern = " ", "_", colname_treatment)
colnames(treatment) <- colname_treatment

### exposure
colname_exposure <- colnames(exposure)
colname_exposure <- gsub(pattern = " ", "_", colname_exposure)
colnames(exposure) <- colname_exposure

### followup
colname_followup <- colnames(followup)
colname_followup <- gsub(pattern = " ", "_", colname_followup)
colnames(followup) <- colname_followup

### rename column "Case_Id", "Case_Submitter_Id", "Follow_Up_Submitter_Id", "Follow_Up_Id" for merge with other data frame.
followup <- followup %>% dplyr::rename("Case_ID" = "Case_Id") 
followup <- followup %>% dplyr::rename("Case_Submitter_ID" = "Case_Submitter_Id") 
followup <- followup %>% dplyr::rename("Follow_Up_Submitter_ID" = "Follow_Up_Submitter_Id") 
followup <- followup %>% dplyr::rename("Follow_Up_ID" = "Follow_Up_Id") 

### clinical
colname_clinical <- colnames(clinical)
colname_clinical <- gsub(pattern = " ", "_", colname_clinical)
colnames(clinical) <- colname_clinical

## study
colname_study <- colnames(study)
colname_study <- gsub(pattern = " ", "_", colname_study)
colnames(study) <- colname_study
study <- study %>% dplyr::rename("Case_number" = "Cases_#") # rename column "Case_Id" to "Case_ID" for merge with other data frame.

clinical_info_list <- list(treatment, exposure, followup, clinical)
names(clinical_info_list) <- c("treatment", "exposure", "followup", "clinical")

# Retrieve column name again. Becuase the "_" was removed from the columns at the above.
colname_biospeciemen_df <- data.frame(column = colnames(biospeciemen))
colname_clinical_df <- data.frame(column = colnames(clinical))
colname_treatment_df <- data.frame(column = colnames(treatment))
colname_exposure_df <- data.frame(column = colnames(exposure))
colname_followup_df <- data.frame(column = colnames(followup))
colname_study_df <- data.frame(column = colnames(study))

# this data frame will be used for checking the data. 
clinical_info_merged <- full_join(biospeciemen, clinical, by = "Case_ID", suffix = c("", "_clinical"))
nrow(clinical_info_merged) # 5306

clinical_info_merged <- full_join(clinical_info_merged, followup, by = "Case_ID", suffix = c("", "_followup"))
nrow(clinical_info_merged) # 5306

clinical_info_merged <- full_join(clinical_info_merged, treatment, by = "Case_ID", suffix = c("", "_treatment"))
nrow(clinical_info_merged) # 5339

clinical_info_merged <- full_join(clinical_info_merged, exposure, by = "Case_ID", suffix = c("", "_exposure"))
nrow(clinical_info_merged) # 5339