In vivekkohar/test: Systems biology tool to simulate gene regulatory circuits

knitr::opts_chunk$set(echo = TRUE, eval = FALSE)

Evaluate the networks

Once the network simulations are done, one can compare the simulations to the expression data.

Comparing data

fileNames <- list.files("simAll/")
nClusters <- 2
networkMetrics <- data.frame(Nodes=integer(),Interactions=integer(),PosInt = integer(),ClusterA = integer(), ClusterB = integer(),Connected = logical(), ModelCount = integer(), transitivity = numeric(), MeanDistance = numeric(), Accuracy = numeric(), Filename = character(), MiPos = numeric(), MiNeg = numeric())

for(fileCounter in seq_along(fileNames)){

  fileNameTmp <- fileNames[fileCounter]
  print(fileNameTmp)

posMi <- as.numeric(gsub("^(?:[^_]+_){1}([^_]+).*", "\\1", fileNameTmp))
negMi <- gsub("^(?:[^_]+_){2}([^_]+).*", "\\1", fileNameTmp)
negMi <- as.numeric(substr(negMi,0,nchar(negMi)-4))

racipe <- readRDS(paste0("simAll/",fileNameTmp))
racipe <- sracipeNormalize(racipe)
expData <- assay(racipe,1)

dataRef <- # your exp data
dataSim <- expData
hS <- sracipeHeatmapSimilarity(dataReference = tmp1,dataSimulation = tmp2,nClusters = nClusters)
ClusterA <- hS$simulated.cluster.freq[2]
ClusterB <- hS$simulated.cluster.freq[3]
circuit <-  as.data.frame(sracipeCircuit(racipe))
g <- graph_from_data_frame(circuit, directed = TRUE, vertices = NULL)
expDataTmp <- as.matrix(expData)
storage.mode(expDataTmp) <- "logical"
networkMetricsTmp <- data.frame(Nodes= length(union(circuit$Source,circuit$Target)),Interactions=length(circuit$Source),PosInt = length(which(circuit$Type ==1)),ClusterA = ClusterA, ClusterB = ClusterB,Connected = igraph::is_connected(g), ModelCount = sracipeConfig(racipe)$simParams["numModels"], Transitivity = igraph::transitivity(g),MeanDistance = igraph::mean_distance(g, directed = TRUE, unconnected = FALSE), Filename=fileNameTmp, MiPos = posMi, MiNeg = negMi)

networkMetricsTmp$Accuracy <- (networkMetricsTmp$ClusterA + networkMetricsTmp$ClusterB)
networkMetricsTmp$Accuracy2 <- hS$KL
networkMetricsTmp$ClusterA2 <- hS$cluster.similarity[2]
networkMetricsTmp$ClusterB2 <- hS$cluster.similarity[3]
networkMetrics <- rbind(networkMetrics, networkMetricsTmp)
}
saveRDS(networkMetrics, file = "networkMetrics.rds")

Comparing binarized data

EMClassCondition <- readRDS("EMClassCondition.rds")
hammingTable <- c(0,7,12,17,22,28,35,42,49,56,64,71,79)
fileNames <- list.files("simAll/")
networkMetrics <- data.frame(Nodes=integer(),Interactions=integer(),PosInt = integer(),EModels = integer(), MModels = integer(),Connected = logical(),MGenes = integer(), Egenes = integer(), ModelCount = integer(), CellLine = character(), Signal = character(), Dataset = character(), transitivity = numeric(), MeanDistance = numeric(), Accuracy = numeric(), Filename = character(), MiPos = numeric(), MiNeg = numeric())

for(fileCounter in seq_along(fileNames)){

  fileNameTmp <- fileNames[fileCounter]
  print(fileNameTmp)

cellLine <- str_match(fileNameTmp, "(.*?)_")[,2]
signal <- str_match(fileNameTmp, "_(.*?)_")[,2]
posMi <- as.numeric(gsub("^(?:[^_]+_){2}([^_]+).*", "\\1", fileNameTmp))
negMi <- gsub("^(?:[^_]+_){3}([^_]+).*", "\\1", fileNameTmp)
negMi <- as.numeric(substr(negMi,0,nchar(negMi)-4))

condition <- paste0(cellLine, "_",signal)

racipe <- readRDS(paste0("simAll/",fileNameTmp))
# sracipePlotCircuit(racipe,plotToFile = FALSE)
racipe <- sracipeNormalize(racipe)
expData <- assay(racipe,1)

EMGenes <- EMClassCondition[,condition]
MState <- EMGenes[!is.na(EMGenes)]
circuit <-  as.data.frame(sracipeCircuit(racipe))
g <- graph_from_data_frame(circuit, directed = TRUE, vertices = NULL)
# igraph::is_connected(g)



commonGenes <- rownames(expData)
commonGenes <- commonGenes[which(commonGenes %in% names(MState))]
expData <- expData[commonGenes,]
expData <- Binarize::binarizeMatrix(expData, method = "kMeans")[,1:dim(expData)[2]]

MState <- MState[commonGenes]
EState <- !MState
expDataTmp <- as.matrix(expData)
storage.mode(expDataTmp) <- "logical"
if(length(MState)>85) message("Error")
hammingCutoff <- findInterval(length(MState),hammingTable)
hamDistE <- apply((expDataTmp), 2, function(x) sum(x != EState))
hamDistM <- apply((expDataTmp), 2, function(x) sum(x != MState))
networkMetricsTmp <- data.frame(Nodes= length(union(circuit$Source,circuit$Target)),Interactions=length(circuit$Source),PosInt = length(which(circuit$Type ==1)),EModels = length(which(hamDistE < hammingCutoff)), MModels = length(which(hamDistM < hammingCutoff)),Connected = igraph::is_connected(g),MGenes = sum(MState), Egenes = sum(EState), ModelCount = sracipeConfig(racipe)$simParams["numModels"], CellLine = cellLine, Signal = signal, Dataset = condition, Transitivity = igraph::transitivity(g),MeanDistance = igraph::mean_distance(g, directed = TRUE, unconnected = FALSE), Filename=fileNameTmp, MiPos = posMi, MiNeg = negMi)
networkMetricsTmp$Accuracy <- sum(networkMetricsTmp$EModels + networkMetricsTmp$MModels)/networkMetricsTmp$ModelCount
rownames(networkMetricsTmp) <- condition
networkMetrics <- rbind(networkMetrics, networkMetricsTmp)
}
saveRDS(networkMetrics, file = "networkMetrics.rds")