Add Omnipath as source for CellComm

DESCRIPTION

@@ -50,7 +50,8 @@ Imports:
     gridExtra,
     Rcpp,
     RcppEigen,
-    sparseMatrixStats
+    sparseMatrixStats,
+    OmnipathR
 SystemRequirements: Java JDK Oracle 1.8
 RoxygenNote: 7.1.2
 Collate: 

R/Omnipath_LR_network.R

@@ -0,0 +1,10 @@
+#' Builds ligand-receptor network prior knowledge for fusca using Omnipath.
+#'
+#' @return data frame containing ligand and receptor.
+#'
+#' @export
+Omnipath_LR_network <- function(){
+  pairs <- nichenet_lr_network(ramilowski = NULL)
+  pairs$Pair.Name <- paste(pairs$from, pairs$to, sep = "_")
+  return(pairs)
+}

R/cellnetworksPPI_CellComm.R

@@ -116,32 +116,32 @@ cellnetworksPPI2 <- function(ppi, expr, samples, corThr1, corThr2, column,
                                                  weighted = 'correlation')
     cor.m <- igraph::as_data_frame(cor_g, 'edges')
     rownames(cor.m) <- paste(as.vector(cor.m$from), as.vector(cor.m$to), sep='_')
-    #cor_g <- graph_from_adjacency_matrix(cor, mode='undirected')
-    #cor.m <- melt(cor)
-    #network1 <- igraph::graph.data.frame(cor.m, directed=FALSE);
     #overlap of edges
-    #tmp.net <- (cor_g %s% ppi)
     el <- igraph::get.edgelist(ppi)
     rownames(el) <- paste(el[,1], el[,2], sep='_')
     keep <- intersect(rownames(cor.m), rownames(el))
     cor.m2 <- cor.m[keep,]
-    #cor.m2 <- cor.m2[which(cor.m2$value > 0),]
-    #cor.m2 <- cor.m2[which(cor.m2$correlation > 0),]
     cor.m2 <- cor.m2[which(cor.m2$correlation > corThr2),]
-    network <- igraph::graph.data.frame(cor.m2, directed=FALSE);
-    igraph::E(network)$weight  <- as.vector(cor.m2$correlation)
-    network <- igraph::simplify(network, remove.multiple = TRUE,
-                                remove.loops = TRUE)
-    edges <- as.data.frame(igraph::get.edgelist(network))
-    edges$weight <- abs(igraph::E(network)$weight)
-    edges$corr <- igraph::E(network)$weight
-    celltype <- gsub(" ", ".", celltype)
-    filename <- paste(dir.prefix, celltype, '.txt',sep='')
-    write.table(edges, file=filename, sep='\t', row.names=FALSE,
-                col.names = FALSE, quote=FALSE) #input network
-    gc()
-    corTables[[celltype]] <- network
-    gc(verbose=TRUE)
+    if (nrow(cor.m2) == 0) {
+        print(paste0("No correlation found for ", celltype, " cell type."))
+        gc()
+    }
+    else {
+        network <- igraph::graph.data.frame(cor.m2, directed=FALSE);
+        igraph::E(network)$weight  <- as.vector(cor.m2$correlation)
+        network <- igraph::simplify(network, remove.multiple = TRUE,
+                                    remove.loops = TRUE)
+        edges <- as.data.frame(igraph::get.edgelist(network))
+        edges$weight <- abs(igraph::E(network)$weight)
+        edges$corr <- igraph::E(network)$weight
+        celltype <- gsub(" ", ".", celltype)
+        filename <- paste(dir.prefix, celltype, '.txt',sep='')
+        write.table(edges, file=filename, sep='\t', row.names=FALSE,
+                    col.names = FALSE, quote=FALSE) #input network
+        gc()
+        corTables[[celltype]] <- network
+        gc(verbose=TRUE)
+    }
   }
   return(corTables)
 }

R/clusterPermutation_CellComm.R

@@ -11,22 +11,25 @@
 #' population.pairing function or the calculateObservedMean function.
 #' @param cluster.label character; column in the metadata table
 #' (CellRouter@@sampTab) corresponding to the population.
+#' @param verbose boolean; verbosity of the function.
 #'
 #' @return list; the interaction data frame for each permutation, including the
 #' columns of mean expression of ligands, receptors and pairs.
 #'
 #' @export
 clusterPermutation <- function(cellrouter, assay.type='RNA', genelist,
-                               nPerm, interactions, cluster.label){
+                               nPerm, interactions, cluster.label, verbose=TRUE){
   pcellrouter <- cellrouter
   mean.pairs <- list()
   for(j in 1:nPerm){
-    cat(j, '______________________________\n')
+    if(verbose == TRUE){
+      cat(j, '______________________________\n')
+    }
     pclusters <- as.vector(slot(pcellrouter, 'assays')[[assay.type]]@sampTab[[cluster.label]])
     pclusters <- pclusters[sample(length(pclusters))]
     slot(pcellrouter, 'assays')[[assay.type]]@sampTab[[cluster.label]] <- pclusters
     mean.expr <- computeValue(pcellrouter, assay.type, genelist, cluster.label,
-                              fun='mean')
+                              fun='mean', verbose)
     # interactions <- population.pairing(mean.expr = mean.expr, ligands=ligands,
     #                                    receptors=receptors, threshold = 0.25, pairs.m)
     interactions2 <- interactions
@@ -68,6 +71,7 @@ clusterPermutation <- function(cellrouter, assay.type='RNA', genelist,
 #' @param subcluster.column character; the name of the column where the
 #' subclustering information will be stored.
 #' @param clusters character; selected clusters.
+#' @param verbose boolean; verbosity of the function.
 #'
 #' @return list; the interaction data frame for each permutation, including the
 #' columns of mean expression of ligands, receptors and pairs.
@@ -76,7 +80,7 @@ clusterPermutation <- function(cellrouter, assay.type='RNA', genelist,
 clusterPermutationSubcluster <- function(cellrouter, assay.type='RNA', genelist,
                                          nPerm, interactions, cluster.label,
                                          subcluster.column='Subpopulation',
-                                         clusters){
+                                         clusters, verbose=TRUE){
   pcellrouter <- cellrouter
   sampTab <- slot(pcellrouter, 'assays')[[assay.type]]@sampTab[
     slot(pcellrouter, 'assays')[[assay.type]]@sampTab[[cluster.label]] %in% clusters, ]
@@ -85,7 +89,9 @@ clusterPermutationSubcluster <- function(cellrouter, assay.type='RNA', genelist,
   slot(pcellrouter, 'assays')[[assay.type]]@ndata <- expDat
   mean.pairs <- list()
   for(j in 1:nPerm){
-    cat(j, '______________________________\n')
+    if(verbose == TRUE){
+      cat(j, '______________________________\n')
+    }
     pclusters <- as.vector(slot(pcellrouter, 'assays')[[assay.type]]@sampTab[[subcluster.column]])
     pclusters <- pclusters[sample(length(pclusters))]
     slot(pcellrouter, 'assays')[[assay.type]]@sampTab[[subcluster.column]] <- pclusters
@@ -94,7 +100,7 @@ clusterPermutationSubcluster <- function(cellrouter, assay.type='RNA', genelist,
                                          genelist = genelist,
                                          column = cluster.label,
                                          subcluster.column = subcluster.column,
-                                         clusters = clusters, fun='mean')
+                                         clusters = clusters, fun='mean', verbose)
     # interactions <- population.pairing(mean.expr = mean.expr, ligands=ligands,
     #                                    receptors=receptors, threshold = 0.25, pairs.m)
     interactions2 <- interactions

R/computeValue.R

@@ -9,6 +9,7 @@
 #' differential expression. For example, if 'population' is specified,
 #' population-specific gene signatures will be identified.
 #' @param fun character; statistical function to summary the gene expression.
+#' @param verbose boolean; verbosity of the function.
 #'
 #' @return list; the statistics of the expressed genes in the population (p
 #' column), not in the population (np), and the percentage (percent).
@@ -17,20 +18,28 @@
 #' @docType methods
 #' @rdname computeValue-methods
 setGeneric("computeValue", function(object, assay.type='RNA',
-                                    genelist, column='population', fun='max')
+                                    genelist, column='population', fun='max', verbose = TRUE)
   standardGeneric("computeValue"))
 #' @rdname computeValue-methods
 #' @aliases computeValue
 setMethod("computeValue",
           signature = "CellRouter",
-          definition = function(object, assay.type, genelist, column, fun){
-            print('discovering subpopulation-specific gene signatures')
+          definition = function(object, assay.type, genelist, column, fun, verbose){
+            if(verbose == TRUE){
+              print('discovering subpopulation-specific gene signatures')
+            }
+            genelist <- unique(intersect(genelist, rownames(slot(object, 'assays')[[assay.type]]@ndata)))
+            if (length(genelist) == 0) {
+              stop("Error no matching between the genelist and the assay")
+            }
             expDat <- slot(object, 'assays')[[assay.type]]@ndata[genelist,]
             membs <- as.vector(slot(object, 'assays')[[assay.type]]@sampTab[[column]])
             membs_df <- as.vector(slot(object, 'assays')[[assay.type]]@sampTab[ , c('sample_id', column), drop=FALSE])
             diffs <- list()
             for(i in unique(membs)){
-              cat('cluster ', i, '\n')
+              if(verbose == TRUE){
+                cat('cluster ', i, '\n')
+              }
               if(sum(membs == i) == 0) next
               m_indexes <- membs_df[which(membs_df[[column]] != i), 'sample_id']
               n_indexes <- membs_df[which(membs_df[[column]] == i), 'sample_id']
@@ -73,6 +82,7 @@ detectGenes <- function(expr, min_expr = 0.1){
 #' subclustering information will be stored.
 #' @param clusters character; selected clusters.
 #' @param fun character; statistical function to summary the gene expression.
+#' @param verbose boolean; verbosity of the function.
 #'
 #' @return list; the statistics of the expressed genes in the population (p
 #' column), not in the population (np), and the percentage (percent).
@@ -85,23 +95,28 @@ setGeneric("computeValueSubclusters", function(object, assay.type='RNA',
                                                column='population',
                                                subcluster.column='Subpopulation',
                                                clusters,
-                                               fun='max')
+                                               fun='max',
+                                               verbose = TRUE)
   standardGeneric("computeValueSubclusters"))
 #' @rdname computeValueSubclusters-methods
 #' @aliases computeValueSubclusters
 setMethod("computeValueSubclusters",
           signature = "CellRouter",
           definition = function(object, assay.type, genelist, column,
-                                subcluster.column, clusters, fun){
-            print('discovering subpopulation-specific gene signatures')
+                                subcluster.column, clusters, fun, verbose){
+            if(verbose == TRUE){
+              print('discovering subpopulation-specific gene signatures')
+            }
             sampTab <- slot(object, 'assays')[[assay.type]]@sampTab[
               slot(object, 'assays')[[assay.type]]@sampTab[[column]] %in% clusters,]
             expDat <- slot(object, 'assays')[[assay.type]]@ndata[genelist, rownames(sampTab)]
             #membs <- as.vector(object@sampTab$population)
             membs <- as.vector(sampTab[[subcluster.column]])
             diffs <- list()
             for(i in unique(membs)){
-              cat('cluster ', i, '\n')
+              if(verbose == TRUE){
+                cat('cluster ', i, '\n')
+              }
               if(sum(membs == i) == 0) next
               if(fun == 'max'){
                 m <- if(sum(membs != i) > 1) apply(expDat[, membs != i], 1, max) else expDat[, membs != i]

R/createPPI_OmniPath.R

@@ -0,0 +1,51 @@
+#' Create protein-protein interactions.
+#'
+#' Create a protein-protein interaction network for the reconstruction of cell
+#' signaling pathways from OmniPath.
+#'
+#' @param expr dgCMatrix; gene expression in each cell, found in
+#' CellRouter@@ndata.
+#' @param species character; defining the species: Hs for Homo Sapiens, Mm for Mus Musculus or Rn for Rattus norvegicus.
+#' @param verbose boolean; verbosity of the function.
+#'
+#' @return list; the network (igraph object), the ctable (data frame) with the
+#' edge list, the table (data frame) with the information of the ids of gene A
+#' and B and the scores, the genes in total (character vector), genesPPI
+#' (character vector) with the genes in the interactions, and the remove
+#' (character vector) with the genes not used in the interactions.
+#'
+#' @export
+createPPIOmniPath <- function(expr, species = "Hs", verbose = TRUE){
+  a <- expr
+  if (species = "Hs") {
+    organism <- 9606
+  }
+  else if (species = "Mm") {
+    organism <- 10090
+  }
+  else if (species = "Rn") {
+    organism <- 10116
+  }
+  data <- as.data.frame(import_omnipath_interactions(datasets = 'omnipath', entity_types = 'protein', organism = organism))
+  if(verbose == TRUE){
+    cat(dim(ppi))
+    cat('\nProtein interaction data loaded!')
+  }
+  # The columns names are usually fixed.
+  ppi <- data[, c('source_genesymbol','target_genesymbol', 'curation_effort')] #geneA and geneB, score
+  idmap <- data.frame(idA=ppi$source_genesymbol, idB=ppi$target_genesymbol, score=as.numeric(ppi$curation_effort))
+  allgenes <- unique(c(as.vector(idmap$idA), as.vector(idmap$idB)))
+  inetwork <- igraph::graph.data.frame(idmap, directed=FALSE);
+  inetwork <- igraph::simplify(inetwork, remove.multiple = TRUE,
+                               remove.loops = TRUE)
+  genesPPI <- intersect(rownames(expr), allgenes)
+  remove <- setdiff(igraph::V(inetwork)$name, rownames(expr))
+  inetwork <- igraph::delete.vertices(inetwork, remove)
+  ctable <- as.data.frame(igraph::get.edgelist(inetwork))
+  colnames(ctable) <- c('from', 'to')
+  names <- paste(as.vector(ctable$from), as.vector(ctable$to), sep="_")
+  rownames(ctable) <- names
+  network <- list(network=inetwork, ctable=ctable, table=idmap, genes=allgenes,
+                  genesPPI=genesPPI, removed=remove)
+  return(network)
+}

R/findpaths.simpleRJava_CellComm.R

@@ -127,7 +127,7 @@ setMethod("findpaths.simpleRJava",
                                   'FlowNetwork_all_paths_subnet.gml',
                                   format = 'gml');
             } else {
-              cat('No cells in the path', dir,'. The gaph was not created.\n')
+              cat('No cells in the path', dir,'. The graph was not created.\n')
             }
 
             # Change directory back to the original one.