#!/usr/bin/env Rscript
# reproduce_hcutoff_from_expression.R
# Template for re-running h-cutoff on a real expression matrix.
# Input:
#   1. expression matrix CSV, genes as rows, samples as columns
#   2. gene list CSV with one column named gene
# Output:
#   edge weights, h-subnet, h-backbone and two-order Hub genes.
args <- commandArgs(trailingOnly=TRUE)
if (length(args) < 2) {
  stop("Usage: Rscript reproduce_hcutoff_from_expression.R expression.csv gene_list.csv")
}
expr <- read.csv(args[1], row.names=1, check.names=FALSE)
genes <- read.csv(args[2], stringsAsFactors=FALSE)$gene
genes <- intersect(genes, rownames(expr))
x <- t(expr[genes,,drop=FALSE])
cor_mat <- abs(cor(x, method="pearson", use="pairwise.complete.obs"))
diag(cor_mat) <- NA

continuous_h <- function(v, scale=100) {
  z <- sort(v[is.finite(v)] * scale, decreasing=TRUE)
  h <- 0
  for (i in seq_along(z)) {
    if (z[i] >= i) h <- i else break
  }
  return(list(h=h, cutoff=h/scale))
}

edge_list <- data.frame()
for (i in seq_len(nrow(cor_mat)-1)) {
  for (j in (i+1):ncol(cor_mat)) {
    edge_list <- rbind(edge_list, data.frame(source=rownames(cor_mat)[i], target=colnames(cor_mat)[j], weight=cor_mat[i,j]))
  }
}
h1 <- continuous_h(edge_list$weight, scale=100)
edge_list$one_order_status <- ifelse(edge_list$weight >= h1$cutoff, "h-subnet", "candidate-weak-bridge")
write.csv(edge_list, "edge_weights_from_expression.csv", row.names=FALSE)
message("one-order h = ", h1$h, "; edge cutoff = ", h1$cutoff)
# For weak-bridge and two-order hub extraction, use the Python/networkx script in this package,
# or calculate edge betweenness using igraph and repeat the node-level continuous_h procedure.
