Add initial R project file for xMap biomarker analysis

xmap_biomarkers/scripts/backend/data_handling.R

@@ -0,0 +1,118 @@
+#fill out useful info
+fill_analyte_info <- function(df, antigens=I00_Antigens){
+  lookup_df <- bind_rows(antigens)[-1]
+  # Remove duplicates from lookup_df based on Antigen.name
+  lookup_df <- lookup_df %>%
+    distinct(Antigen.name, .keep_all = TRUE)
+  # Create a new column with the row names in main_df
+  main_df <- df %>%
+    rownames_to_column(var = "RowName")
+  # Check if most RowName values are not found in Antigen.name
+  rowname_matches <- main_df$RowName %in% lookup_df$Antigen.name
+  proportion_matches <- sum(rowname_matches) / length(main_df$RowName)
+  
+  if (proportion_matches < 0.5) {
+    warning("Most RowName values are not found in Antigen.name")
+  } else {
+    message("Most RowName values are found in Antigen.name")
+  }
+  # Merge the main dataframe with the lookup dataframe based on the Antigen.name column
+  result_df <- main_df %>%
+    left_join(lookup_df, by = c("RowName" = "Antigen.name"))
+  return(result_df)
+}
+
+# Function to subset the dataframe based on the lowest (n) values
+subset_top_n <- function(df, n) {
+  if ("Q.Value" %in% colnames(df)) {
+    subset_df <- df %>% top_n(-n, Q.Value)
+  } else if ("adj.P.Val" %in% colnames(df)) {
+    subset_df <- df %>% top_n(-n, adj.P.Val)
+  } else {
+    return(df)
+  }
+  #subset_df <- rownames_to_column(subset_df, var = "RowName")
+  return(subset_df)
+}
+
+excel_subset_export <- function(dflist, n=15){
+  # Apply the function to the dflist
+  subsetted_dflist <- lapply(dflist, subset_top_n, n)
+  subsetted_dflist <- lapply(subsetted_dflist,fill_analyte_info)
+  # Get the name of the dflist
+  dflist_name <- deparse(substitute(dflist))
+  # Write the list of subsetted dataframes to an Excel file
+  excel_file_name <- paste0(dflist_name, "_Top_", n, ".xlsx")
+  # Create a new workbook
+  wb <- createWorkbook()
+  # Add worksheets for each subsetted dataframe and write data
+  for (i in seq_along(subsetted_dflist)) {
+    sheet_name <- names(dflist)[i]
+    addWorksheet(wb, sheet_name)
+    writeData(wb, sheet_name, subsetted_dflist[[i]])
+  }
+  
+  # Save the workbook
+  saveWorkbook(wb, file = file.path("stats", excel_file_name), overwrite = TRUE)
+}
+
+#Final Modifications
+full_analyte_info <- function(df, antigens=I00_Antigens){
+  lookup_df <- bind_rows(antigens)[-1]
+  # Remove duplicates from lookup_df based on Antigen.name
+  lookup_df <- lookup_df %>%
+    distinct(Antigen.name, .keep_all = TRUE)
+  # Create a new column with the row names in main_df
+  main_df <- df %>%
+    rownames_to_column(var = "RowName")
+  # Check if most RowName values are not found in Antigen.name
+  rowname_matches <- main_df$RowName %in% lookup_df$Antigen.name
+  proportion_matches <- sum(rowname_matches) / length(main_df$RowName)
+  
+  if (proportion_matches < 0.5) {
+    warning("Most RowName values are not found in Antigen.name")
+  }
+  # else {
+  #   message("Most RowName values are found in Antigen.name")
+  # }
+  # Merge the main dataframe with the lookup dataframe based on the Antigen.name column
+  result_df <- main_df %>%
+    left_join(lookup_df, by = c("RowName" = "Antigen.name"))
+  return(result_df)
+}
+
+excel_export <- function(dflist, name = deparse(substitute(dflist))) {
+  file_path <- file.path(getwd(), "reports", paste0(name, ".xlsx"))
+  
+  # Create a new workbook
+  wb <- createWorkbook()
+  
+  # Iterate through the list of dataframes and add worksheets
+  for (setid in seq_along(dflist)) {
+    sheet_name <- paste0("Dataset_", setid)
+    addWorksheet(wb, sheet_name)
+    writeData(wb, sheet_name, dflist[[setid]])
+  }
+  
+  # Save the workbook
+  saveWorkbook(wb, file = file_path, overwrite = TRUE)
+}
+
+replace_antigen_names<- function(vector, antigens=I00_Antigens, replacement_col = "Gene.name"){
+  lookup_df <- bind_rows(antigens)[-1]
+  # Remove duplicates from lookup_df based on Antigen.name
+  lookup_df <- lookup_df %>%
+    distinct(Antigen.name, .keep_all = TRUE)
+  
+  # Create a lookup dictionary with Antigen.name as the key and the specified column as the value
+  lookup_dict <- setNames(lookup_df[[replacement_col]], lookup_df$Antigen.name)
+  replaced_vector <- lookup_dict[vector]
+  return(replaced_vector)
+}
+
+untransform_subset <- function(restriction_df, df, subset_method = limma_subset, ...){
+  subset <- do.call(subset_method, c(list(restriction_df), list(...)))
+  df_subset <- df[,colnames(df) %in% colnames(subset)]
+  colnames(df_subset)[-1:-2]=replace_antigen_names(colnames(df_subset[-1:-2]))
+  return(df_subset)
+}

xmap_biomarkers/scripts/backend/data_import.r

@@ -0,0 +1,148 @@
+#Takes a list of package names and automates the process of checking, installing, and loading them. 
+#Handles both CRAN and Bioconductor packages, ensuring that the necessary packages are available in the user's R environment.
+load_dependencies <- function(pkg_list) {
+  # Load or install packages from list
+  for (pkg in pkg_list) {
+    if (substr(pkg, 1, 11) == "BiocManager") {
+      pkg = substr(pkg, 14, nchar(pkg))
+      if (!requireNamespace("BiocManager", quietly = TRUE)) {
+        install.packages("BiocManager")
+      }
+      if (!require(pkg, character.only = TRUE)) {
+        BiocManager::install(pkg)
+      }
+    }
+    else{
+      if (!require(pkg, character.only = TRUE)) {
+        install.packages(pkg)
+      }
+    }
+    library(pkg, character.only = TRUE)
+  }
+}
+# This function checks the integrity of a set of dataset files in a specified directory, filtering them by keyword, 
+# and ensures that each Data Intensity file has a corresponding Antigen List file. 
+# If any files are found to be missing an Antigen List, the function returns a list of these files, otherwise it returns TRUE.
+  dataset_file_integrity_check <- function(keyword) {
+    # Get the list of files in the "sample_data" subdirectory
+    file_list <- list.files("data", pattern = "\\.xlsx", full.names = TRUE)
+    
+    # Filter the file list to include only files containing the specified keyword
+    keyword_files <- grep(keyword, file_list, value = TRUE, ignore.case = TRUE)
+    
+    # Replace any space or hyphen in each of the names in keyword files with an underscore
+    keyword_files <- gsub(" |-", "_", keyword_files)
+    
+    # Remove characters in the file names up to one after the dataset keyword
+    keyword_files <- gsub(paste0(".*", keyword, "[ _-](.*)\\.xlsx$"), "\\1.xlsx", keyword_files)
+    
+    # Find the non-dictionary tags present in the names of several files
+    non_dict_tags <- unique(gsub("(_Antigen_list.xlsx|_Data_Intensity.xlsx)", "", keyword_files))
+    
+    # Identify Data Intensity files that don't have a corresponding Antigen List file
+    missing_antigen_list_files <- character()
+    for (tag in non_dict_tags) {
+      antigen_file <- paste0(tag, "_Antigen_list.xlsx")
+      intensity_file <- paste0(tag, "_Data_Intensity.xlsx")
+      if (intensity_file %in% keyword_files && !(antigen_file %in% keyword_files)) {
+        missing_antigen_list_files <- c(missing_antigen_list_files, intensity_file)
+      }
+    }
+    
+    # Check if the list of data files with no antigen file is empty
+    if (length(missing_antigen_list_files) == 0) {
+      return(TRUE)
+    } else {
+      #print(paste("Some files are missing antigen lists:", missing_antigen_list_files))
+      return(missing_antigen_list_files)
+    }
+  }
+# Function to filter out data frames with a given keyword in their name
+# keyword: the keyword to search for in the data frame names
+# e: the environment to search in (default is parent frame)
+  filterclean <- function(keyword, e = parent.frame()) {
+  
+  # Get list of data frames in the specified environment
+  dflist = Filter(function(x) is (x, "data.frame"),
+                  mget(ls(e),envir= e))
+  # print(dflist)
+  # Filter out data frames without the specified keyword
+  dflist = dflist[grepl(keyword,ls(dflist))]
+  # print(dflist)
+  # Return filtered list of data frames
+  return (dflist)
+}
+## This import function first brings all excel documents in the sample_data directory which include the dataset string and imports them as dataframes
+#  The dataframes are then grouped according to _Data and _Antigen keywords in the file naming convention (filterclean function)
+  import <- function(dataset){
+    ## Import all libraries needed for downstream
+    pkg_list = c("rlang",
+                 "gridExtra", 
+                 "ggplot2", 
+                 "ggfortify", 
+                 "MASS", 
+                 "BiocManager::lumi",
+                 "BiocManager::limma", 
+                 "readxl", 
+                 "dplyr",
+                 "broom",
+                 "BiocManager::qvalue",
+                 "openxlsx",
+                 "tibble",
+                 "pheatmap",
+                 "pROC",
+                 "tidyverse",
+                 "msigdbr",
+                 "BiocManager::clusterProfiler")
+    load_dependencies(pkg_list)
+    # set path to sample data directory
+    sample_data=paste(getwd(),"/data/", sep="")
+    # get names of all files with .xlsx extension in sample_data directory
+    names = list.files(path=sample_data, pattern = ".xlsx", recursive=TRUE)
+    
+    # check if any required files are missing
+    no_missing_files <- dataset_file_integrity_check(dataset)
+    
+    # if files are missing
+    if (no_missing_files != TRUE) {
+      # if files are missing
+      # ask user if they want to halt the function or continue without missing files
+      stop_message <- paste0("Some files are missing antigen lists: ", no_missing_files, "\n")
+      cat(stop_message)
+      choice <- readline("Do you want to halt the function? (y/n) ")
+      # if user chooses to halt the function, stop and print error message
+      if (choice == "y") {
+        stop(stop_message)
+      }}
+    # read all files in sample_data directory and assign to variables with shortened names
+    for (file in names){
+      shortindex = gregexpr(pattern=dataset,file)[[1]][1]
+      shortname = substring(file,shortindex+6)
+      shortname = gsub(" |-", "_", shortname)
+      assign(paste0(shortname), read_xlsx(paste(sample_data,file,sep="")))
+    }
+    # assign cleaned data and antigen data to variables with dataset name and appropriate suffixes
+    assign(paste(dataset, "_A00", sep=""),filterclean("_Data"))
+    assign(paste(dataset, "_Antigens", sep=""), filterclean("_Antigen"))
+    # clean and format data
+    A01_Input = lapply(get(paste(dataset, "_A00", sep="")),function(df){
+      df = data.frame(df)
+      names(df)[2] = "group"
+      # assign group value of 1 if it contains "GC", 0 if it contains "HD", and "NA" if neither
+      df$group =ifelse(grepl("GC",df$group), 1,
+                       ifelse(grepl("HD",df$group), 0, "NA"))
+      df
+    })
+    # format antigen data
+    AA_Antigens = lapply(get(paste(dataset, "_Antigens", sep="")),function(df){
+      df = data.frame(df)
+      names(df)[1] = "analyte"
+      df
+    })
+    # assign cleaned and formatted data to global environment variables
+    assign("I01_Import", A01_Input, env=globalenv())
+    assign("I00_Antigens", AA_Antigens, env=globalenv())
+  }
+#   
+  
+

xmap_biomarkers/scripts/backend/display_division.R

@@ -0,0 +1,54 @@
+# These functions are used to calculate the optimal number of columns and rows for a display of n items. 
+# The display_division function takes an argument n 
+# which represents the total number of items to be displayed.
+# It then calls two other functions h_display_division and v_display_division 
+# to calculate the optimal number of columns (h) and rows (v) for displaying the items.
+
+h_display_division <- function(n, max_div = 5){
+  divisors = 4:max_div
+  if (n>max_div){
+    valid_divisors = divisors[n %% divisors ==0]
+    if (length(valid_divisors) > 0){
+      return(max(valid_divisors))
+    }
+    else{
+      remains = list()
+      for(i in seq_along(divisors)){
+        remains[i] = n %% divisors[i]}
+      return(divisors[length(divisors)-which.max(rev(remains))+1])
+    }
+  }
+  else{
+    return(n)
+  }
+}
+
+v_display_division <- function(n, h, max_div = 4){
+  divisors = 3:max_div
+  if (n/h< max_div){
+    return(ceiling(n/h))
+  }
+  else {
+    valid_divisors = divisors[n %% divisors*h ==0]
+    if (length(valid_divisors) > 0){
+      return(max(valid_divisors))
+    }
+    else{
+      remains = list()
+      for(i in seq_along(divisors)){
+        remains[i] = n %% divisors[i]*h}
+      return(divisors[length(divisors)-which.max(rev(remains))+1])
+    }
+  }
+}
+
+
+
+display_division <- function(n){
+  h = h_display_division(n)
+  v = v_display_division(n,h)
+  d = ceiling(n/(h*v))
+  return(c(h,v,d))
+}
+
+

xmap_biomarkers/scripts/backend/fancy_roc.R

@@ -0,0 +1,81 @@
+fancy_roc <- function(df, report, P.Val = 0.05, validation = "none") {
+  library(pROC)
+  library(ggplot2)
+  lrep_sigs <- report[report$adj.P.Val < P.Val,]
+  analytes <- row.names(lrep_sigs)
+  df_selected <- df[, c("group", analytes)]
+  df_selected$group = as.numeric(df_selected$group)
+  logistic_model <- glm(group ~ ., data = df_selected, family = "binomial")
+  
+  logistic_model_stepwise <- NULL
+  auc_stepwise <- NULL
+  predicted_probabilities_stepwise <- NULL
+  if (length(analytes) > 1) {
+    logistic_model_stepwise <- step(logistic_model, direction = "backward")
+    predicted_probabilities_stepwise <- predict(logistic_model_stepwise, type = "response")
+    roc_obj_stepwise <- roc(df_selected$group, predicted_probabilities_stepwise)
+    auc_stepwise <- auc(roc_obj_stepwise)
+  }
+  
+  predicted_probabilities <- predict(logistic_model, type = "response")
+  roc_obj <- roc(df_selected$group, predicted_probabilities)
+  auc <- auc(roc_obj)
+  
+  roc_data <- data.frame(
+    FPR = roc_obj$specificities,
+    TPR = roc_obj$sensitivities,
+    Model = "Logistic Regression"
+  )
+  
+  if (!is.null(logistic_model_stepwise)) {
+    roc_data_stepwise <- data.frame(
+      FPR = roc_obj_stepwise$specificities,
+      TPR = roc_obj_stepwise$sensitivities,
+      Model = "Logistic Regression (Backwards Step)"
+    )
+    roc_data <- rbind(roc_data, roc_data_stepwise)
+  }
+  
+    plot <- ggplot(data = roc_data, aes(x = FPR, y = TPR, color = Model)) +
+      geom_line(size = 1) +
+      labs(
+        x = "1 - Specificity",
+        y = "Sensitivity",
+        title = ""
+      ) 
+    
+    +
+      theme(
+        panel.background = element_rect("white"),
+        legend.position = c(0.75,0.15), # Change legend location
+        #legend.justification = "bottom", 
+        plot.title = element_text(hjust = 0.5),
+        legend.text = element_text(size=14), # Change legend text size
+        legend.background = element_rect(color="black"),
+        axis.text = element_text(size = 14) # Change axis tick text size
+      ) +
+      coord_cartesian(xlim = c(1, 0), ylim = c(0, 1)) +
+      scale_color_manual(values = c("red", "blue")) +
+      annotate(
+        "text",
+        x = 0.4,
+        y = 0.3,
+        label = paste("AUC: ", paste0(round(auc, 3))),
+        color = "red",
+        size = 6 #AUC text size
+      )
+    
+    if (!is.null(logistic_model_stepwise)) {
+      plot <- plot +
+        annotate(
+          "text",
+          x = 0.4,
+          y = 0.25,
+          label = paste("AUC: ", paste0(round(auc_stepwise, 3))),
+          color = "blue",
+          size = 6 #AUC text size
+        )
+    }
+    
+    print(plot)
+}

xmap_biomarkers/scripts/backend/gsea.R

@@ -0,0 +1,33 @@
+msigdb_workflow <- function(signif.genes, category = "C2") {
+  library(msigdbr)
+  library(ggplot2)
+
+  msigdb_data <- msigdbr(species = "Homo sapiens", category = category)
+  # signif.genes <- reports_P03_Transformed_bcrsn$Set_3_limma
+  # top_genes <- head(signif.genes[order(signif.genes$adj.P.Val),], 25)
+  top_genes <- signif.genes[order(signif.genes$adj.P.Val), ]
+  top_genes <- fill_analyte_info(top_genes)
+  top_genes <- top_genes$Gene.name
+
+  msigdb_genes <- select(msigdb_data, gs_name, gene_symbol)
+  enrich_msigdb <- enricher(top_genes, TERM2GENE = msigdb_genes)
+
+  enrich_msigdb_df <- enrich_msigdb@result %>%
+    separate(BgRatio, into = c("size.term", "size.category"), sep = "/") %>%
+    separate(GeneRatio, into = c("size.overlap.term", "size.overlap.category"), sep = "/") %>%
+    mutate_at(vars("size.term", "size.category", "size.overlap.term", "size.overlap.category"), as.numeric) %>%
+    mutate("k.K" = size.overlap.term / size.term)
+
+  enrich_plot <- enrich_msigdb_df %>%
+    filter(p.adjust <= 0.05) %>%
+    ggplot(aes(x = reorder(Description, k.K), y = k.K)) +
+    geom_col() +
+    theme_classic() +
+    coord_flip() +
+    labs(
+      y = "Significant genes in set / Total genes in set \nk/K", x = "Gene set",
+      title = paste("Differentially expressed genes enriched in", category, "Gene sets (P<0.05)")
+    )
+
+  return(list(data = enrich_msigdb_df, plot = enrich_plot))
+}

xmap_biomarkers/scripts/backend/modeling.R

@@ -0,0 +1,95 @@
+loocv_validation <- function(df, logistic_model) {
+  n_samples <- nrow(df)
+  predicted_probabilities <- numeric(n_samples)
+  
+  for (i in 1:n_samples) {
+    test_set <- df[i,]
+    test_prob <- predict(logistic_model, newdata = test_set[-1], type = "response")
+    predicted_probabilities[i] <- test_prob
+  }
+  
+  return(predicted_probabilities)
+  
+}
+
+roc_curve <- function(df, report, P.Val = 0.05, stepwise = FALSE, validation = "none") {
+  lrep_sigs <- report[report$adj.P.Val < P.Val,]
+  analytes <- row.names(lrep_sigs)
+  df_selected <- df[, c("group", analytes)]
+  df_selected$group = as.numeric(df_selected$group)
+  logistic_model <- glm(group ~ ., data = df_selected, family = "binomial")
+  
+  if (stepwise) {
+    logistic_model <- step(logistic_model, direction = "backward")
+  }
+  
+  model_summary <- summary(logistic_model)
+  accuracy <- NULL
+  if (validation == "LOOCV") {
+    predicted_probabilities <- loocv_validation(df, logistic_model)
+    true_labels <- as.numeric(df$group)
+    threshold <- 0.5
+    predicted_labels <- ifelse(predicted_probabilities >= threshold, 1, 0)
+    correct_predictions <- predicted_labels == true_labels
+    accuracy <- mean(correct_predictions)
+    roc_obj <- roc(true_labels, predicted_probabilities)
+  } else {
+    predicted_probabilities <- predict(logistic_model, type = "response")
+    roc_obj <- roc(df_selected$group, predicted_probabilities)
+  }
+  
+  auc <- auc(roc_obj)
+  plot(roc_obj, main = paste("ROC Curve (AUC =", round(auc, 3), ")"),asp=1)
+  roc_plot <- recordPlot()
+  
+  logistic_regression <- list(model_summary = model_summary, logistic_model = logistic_model)
+  validation_results <- list(validation_method = validation, predicted_probabilities = predicted_probabilities, accuracy = accuracy)
+  roc_results <- list(roc_obj = roc_obj, auc = auc, plot = roc_plot)
+  
+  if (validation == "none") {
+    results <- list(logistic_regression = logistic_regression, roc = roc_results)
+  } else {
+    results <- list(logistic_regression = logistic_regression, validation = validation_results, roc = roc_results)
+  }
+  
+  return(results)
+}
+
+fancy_roc <- function(df, report, P.Val = 0.05, validation = "none") {
+  library(pROC)
+  lrep_sigs <- report[report$adj.P.Val < P.Val,]
+  analytes <- row.names(lrep_sigs)
+  df_selected <- df[, c("group", analytes)]
+  df_selected$group = as.numeric(df_selected$group)
+  logistic_model <- glm(group ~ ., data = df_selected, family = "binomial")
+  
+  logistic_model_stepwise <- NULL
+  auc_stepwise <- NULL
+  predicted_probabilities_stepwise <- NULL
+  if (length(analytes) > 1) {
+    logistic_model_stepwise <- step(logistic_model, direction = "backward")
+    predicted_probabilities_stepwise <- predict(logistic_model_stepwise, type = "response")
+    roc_obj_stepwise <- roc(df_selected$group, predicted_probabilities_stepwise)
+    auc_stepwise <- auc(roc_obj_stepwise)
+  }
+  
+  predicted_probabilities <- predict(logistic_model, type = "response")
+  roc_obj <- roc(df_selected$group, predicted_probabilities)
+  auc <- auc(roc_obj)
+  
+  par(cex.axis = 1.5)
+  plot(roc_obj, col="red", main = "", xlim=c(1,0), ylim=c(0,1))
+
+  if (!is.null(logistic_model_stepwise)) {
+    predicted_probabilities_stepwise <- predict(logistic_model_stepwise, type = "response")
+    roc_obj_stepwise <- roc(df_selected$group, predicted_probabilities_stepwise)
+    auc_stepwise <- auc(roc_obj_stepwise)
+    lines(roc_obj_stepwise, col="blue")
+    legend("bottomright", legend=c("Logistic Regression", "Logistic Regression (Backwards Step)"), col=c("red", "blue"), lty=1, cex=0.8)
+  }
+  else {
+    legend("bottomright", legend="Logistic Regression", col="red", lty=1, cex=0.8)
+  }
+  text(x=0.4, y=0.4, labels=paste("AUC: ", paste0(round(auc, 3))), pos=1, cex=1, col="red")
+  text(x=0.4, y=0.35, labels=paste("AUC: ", paste0(round(auc_stepwise, 3))), pos=1, cex=1, col="blue")}
+

xmap_biomarkers/scripts/backend/preprocessing_stage_1.R

@@ -0,0 +1,144 @@
+
+detect_background_noise <- function(set, empty_id = NULL){
+  # Get list of unique Internal.LIMS.ID values containing the word "empty"
+  empty_ids <- unique(grep("empty", set$Internal.LIMS.ID, value = TRUE, ignore.case = TRUE))
+  
+  # If empty_id is not specified by the user, prompt the user to choose from the available options
+  if (is.null(empty_id)){
+    if (length(empty_ids) == 0){
+      stop("No empty samples found in dataset")
+    } else if (length(empty_ids) == 1){
+      empty_id <- empty_ids
+      message(paste0("Using empty ID: ", empty_id))
+    } else {
+      message("Multiple empty IDs found in dataset:")
+      for (i in seq_along(empty_ids)){
+        message(paste0(i, ": ", empty_ids[i]))
+      }
+      empty_id <- readline(prompt = "Enter the number corresponding to the desired empty ID: ")
+      if (!as.numeric(empty_id) %in% seq_along(empty_ids)){
+        stop("Invalid input. Aborting.")
+      } else {
+        empty_id <- empty_ids[as.numeric(empty_id)]
+      }
+    }
+  }
+  # Filter the data frame to include only the chosen empty ID
+  emptyset <- set[set$Internal.LIMS.ID == empty_id, ]
+  
+  # Calculate the background noise
+  inset = emptyset[-1:-2][colMeans(emptyset[-1:-2])<(median(colMeans(emptyset[-1:-2]))+(1*sd(colMeans(emptyset[-1:-2]))))]
+  if (length(inset) < 0.95*length(emptyset)){
+    calset = emptyset[-1:-2]
+  }else{
+    calset = inset
+  }
+  return(max(colMeans(calset))+sd(colMeans(calset)))
+}
+
+purge_background <- function(set, cutoff){
+  bgmap = set[-1:-2]> cutoff
+  above_cutoff = set[-1:-2][,colSums(bgmap)>0]
+  below_cutoff = set[-1:-2][,colSums(bgmap)==0]
+  keep_set = cbind(set[1:2],above_cutoff)
+  remove_set = cbind(set[1:2],below_cutoff)
+  return(list(keep_set, remove_set))
+}
+
+# The handle_background function takes a list of dataframes dflist,
+# applies the purge_background function to each dataframe using detect_background_noise to determine the cutoff, 
+# and replaces the original dataframe with the keep_set. 
+# The remove_set is added to a new dflist in the global environment called X01_bg_purge. 
+# The function returns the new dflist.
+handle_background <- function(dflist) {
+  new_dflist <- list()
+  scrap_dflist <-list()
+  for (i in seq_along(dflist)) {
+    df <- dflist[[i]]
+    cutoff <- detect_background_noise(df)
+    keep_set <- purge_background(df, cutoff)[[1]]
+    remove_set <- purge_background(df, cutoff)[[2]]
+    new_dflist[[i]] <- keep_set
+    scrap_dflist[[i]] <- remove_set
+  }
+  names(new_dflist) <- names(dflist)
+  names(scrap_dflist) <- names(dflist)
+  assign("X01_bg_purge", scrap_dflist, envir = .GlobalEnv)
+  return(new_dflist)
+}
+
+# This function takes a dataframe "set" and adjusts it by subtracting the average value of an "EMPTY-0001" sample from all other samples. 
+# It then sets any values less than 0 to 0, and adds 1 to all values to avoid breaking the log() function. The adjusted dataframe is returned.
+emptyadjust <- function(set){
+  emptyset = set[set$Internal.LIMS.ID=="EMPTY-0001",]
+  emptyvector = colMeans(emptyset[-1:-2])
+  set = set[set$Internal.LIMS.ID!="EMPTY-0001" & set$Internal.LIMS.ID!="MIX_2-0029",]
+  labels = set[1:2]
+  set = cbind(set[1:2],sweep(set[-1:-2],2,FUN="-",emptyvector))
+  ### This sets anything with less read than the empty to 0, then adds one to everything, so that it doesn't break the log()
+  set[-1:-2][set[-1:-2]<0] <- 0
+  set[-1:-2] = set[-1:-2]+1
+  return(set)
+}
+
+# The compress_duplicates function takes a dataframe "set" and an excel file "layout" as input
+# and combines rows in the dataframe based on a shared identifier in the "layout" file. 
+# It identifies matching rows, averages their values, and keeps the first row while removing the rest. 
+# If any technical replicates have a deviation greater than either item, it prints notifications.
+compress_duplicates <- function(set, layout){
+  outlist = c()
+  ### Reads the layout from excel file
+  slayout <- read_excel(layout)
+  ## Looks through the column named Tube Label for any items containing a hyphen, then uses any name preceding a hyphen as the for item
+  for (n in strsplit(slayout$`Tube label`[grepl("-",slayout$`Tube label`)],"-")){
+    # Filters the layout to only hold items either ending in or containing the for item immediately before the hyphen
+    mergerows = ((slayout %>% filter(grepl(paste0(n[1],"$|",n[1],"-"),`Tube label`)))$`Sample id_LIMS`)
+    # Returns the sample id for those samples, which is present in the main data set
+    mergeset = set[grepl(paste(mergerows, collapse = "|"), set$Internal.LIMS.ID),]
+    # Checks that both technical replicates are close enough together that their deviation is not greater than either item (like if one was 1000 and one was 10, this would print notifications)
+    outlierflag = colSums(sweep(mergeset[-1:-2],2,apply(mergeset[-1:-2], 2, sd ), '-')<0)
+    if (sum(outlierflag)>0){
+      #print(mergerows)
+      #print(outlierflag[outlierflag>0])
+      #tupsum = c(mergerows, colnames(outlierflag[outlierflag>0]))
+      #print(tupsum)
+    }
+    # Reassigns the first row in the set of matches so that it is equal to the mean of all matching sets, for each variable
+    set[grepl(paste(mergerows, collapse = "|"), set$Internal.LIMS.ID),][1,][-1:-2] = colMeans(set[grepl(paste(mergerows, collapse = "|"), set$Internal.LIMS.ID),][-1:-2])
+    # Eliminates all matching samples except the first (which is now the average of all matching) from the dataset
+    for (i in 2:length(mergerows)){
+      set = set[row.names(set) != row.names(set[grepl(paste(mergerows, collapse = "|"), set$Internal.LIMS.ID),][i,]),]
+    }}
+  return(set)
+}
+
+# This function takes a list of data frames dflist and a list of antigen names antigens. 
+# It loops over each data frame in dflist and renames the column names of the data frames according to the antigen list. 
+# If mode=1, the column names are set to the Antigen name 
+# If mode=2, the column names are set to the Gene name. 
+# The function then returns a list of data frames with updated column names.
+set_colname_adapter <- function (dflist, antigens = I00_Antigens, mode=1, controls=get("controls", globalenv())){
+  lapply(1:length(dflist), function(setid){
+    set = dflist[[setid]]
+    antigens = antigens [[setid]]
+    for (i in 3:length(set)){
+      analytenum = as.numeric(strsplit(colnames(set[i]), split='.', fixed = TRUE)[[1]][2])
+      if (mode == 1){
+        colnames(set)[i] = antigens[antigens$analyte == analytenum,]$Antigen.name
+      } else if (mode == 2){
+        colnames(set)[i] = antigens[antigens$analyte == analytenum,]$Gene.name
+      }}
+    set = set[,!(names(set) %in% controls)]
+    return(set)
+  })
+} 
+
+
+### Automatic Processing
+stage_1 <- function(dflist, name){
+  T01 = handle_background(dflist)
+  T02 = lapply(T01,emptyadjust)
+  T03 = lapply(T02,compress_duplicates, layout="data/layout.xlsx")
+  T04 = set_colname_adapter(T03)
+  assign(name, T04, envir = .GlobalEnv)
+}

xmap_biomarkers/scripts/backend/preprocessing_stage_2.R

@@ -0,0 +1,40 @@
+# The function adds a prefix "Set_" to the index of each data frame to create the name for that data frame. 
+# The function returns a list of these generated names.
+simple_names <- function (dflist){
+  namelist = list()
+  for (i in 1:length(dflist)){
+    namelist[i] = paste0("Set_", i)
+  }
+  return(namelist)
+}
+
+# The function mergedown takes a list of data frames dflist and merges them into a single data frame by 
+# taking the row mean of columns with the same name in each data frame. It returns the merged data frame. 
+# If there are columns in a data frame that are not present in any other data frame, they are included in the merged data frame as is.
+mergedown <- function (dflist){
+  outputdf = dflist[[1]]
+  for (setid in 2:length(dflist)){
+    originlist = colnames(dflist[[1]])[-1:-2]
+    mergelist = colnames(dflist[[setid]])[-1:-2]
+    uniquelist = mergelist[!(mergelist %in% originlist)]
+    mergelist = mergelist[mergelist %in% originlist]}
+  for (name in mergelist) {
+    #print(name)
+    outputdf[,name] = rowMeans(data.frame(dflist[[1]][,name], dflist[[setid]][,name]), na.rm = TRUE)
+  }
+  if (length(uniquelist) != 0){
+    for (name in uniquelist) {
+      #print(name)
+      outputdf[name] = dflist[[setid]][,name]
+    }
+  }
+  return(outputdf)
+}
+
+#Automatic Processing
+stage_2 <- function(dflist, name){
+  set3 <- mergedown(dflist)
+  dflist <- c(dflist, list(set3))
+  names(dflist) <- simple_names(dflist)
+  assign(name, dflist, envir = .GlobalEnv)
+}

xmap_biomarkers/scripts/backend/pretty_pictures.R

@@ -0,0 +1,158 @@
+pp_dflist_wrapper <- function(dflist, pp_function){
+  dflist_name = deparse(substitute(dflist))
+  plot_name = paste(strsplit(deparse(substitute(pp_function)), "_")[[1]][-1], collapse = "")
+  dir_path = file.path(getwd(),"plots",dflist_name,plot_name)
+  dir.create(dir_path, recursive=TRUE, showWarnings = FALSE)
+  for (setid in seq_along(dflist)){
+    set_name = paste("Set",setid,sep="_")
+    name = (file.path(dir_path,set_name))
+    pp_function(dflist[[setid]], name)
+  }
+}
+
+pca_plot <- function(data, show_ellipse = TRUE) {
+  library(ggplot2)
+  library(ggfortify)
+  
+  data$group <- factor(data$group, levels = c(0, 1), labels = c("Healthy", "Diseased"))
+  pca_data <- prcomp(data[-1:-2], center = TRUE, scale = TRUE)
+  
+  # Extract PCA scores
+  pca_scores <- as.data.frame(pca_data$x)
+  pca_scores$group <- data$group
+  # Calculate percentage of variance explained by each PC
+  var_exp <- round(pca_data$sdev^2 / sum(pca_data$sdev^2) * 100, 2)
+  
+  # Update axis labels with percentage of variance explained
+  x_label <- paste0("PC1 (", var_exp[1], "%)")
+  y_label <- paste0("PC2 (", var_exp[2], "%)")
+  
+  plot <- ggplot(pca_scores, aes(x = PC1, y = PC2, color = group)) +
+    geom_point() +
+    theme_classic() +
+    labs(x = x_label, y = y_label, title = "PCA Plot") + 
+    scale_color_manual(values = c("Healthy" = "blue", "Diseased"="red"))
+  
+  if (show_ellipse) {
+    plot <- plot + stat_ellipse(aes(fill = group), geom = "polygon", level = 0.95, alpha = 0.2) +
+      labs(title = "") + 
+      scale_fill_manual(values = c("Healthy" = "palegreen", "Diseased"="palegoldenrod"))
+  }
+  
+  return(plot)
+}
+
+
+
+pp_box_plot_multi <- function(data, name){
+  # Set up plot area and device
+  dim = display_division(ncol(data)-2)
+  if ((ncol(data)-2)<=(dim[1]*dim[2])){
+    filename = name
+    col_range = 3:ncol(data)
+    png(filename = paste0(filename,".png"), width = 1200+300*dim[1], height = 900+100*dim[2], res=250)
+    par(mfrow = c(dim[2],dim[1]))
+    for (i in col_range) {
+      plot = boxplot(data[,i] ~ data$group, main = colnames(data)[i], xlab = "Group", ylab="")
+    }
+    dev.off()
+  }
+  else{
+    modifier = dim[1]*dim[2]
+    for (d in seq(dim[3])){
+      filename = paste(name,d,sep="_")
+      png(filename = paste0(filename,".png"), width = 1200+300*dim[1], height = 900+300*dim[2], res=250)
+      col_range = 3:(modifier+2)+(modifier*(d-1))
+      col_range = col_range[col_range<=ncol(data)]
+      par(mfrow = c(dim[2],dim[1]))
+      print(col_range)
+      for (i in col_range) {
+        plot = boxplot(data[,i] ~ data$group, main = colnames(data)[i], xlab = "Group", ylab="")
+      }
+      dev.off()
+    }
+  }
+}
+
+library(ggplot2)
+#library(ggbreak)
+
+pp_gg_box_plot_multi <- function(data, name){
+  
+  data <- data[,-1]
+  # Transform data to long format
+  data_long <- tidyr::pivot_longer(data, -group, names_to="Variable", values_to="Value")
+  
+  # Calculate upper limit for normal scale before break
+  upper_limit <- quantile(data_long$Value, 0.95)
+  
+  # Create the plot
+  p <- ggplot(data_long, aes(x=group, y=Value)) + 
+    geom_boxplot(aes(group=group, fill=factor(group)), outlier.shape = NA, color="black") +  # Set the boxes to neutral color and outline them in black
+    geom_point(aes(color=factor(group)), position = position_jitter(width = 0.3), alpha=0.7) +  # Keep the points colored
+    facet_wrap(~Variable, scales="free_y") + 
+    coord_trans(y="log10") +  # Implement log transform
+    theme_bw() + 
+    scale_fill_manual(values = c("white", "white"), guide=FALSE) +  # Use white color for box fill and disable its legend
+    scale_color_manual(values = c("#E57373", "#4DB6AC"), 
+                       labels = c("Healthy Controls", "XFG Patients"),
+                       name = "") +  # Return to the preferred colors
+    theme(strip.background = element_blank(),
+          strip.text = element_text(size=12, face="bold"),
+          axis.title.x = element_blank(),
+          axis.title.y = element_blank(),
+          axis.ticks.x = element_blank(),
+          axis.text.x = element_blank(),
+          legend.position = c(0.9, 0.1),
+          legend.justification = c(1, 0),
+          legend.background = element_blank(),
+          legend.key = element_blank())
+  
+  # Save the plot to a file
+  ggsave(filename = paste0(name, ".png"), plot = p, width = 10, height = 6)
+}
+
+
+
+pp_heatmap <- function(df, subset = "full", show_row_names = TRUE, show_col_names = TRUE){
+  # Apply the appropriate subset function based on the 'subset' parameter
+  if(subset == "limma"){
+    df = limma_subset(df)
+  } else if(subset == "compstat"){
+    df = compstat_subset(df)
+  }
+  # Remove the first two columns to create the 'subset' dataframe
+  subset = df[-1:-2]
+  # Standardize the columns of the 'subset' dataframe
+  subset = apply(subset, 2, function(x) (x - mean(x)) / sd(x))
+  # Create row names for the 'subset' dataframe based on the 'group' and 'Internal.LIMS.ID' columns
+  rownames(subset) = paste0(ifelse(df$group == 0, "Healthy", "Diseased"), "-", substr(df$Internal.LIMS.ID, start = nchar(df$Internal.LIMS.ID)-3, stop = nchar(df$Internal.LIMS.ID)))
+  # Order the columns of the 'subset' dataframe by decreasing column means
+  subset = subset[, order(colMeans(subset), decreasing = TRUE)]
+  # Transpose the 'subset' dataframe
+  subset = t(subset)
+  
+  pheatmap(subset, 
+           scale = "none", 
+           cluster_rows = TRUE, 
+           cluster_cols = TRUE, 
+           show_rownames = show_row_names, 
+           show_colnames = show_col_names, 
+           treeheight_row = 0,
+           clustering_distance_cols = "euclidean", 
+           clustering_distance_rows = "euclidean",
+           clustering_method = "complete")
+}
+
+pp_multipca <- function(df, name){
+  datalist <- list(base_data=df,limma_data=limma_subset(df),compstat_data = compstat_subset(df))
+  plots <- lapply(names(datalist), function(name) {
+    plot <- pca_plot(datalist[[name]])
+    plot + ggtitle(name)
+  })
+  combined_plot <- do.call(grid.arrange, c(plots, ncol = 3))
+  ggsave(paste(name, "multipca.png", sep="_"), combined_plot, width = 12, height = 4, dpi = 300)
+  }
+
+
+  

xmap_biomarkers/scripts/backend/transformation.R

@@ -0,0 +1,171 @@
+##Basic-ish Maths
+geometric_mean <- function(numbers){
+  gm = prod(numbers)^(1/length(numbers))
+  return(gm)
+}
+##Diff Analysis
+limma_funct <- function(data) {
+  t_set=t(data[-1:-2])
+  design <- model.matrix(~0 + group, data=data)
+  colnames(design) <- c("case", "control")
+  contrasts = makeContrasts(Diff= control - case, levels=design)
+  fit<-lmFit(t_set, design, method="robust", maxit=1000)
+  contrast_fit <- contrasts.fit(fit,contrasts)
+  ebay_fit <- eBayes(contrast_fit)
+  DE_results <- topTable(ebay_fit, n=ncol(data), adjust.method = "fdr", confint = TRUE)
+  print(summary(decideTests(ebay_fit)))
+  return(DE_results)
+}
+sig_test <- function(data){
+  sigtestlist = data.frame()
+  for (name in colnames(data[-1:-2])){
+    #print(name)
+    # print(head(data))
+    setC = data[data$group==0,][[name]]
+    setE = data[data$group==1,][[name]]
+    ShapE = shapiro.test(setE)
+    ShapC = shapiro.test(setC)
+    if (ShapE$p.value >0.05 & ShapC$p.value >0.05){
+      testrow = cbind(tidy(t.test(setC, setE))[c("statistic", "p.value", "method")],ShapC$p.value, ShapE$p.value)
+    } else {
+      testrow = cbind(tidy(wilcox.test(setC, setE))[c("statistic", "p.value", "method")],ShapC$p.value, ShapE$p.value)
+    }
+    testrow$analyte = name
+    sigtestlist = rbind(sigtestlist, testrow)
+  }
+  return(sigtestlist)
+}
+comparative_statistics <- function(df) {
+  sigtest <- sig_test(df)
+  sigtest$bh_p.value <- p.adjust(sigtest$p.value, method = "BH")
+  qobj <- qvalue(p = sigtest$p.value)
+  sigtest$q.value <- qobj$qvalues
+  setC <- df[df$group == 0, ]
+  setE <- df[df$group == 1, ]
+  FC <- apply(setE[, -c(1:2)], 2, function(x) mean(x, na.rm = TRUE)) / 
+    apply(setC[, -c(1:2)], 2, function(x) mean(x, na.rm = TRUE))
+  log2FC <- log2(FC)
+  P.Value <- sigtest$p.value
+  Method <- sigtest$method
+  BH_P.Value <- sigtest$bh_p.value
+  Q.Value <- sigtest$q.value
+  ShapE <- sigtest$'ShapE$p.value'
+  ShapC <- sigtest$'ShapC$p.value'
+  comp <- data.frame(FC, log2FC, P.Value, Method, BH_P.Value, Q.Value, ShapE, ShapC)
+  row.names(comp) <- colnames(df)[-c(1:2)]
+  return(comp)
+}
+## Wrappers
+differential_reports <- function(dflist){
+  output <- list()
+  for (i in seq_along(dflist)) {
+    #print(head(dflist[[i]]))
+    df_name = paste0("Set_", i)
+    limma_output <- limma_funct(dflist[[i]])
+    limma_name <- paste0(df_name, "_limma")
+    output[[limma_name]] <- limma_output
+    comp_output <- comparative_statistics(dflist[[i]])
+    comp_name <- paste0(df_name, "_comparative_stats")
+    output[[comp_name]] <- comp_output
+  }
+  return(output)
+}
+limma_subset <- function(df, mode = "default", n=15, P=0.05){
+  lim <- limma_funct(df)
+  
+  if (mode == "default"){
+    lim_sigs <- row.names(lim[lim$adj.P.Val<P,])
+    lim_data <- cbind(df[1:2],df[lim_sigs])
+  }
+  else if (mode == "raw"){
+    lim_sigs <- row.names(lim[lim$P.Val<P,])
+    lim_data <- cbind(df[1:2],df[lim_sigs])
+  }
+  else if (mode == "top"){
+    lim_sigs <- lim %>% arrange(adj.P.Val) %>% head(n) %>% row.names
+    lim_data <- cbind(df[1:2], df[lim_sigs])
+  }
+  return(lim_data)
+}
+compstat_subset <- function(df){
+  comp <- comparative_statistics(df)
+  comp_sigs <- row.names(comp[comp$Q.Value<0.05,])
+  comp_data <- cbind(df[1:2],df[comp_sigs])
+  return(comp_data)
+}
+clustering_dflist_wrapper<- function(dflist, clust_function){
+  dflist_name = deparse(substitute(dflist))
+  plot_name = deparse(substitute(clust_function))
+  dir_path = file.path(getwd(),"plots",dflist_name,plot_name)
+  dir.create(dir_path, recursive=TRUE, showWarnings = FALSE)
+  for (setid in seq_along(dflist)){
+    set_name = paste("Set",setid,sep="_")
+    name = (file.path(dir_path,set_name))
+    clust_function(dflist[[setid]], name)
+  }
+}
+##Transformers
+#BoxCox
+transformer_boxcox <- function(df, weighted = TRUE){
+  t_df = df
+  for (colname in colnames(df)[-1:-2]){
+    #print(colname)
+    lambda = determine_lambda(colname, df)
+    if (weighted == TRUE){
+      t_df[[colname]]=bc_weighted_transform(df[[colname]],lambda)        
+    }
+    else{
+      t_df[[colname]]=bc_transform(df[[colname]],lambda)
+    }
+    
+  }
+  rm(dataf,column,envir=globalenv())
+  return(t_df)
+}
+##Boxcox Modules
+determine_lambda <- function(colname, df) {
+  dataf <<- as.data.frame(df)
+  column <<- df[,colname]
+  # column <<- column
+  model <- lm(column ~ group, data = dataf)
+  bc <- boxcox(model, lambda = seq(-5, 5))
+  lambda <- bc$x[which(bc$y==max(bc$y))]
+  return(lambda)
+}
+bc_transform <- function(y, lambda=0) {
+  if (lambda == 0L) { log(y) }
+  else { (y^lambda - 1) / lambda }
+}
+bc_weighted_transform <- function(y, lambda=0) {
+  geom = geometric_mean(y)
+  if (lambda == 0L) { log(y) }
+  else { (y^lambda - 1) / lambda*geom^(lambda-1)}
+}
+#RSN
+transformer_rsn <- function(df){
+  subset = as.matrix(df[-1:-2])
+  sink(nullfile <- tempfile())
+  rsn_transform = lumiN(subset, method= "rsn")
+  sink(NULL)
+  na_cols <- which(colSums(is.na(rsn_transform)) > 0)
+  if (length(na_cols) > 0) {
+    cat("Columns with NAs:", colnames(rsn_transform)[na_cols], "\n")
+    rsn_transform <- rsn_transform[, -na_cols]
+  }
+  return(cbind(df[1:2],rsn_transform))
+}
+
+
+
+##Automatic Processing
+stage_3 <- function (dflist, name){
+  trans_list <- lapply(dflist, transformer_boxcox)
+  trans_list_rsn <- lapply(trans_list, transformer_rsn)
+  bc_report_name <- paste0("reports_",name,"_bc")
+  bcrsn_report_name <- paste0("reports_",name,"_bcrsn")
+  assign(bc_report_name, differential_reports(trans_list), envir = .GlobalEnv)
+  assign(bcrsn_report_name, differential_reports(trans_list_rsn), envir = .GlobalEnv)
+  assign(name, trans_list_rsn, envir = .GlobalEnv)
+}
+
+