main.R

library(tercen)
library(dplyr)
library(tidyr)
library(tibble)

# Set verbose flag for detailed logging
verbose <- TRUE

# Function to log messages if verbose is TRUE
log_message <- function(message) {
  if (verbose) {
    cat(paste0("[Hypergate Operator] ", message, "\n"))
  }
}

log_message("Starting Hypergate Operator")

# Load hypergate package
library(hypergate)

# Get Tercen context
# http://127.0.0.1:5400/kumo_test/w/0b4ec8f92d47f8a388859024bb0186e7/ds/4a4a5702-7e70-4266-bba8-a00309029928
# options("tercen.workflowId"= "0b4ec8f92d47f8a388859024bb0186e7")
# options("tercen.stepId"= "4a4a5702-7e70-4266-bba8-a00309029928")
# ctx <- tercenCtx(username = "test", password = "test", serviceUri = "http://127.0.0.1:5400")
ctx <- tercenCtx()


# Get operator properties
#CD3⁺CD4⁺CD45RA⁺CCR7⁺CD62L⁺CD27⁺CD8⁺-
# CD3+ CD4+ CD45RA+ CCR7+ CD27+ CD8-
high_expression_threshold <- ctx$op.value('high_expression', as.numeric, 0.5)
low_expression_threshold <- ctx$op.value('low_expression', as.numeric, 0.5)
expressed_markers_str <- ctx$op.value('expressed_marker', as.character, '')
not_expressed_markers_str <- ctx$op.value('not_expressed_marker', as.character, '')

log_message(paste0("High expression threshold: ", high_expression_threshold))
log_message(paste0("Low expression threshold: ", low_expression_threshold))
log_message(paste0("Expressed markers: ", expressed_markers_str))
log_message(paste0("Not expressed markers: ", not_expressed_markers_str))

# Parse comma-separated marker lists
expressed_markers <- if (expressed_markers_str != "") {
  strsplit(expressed_markers_str, ",")[[1]]
} else {
  character(0)
}

not_expressed_markers <- if (not_expressed_markers_str != "") {
  strsplit(not_expressed_markers_str, ",")[[1]]
} else {
  character(0)
}

# Validate that at least one marker is specified
if (length(expressed_markers) == 0 && length(not_expressed_markers) == 0) {
  stop("At least one marker must be specified in either expressed_marker or not_expressed_marker")
}

# Get data from Tercen context
log_message("Extracting data from Tercen context")

# Get main data with .y values
main_data <- ctx %>% 
  select(.y, .ci, .ri) 

# Get column projections (event IDs)
col_data <- ctx %>% 
  cselect() %>%
  mutate(.ci = seq_len(n()) - 1)  # Add .ci index starting from 0

# Get row projections (variable/channel information)
row_data <- ctx %>% 
  rselect() %>%
  mutate(.ri = seq_len(n()) - 1)  # Add .ri index starting from 0

log_message(paste0("Main data rows: ", nrow(main_data)))
log_message(paste0("Column data rows: ", nrow(col_data)))
log_message(paste0("Row data rows: ", nrow(row_data)))

# Transform data from long to wide format
log_message("Transforming data to wide format")

# Join main data with row data to get marker information
data_with_markers <- main_data %>%
  left_join(row_data, by = ".ri")

# Reshape to wide format where rows are events and columns are markers
# We need to handle potential memory issues with large datasets
# Process in chunks if needed

# Function to process data in chunks
process_in_chunks <- function(data, chunk_size = 100000) {
  log_message(paste0("Processing data in chunks of size ", chunk_size))
  
  # Get unique event IDs (.ci values)
  event_ids <- unique(data$.ci)
  n_events <- length(event_ids)
  
  # Calculate number of chunks
  n_chunks <- ceiling(n_events / chunk_size)
  
  # Initialize result list
  result_list <- list()
  
  for (i in 1:n_chunks) {
    log_message(paste0("Processing chunk ", i, " of ", n_chunks))
    
    # Get event IDs for this chunk
    start_idx <- (i - 1) * chunk_size + 1
    end_idx <- min(i * chunk_size, n_events)
    chunk_event_ids <- event_ids[start_idx:end_idx]
    
    # Filter data for this chunk
    chunk_data <- data %>%
      filter(.ci %in% chunk_event_ids)
    
    # Process chunk
    chunk_result <- chunk_data %>%
      pivot_wider(
        id_cols = .ci,
        names_from = logicle..channel_description,
        values_from = .y,
        values_fn = list(.y = mean)  # Use mean if multiple values per cell
      )
    
    # Add to result list
    result_list[[i]] <- chunk_result
  }
  
  # Combine results
  do.call(rbind, result_list)
}

# Check if we need to process in chunks (based on data size)
if (nrow(data_with_markers) > 1000000) {
  wide_data <- process_in_chunks(data_with_markers)
} else {
  wide_data <- data_with_markers %>%
    pivot_wider(
      id_cols = .ci,
      names_from = logicle..channel_description,
      values_from = .y,
      values_fn = list(.y = mean)  # Use mean if multiple values per cell
    )

}

log_message(paste0("Wide data dimensions: ", nrow(wide_data), " x ", ncol(wide_data)))

# Identify population of interest based on marker expression
log_message("Identifying population of interest")

# Initialize a logical vector for population membership
population_membership <- rep(TRUE, nrow(wide_data))

# Apply high expression thresholds
for (marker in expressed_markers) {
  if (marker %in% colnames(wide_data)) {
    marker_values <- wide_data[[marker]]
    threshold <- quantile(marker_values, high_expression_threshold, na.rm = TRUE)
    population_membership <- population_membership & (marker_values > threshold)
    log_message(paste0("Applied high expression threshold for ", marker, ": ", threshold))
  } else {
    log_message(paste0("Warning: Marker ", marker, " not found in data"))
  }
}

# Apply low expression thresholds
for (marker in not_expressed_markers) {
  if (marker %in% colnames(wide_data)) {
    marker_values <- wide_data[[marker]]
    threshold <- quantile(marker_values, low_expression_threshold, na.rm = TRUE)
    population_membership <- population_membership & (marker_values < threshold)
    log_message(paste0("Applied low expression threshold for ", marker, ": ", threshold))
  } else {
    log_message(paste0("Warning: Marker ", marker, " not found in data"))
  }
}

# Count events in population of interest
n_population <- sum(population_membership, na.rm = TRUE)
log_message(paste0("Population of interest contains ", n_population, " events"))

# Prepare data for hypergate
log_message("Preparing data for hypergate")

# Remove .ci column for hypergate input
xp_data <- wide_data %>%
  select(-(.ci)) %>%
  as.matrix()

# Create gate vector (1 for population of interest, 0 for others)
gate_vector <- ifelse(population_membership, 1, 0)

# Apply hypergate algorithm
log_message("Applying hypergate algorithm")

# Handle potential errors in hypergate
tryCatch({
  # Run hypergate
  hg_result <- hypergate(
    xp = xp_data,
    gate_vector = gate_vector,
    level = 1,  # Our population of interest is labeled as 1
    verbose = verbose,
    beta = 100,  # Equal weight to precision and recall
    delta_add = 0.0001 #
  )
  
  # Apply the gate to get membership
  gate_membership <- subset_matrix_hg(hg_result, xp_data)
  
  log_message(paste0("Hypergate identified ", sum(gate_membership), " events in the gate"))
  
  # Prepare output data
  output_data <- tibble(
    .ci = wide_data$.ci,
    hypergate_member = as.integer(gate_membership)
  )
  
}, error = function(e) {
  # If hypergate fails, use the initial population membership
  log_message(paste0("Hypergate error: ", e$message))
  log_message("Using initial population membership as fallback")
  
  output_data <- tibble(
    .ci = wide_data$.ci,
    hypergate_member = as.integer(population_membership)
  )
})

# Verify output has same number of rows as column projection
if (nrow(output_data) != nrow(col_data)) {
  log_message("Warning: Output rows don't match column projection rows")
}

# Save results to Tercen
log_message("Saving results to Tercen")
output_data %>%
  ctx$addNamespace() %>%
  ctx$save()

log_message("Hypergate operator completed successfully")