Extract Z-scores from model.interpretation.plot

[flowerpotexpt]

Hey!

Not really an issue, but I wondered if I could extract the model Z-score values for each feature and sample as shown in the plot's heat map? That would be a lot of help!

Thanks in advance!

[scilexenko]

Hi,
It may help

extract_siamcat_plot_data <- function(siamcat, consens.thres = 0.5, heatmap.type = "zscore", 
                                      max.show = 50, log.n0 = 1e-06) {
  # Get basic components
  feature.weights <- SIAMCAT::feature_weights(siamcat)
  weight.matrix <- SIAMCAT::weight_matrix(siamcat)
  model.type <- SIAMCAT::model_type(siamcat)
  label <- SIAMCAT::label(siamcat)
  
  # Since we can't directly call model.interpretation.select.features,
  # let's implement a simplified version of feature selection
  if (label$type == "BINARY") {
    # For binary classification, select features based on their median weights
    median_weights <- apply(weight.matrix, 1, median)
    sorted_weights <- sort(abs(median_weights), decreasing = TRUE)
    # Select top features up to max.show
    sel.idx <- names(sorted_weights)[1:min(max.show, length(sorted_weights))]
  } else {
    # For regression, just select top features by absolute weight
    median_weights <- apply(weight.matrix, 1, median)
    sorted_weights <- sort(abs(median_weights), decreasing = TRUE)
    sel.idx <- names(sorted_weights)[1:min(max.show, length(sorted_weights))]
  }
  
  # Get feature data based on feature type
  feat <- NULL
  if (heatmap.type == "zscore") {
    if (SIAMCAT::feature_type(siamcat) == "original") {
      feat <- SIAMCAT::get.orig_feat.matrix(siamcat)
    } else if (SIAMCAT::feature_type(siamcat) == "filtered") {
      feat <- SIAMCAT::get.filt_feat.matrix(siamcat)
    } else if (SIAMCAT::feature_type(siamcat) == "normalized") {
      feat <- SIAMCAT::get.norm_feat.matrix(siamcat)
    }
  } else {
    feat <- SIAMCAT::get.orig_feat.matrix(siamcat)
  }
  
  # Get prediction data
  pred_matrix <- SIAMCAT::pred_matrix(siamcat)
  mean.agg.pred <- rowMeans(pred_matrix)
  
  # Sample sorting
  if (label$type == "BINARY") {
    srt.idx <- sort(label$label + mean.agg.pred, index.return = TRUE)$ix
  } else if (label$type == "CONTINUOUS") {
    srt.idx <- sort(label$label, index.return = TRUE)$ix
  }
  
  # Since we can't call the specific heatmap preparation functions,
  # let's create simplified versions
  img.data <- NULL
  if (heatmap.type == "zscore") {
    # Z-score normalization
    selected_features <- feat[sel.idx, srt.idx, drop = FALSE]
    img.data <- t(scale(t(selected_features)))
  } else if (heatmap.type == "fc") {
    # Fold change calculation (simplified)
    if (label$type != "BINARY") {
      stop("Fold-change heatmap only works with binary labels")
    }
    
    # Get feature data for selected features
    selected_features <- feat[sel.idx, srt.idx, drop = FALSE]
    
    # Calculate fold change (very simplified approach)
    cases <- which(label$label[srt.idx] == 1)
    controls <- which(label$label[srt.idx] == 0)
    
    # Add pseudo-count and log-transform
    if (is.null(log.n0)) {
      log.n0 <- quantile(feat[feat != 0], 0.05)
    }
    
    # Calculate fold change (this is a simplified version)
    case_means <- rowMeans(selected_features[, cases, drop = FALSE] + log.n0)
    control_means <- rowMeans(selected_features[, controls, drop = FALSE] + log.n0)
    img.data <- log2(case_means / control_means)
  } else if (heatmap.type == "log") {
    # Log transformation
    if (is.null(log.n0)) {
      log.n0 <- quantile(feat[feat != 0], 0.05)
    }
    selected_features <- feat[sel.idx, srt.idx, drop = FALSE]
    img.data <- log10(selected_features + log.n0)
  }
  
  # Create relative weights
  rel.weights <- t(t(weight.matrix)/colSums(abs(weight.matrix), na.rm = TRUE))
  
  # Return data as a list
  return(list(
    feature_data = feat,
    selected_features = sel.idx,
    sorted_samples = srt.idx,
    feature_weights = feature.weights,
    weight_matrix = weight.matrix,
    relative_weights = rel.weights[sel.idx, , drop = FALSE],
    predictions = mean.agg.pred,
    labels = label,
    metadata = SIAMCAT::meta(siamcat),
    heatmap_data = img.data
  ))
}

Ivan