library(dplyr, lib.loc = "/pstore/apps/bioinfo/R/4.0.1-foss")
library(patchwork, lib.loc = "/pstore/home/valdeola/R/x86_64-pc-linux-gnu-library/4.0.1-foss")
library(Seurat, lib.loc = "/pstore/apps/bioinfo/R/4.0.1-foss")
library(Rcpp, lib.loc = "/pstore/apps/bioinfo/R/4.0.1-foss")
library(mistyR, lib.loc = "/pstore/apps/bioinfo/R/4.0.1-foss")
library(purrr, lib.loc = "/pstore/apps/bioinfo/R/4.0.1-foss")
library(tibble, lib.loc = "/pstore/apps/bioinfo/R/4.0.1-foss")
library(readr, lib.loc = "/pstore/apps/bioinfo/R/4.0.1-foss")
library(stringr, lib.loc = "/pstore/apps/bioinfo/R/4.0.1-foss")
library(ggplot2, lib.loc = "/pstore/apps/bioinfo/R/4.0.1-foss")
library(cowplot, lib.loc = "/pstore/apps/bioinfo/R/4.0.1-foss")
library(RColorBrewer, lib.loc = "/pstore/apps/bioinfo/R/4.0.1-foss")
library(liana,  lib.loc = "/pstore/apps/bioinfo/R/4.0.1-foss")
library(kableExtra, lib.loc = "/pstore/apps/bioinfo/R/4.0.1-foss")
library(igraph, lib.loc = "/pstore/apps/bioinfo/R/4.0.1-foss")
library(OmnipathR, lib.loc = "/pstore/apps/bioinfo/R/4.0.1-foss")

data_directory <-  params$data_directory
analysis_name <- params$analysis_name

1 Results

1.1 Reading Misty results

results_misty_folders <- paste0(data_directory, analysis_name,
  "IntermediaryFiles/Misty_Results/results_TF_Ligands_DorotheaClusters/")

misty_results <- 
  collect_results(as.list(list.dirs(results_misty_folders, recursive = FALSE)))

These are the ligand-TF regulations that we consider:

misty_results_toconsider <-
  misty_results$importances.aggregated[["juxta.ligands_2"]] %>%
  tidyr::pivot_longer(!Predictor, names_to="TFs", values_to  = "importance") %>%
  dplyr::filter(importance > 1) 


misty_results_toconsider %>% 
  dplyr::arrange(desc(importance)) %>% 
  kbl() %>% kable_styling()
Predictor TFs importance
RNF43 SPI1 2.168153
MMP1 FOS 1.927795
RNF43 TBX21 1.878013
MMP1 JUN 1.820045
RNF43 ETV4 1.726244
RNF43 PBX3 1.656090
RNF43 SP3 1.620449
RNF43 ESR1 1.582334
RNF43 LYL1 1.446047
MMP1 NFKB1 1.439213
DCN ETV4 1.342204
RNF43 STAT1 1.320641
RNF43 JUN 1.306373
DCN MEIS1 1.252182
CXCL14 MAF 1.216626
MMP1 HIF1A 1.131988
CXCL14 LYL1 1.131581
CXCL14 IKZF1 1.129121
DCN SPI1 1.126616
PLAU BACH2 1.122926
RNF43 EBF1 1.067074
PLAU NFATC1 1.059674
RNF43 TEAD4 1.052134
THBS2 STAT1 1.044323
LUM RUNX2 1.013957
CXCL14 EBF1 1.005273 
ligands_toconsider <- misty_results_toconsider %>% 
  dplyr::pull(Predictor) %>% unique()

1.2 Reading LIANA results

file_name_toread <- paste0(data_directory, 
  analysis_name, "IntermediaryFiles/MergeClustering/lianaResults.rds")
liana_results <- readRDS(file_name_toread)

Ligand receptor interactions to consider

liana_significant_myligands <- liana_results %>% 
  dplyr::filter(ligand %in% ligands_toconsider) %>%
  dplyr::filter(source %in% c("1", "0"), target %in% c("1")) %>% 
  dplyr::filter(aggregate_rank < 0.01) 

liana_significant_myligands %>%
  # top_n(25, desc(aggregate_rank)) %>%
  liana_dotplot(source_groups = c("0", "1"),
                target_groups = c("1"))

1.3 Building a network using OmniPath

SigNet <- import_omnipath_interactions() %>% 
  dplyr::mutate(Effect = ifelse(consensus_inhibition == "1", "inhibition", "stimulation")) %>% 
  dplyr::select(source_genesymbol, target_genesymbol, Effect)

igraph_SigNet <- graph_from_data_frame(SigNet)

omni_resources <- readRDS(system.file(package = "liana", 
  "omni_resources.rds"))

Ligrec_Net <- omni_resources$OmniPath %>% 
  dplyr::mutate(Effect = ifelse(consensus_inhibition == "1", "inhibition", "stimulation")) %>% 
  dplyr::select(source_genesymbol, target_genesymbol, Effect)

all_edges <- c()

for (current_ligand in ligands_toconsider){
  receptors_currentLigand <-  
    dplyr::filter(liana_significant_myligands, ligand == current_ligand) %>%
    dplyr::pull(receptor) %>% unique()
  TFs_currentLigand <-  
    dplyr::filter(misty_results_toconsider, Predictor == current_ligand) %>% 
    dplyr::pull(TFs) %>% unique()
    
  for (current_receptor in receptors_currentLigand){
    vpath_currentRec <- shortest_paths(
      igraph_SigNet,
      from = current_receptor,
      to = TFs_currentLigand,
      output = c("epath"))  %>% unlist() %>% unique()
      all_edges <- c(all_edges, vpath_currentRec)
  }
}

## Signaling graph
Signaling_df <- subgraph.edges(igraph_SigNet, eids= unique(all_edges), 
          delete.vertices = TRUE) %>% 
  igraph::as_data_frame()


## Ligand receptor graph based on LIANA resutls 
ligrec_df <- liana_significant_myligands %>% 
  dplyr::select(ligand, receptor) %>%
  dplyr::inner_join(Ligrec_Net, by = c("ligand"="source_genesymbol", "receptor"="target_genesymbol")) %>%
  dplyr::rename(from = "ligand", to = "receptor")
LigRec_Sign_df <- rbind.data.frame(ligrec_df, Signaling_df)

1.4 Average expression of network nodes

We compute the average expression of the receptors, signaling intermediates and TFs in the cluster 1 (tumor neighborhood). We therefore need to read Seurat objects

file_seuratobject <- 
  paste0(data_directory,analysis_name, "IntermediaryFiles/MergeClustering/SeuratObject_dorothea_res05.rds")

seurat_merge_obj <- readRDS(file_seuratobject)[[2]]

avg_expression_values <- 
  AverageExpression(seurat_merge_obj,
                  features= unique(c(LigRec_Sign_df$from,LigRec_Sign_df$to)),
                  assays = "SCT") %>%
  as.data.frame() %>% dplyr::select(SCT.1) 

## We add the type of nodes to the avg_expression values. 
avg_expression_values <- avg_expression_values %>% 
  tibble::rownames_to_column(var = "gene") %>% 
  dplyr::mutate(
    NodeType = ifelse(gene %in% ligands_toconsider, "Ligand", 
      ifelse(gene %in% unique(liana_significant_myligands$receptor), "Receptor", 
      ifelse(gene %in% unique(misty_results_toconsider$TFs), "TF", "Signaling Intermediary")))) %>% tibble::column_to_rownames(var = "gene")

## We additionally remove the LR interaction affected by this node removal


LigRec_Sign_df <- 
  dplyr::filter(LigRec_Sign_df, from %in% rownames(avg_expression_values) & 
                  to %in% rownames(avg_expression_values)) %>%
  dplyr::filter(to != "IGF2R")

network_towrite <- paste0(data_directory, 
  analysis_name, "IntermediaryFiles/MergeClustering/NetworkResults/", 
  "network_all.csv")

attributes_towrite <-  paste0(data_directory, 
  analysis_name, "IntermediaryFiles/MergeClustering/NetworkResults/", 
  "attributes_all.csv")


## We remove from the network the gene which is not expressed
write.table(LigRec_Sign_df,file = network_towrite,
  quote = FALSE, row.names = FALSE, col.names = FALSE, sep= ";")

write.table(avg_expression_values, file = attributes_towrite,
    quote = FALSE, row.names = TRUE, col.names = FALSE, sep= ";")

1.5 Import network from Cytoscape