From bf4c12f41c45561bcc2010b5123f0c2c322265c8 Mon Sep 17 00:00:00 2001 From: Goultard59 Date: Thu, 5 May 2022 14:26:17 +0200 Subject: [PATCH 01/17] Create createPPI_OmniPath.R --- R/createPPI_OmniPath.R | 51 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 51 insertions(+) create mode 100644 R/createPPI_OmniPath.R diff --git a/R/createPPI_OmniPath.R b/R/createPPI_OmniPath.R new file mode 100644 index 0000000..f80842a --- /dev/null +++ b/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. +#' +#' @param filename character; protein interaction data. +#' @param expr dgCMatrix; gene expression in each cell, found in +#' CellRouter@@ndata. +#' +#' @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 +createPPI <- function(filename, expr){ + a <- expr + data = read.csv(filename, sep='\t', header=FALSE) + cat(dim(data)) + cat('\nProtein interaction data loaded!') + # The columns names are usually fixed. + ppi <- data[, c('V5','V6', 'V15')] #geneA and geneB, score + ppi <- as.data.frame(apply(ppi, 2, function(x)gsub('\\s+', '',x))) + geneA <- as.vector(unlist(lapply(strsplit(as.character(ppi[,'V5']), + split="\\:"), "[", 2))) + geneB <- as.vector(unlist(lapply(strsplit(as.character(ppi[,'V6']), + split="\\:"), "[", 2))) + score <- as.vector(unlist(lapply(strsplit(as.character(ppi[,'V15']), + split="\\:"), "[", 2))) + score <- as.numeric(score) + geneA <- as.vector(unlist(lapply(strsplit(as.character(geneA), + split="\\|"), "[", 1))) + geneB <- as.vector(unlist(lapply(strsplit(as.character(geneB), + split="\\|"), "[", 1))) + idmap <- data.frame(idA=geneA, idB=geneB, score=score) + 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) +} From e5b04c0ea9c3b1bb903af6db69638d44872f5f05 Mon Sep 17 00:00:00 2001 From: Goultard59 Date: Thu, 5 May 2022 14:36:35 +0200 Subject: [PATCH 02/17] Update DESCRIPTION --- DESCRIPTION | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/DESCRIPTION b/DESCRIPTION index 20f4e32..af1876d 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -50,7 +50,8 @@ Imports: gridExtra, Rcpp, RcppEigen, - sparseMatrixStats + sparseMatrixStats, + OmnipathR SystemRequirements: Java JDK Oracle 1.8 RoxygenNote: 7.1.2 Collate: From 8fffcf44f811b6884d23c66f681c3d6646af4b71 Mon Sep 17 00:00:00 2001 From: Goultard59 Date: Thu, 5 May 2022 14:44:12 +0200 Subject: [PATCH 03/17] Update createPPI_OmniPath.R --- R/createPPI_OmniPath.R | 31 ++++++++++++++----------------- 1 file changed, 14 insertions(+), 17 deletions(-) diff --git a/R/createPPI_OmniPath.R b/R/createPPI_OmniPath.R index f80842a..f9368f5 100644 --- a/R/createPPI_OmniPath.R +++ b/R/createPPI_OmniPath.R @@ -14,26 +14,23 @@ #' (character vector) with the genes not used in the interactions. #' #' @export -createPPI <- function(filename, expr){ +createPPI <- function(expr, species){ a <- expr - data = read.csv(filename, sep='\t', header=FALSE) - cat(dim(data)) + 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)) + cat(dim(ppi)) cat('\nProtein interaction data loaded!') # The columns names are usually fixed. - ppi <- data[, c('V5','V6', 'V15')] #geneA and geneB, score - ppi <- as.data.frame(apply(ppi, 2, function(x)gsub('\\s+', '',x))) - geneA <- as.vector(unlist(lapply(strsplit(as.character(ppi[,'V5']), - split="\\:"), "[", 2))) - geneB <- as.vector(unlist(lapply(strsplit(as.character(ppi[,'V6']), - split="\\:"), "[", 2))) - score <- as.vector(unlist(lapply(strsplit(as.character(ppi[,'V15']), - split="\\:"), "[", 2))) - score <- as.numeric(score) - geneA <- as.vector(unlist(lapply(strsplit(as.character(geneA), - split="\\|"), "[", 1))) - geneB <- as.vector(unlist(lapply(strsplit(as.character(geneB), - split="\\|"), "[", 1))) - idmap <- data.frame(idA=geneA, idB=geneB, score=score) + 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, From 634b28cf628b21ccd2765de0cfad1ccbadd4738b Mon Sep 17 00:00:00 2001 From: Goultard59 Date: Thu, 5 May 2022 15:39:50 +0200 Subject: [PATCH 04/17] Update cellnetworksPPI_CellComm.R --- R/cellnetworksPPI_CellComm.R | 40 ++++++++++++++++++------------------ 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/R/cellnetworksPPI_CellComm.R b/R/cellnetworksPPI_CellComm.R index b46bf7d..9135255 100755 --- a/R/cellnetworksPPI_CellComm.R +++ b/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) } From 188c0cb8c1ce610dcf4cea2e2b0d5fc8176dc86a Mon Sep 17 00:00:00 2001 From: Goultard59 Date: Thu, 5 May 2022 15:46:15 +0200 Subject: [PATCH 05/17] Update findpaths.simpleRJava_CellComm.R --- R/findpaths.simpleRJava_CellComm.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/findpaths.simpleRJava_CellComm.R b/R/findpaths.simpleRJava_CellComm.R index 093e7c1..4de85ae 100755 --- a/R/findpaths.simpleRJava_CellComm.R +++ b/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. From 84b20e48b328dac8cd7117dd993c2ff248bfcfa3 Mon Sep 17 00:00:00 2001 From: Goultard59 Date: Thu, 5 May 2022 16:01:46 +0200 Subject: [PATCH 06/17] Update computeValue.R --- R/computeValue.R | 25 +++++++++++++++++-------- 1 file changed, 17 insertions(+), 8 deletions(-) diff --git a/R/computeValue.R b/R/computeValue.R index c15a45d..6d55ce5 100755 --- a/R/computeValue.R +++ b/R/computeValue.R @@ -17,20 +17,24 @@ #' @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') + } 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'] @@ -85,15 +89,18 @@ 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)] @@ -101,7 +108,9 @@ setMethod("computeValueSubclusters", 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] From 2e619bba15ab1b9dec956e5cc9746e7f4c7b3cf2 Mon Sep 17 00:00:00 2001 From: Goultard59 Date: Thu, 5 May 2022 16:05:24 +0200 Subject: [PATCH 07/17] Update clusterPermutation_CellComm.R --- R/clusterPermutation_CellComm.R | 16 ++++++++++------ 1 file changed, 10 insertions(+), 6 deletions(-) diff --git a/R/clusterPermutation_CellComm.R b/R/clusterPermutation_CellComm.R index d034349..a15d0fd 100755 --- a/R/clusterPermutation_CellComm.R +++ b/R/clusterPermutation_CellComm.R @@ -17,16 +17,18 @@ #' #' @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 @@ -76,7 +78,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 +87,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 +98,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 From 2c12a71084ceccab65d6ac944e9c1257a5fe6dfe Mon Sep 17 00:00:00 2001 From: Goultard59 Date: Thu, 5 May 2022 16:22:00 +0200 Subject: [PATCH 08/17] Update computeValue.R --- R/computeValue.R | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/R/computeValue.R b/R/computeValue.R index 6d55ce5..2ffd2e8 100755 --- a/R/computeValue.R +++ b/R/computeValue.R @@ -27,6 +27,10 @@ setMethod("computeValue", 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]) From 4593093503964c0d611974f59b7945a3d6615161 Mon Sep 17 00:00:00 2001 From: Goultard59 Date: Thu, 5 May 2022 16:24:47 +0200 Subject: [PATCH 09/17] Create Omnipath_LR_network.R --- R/Omnipath_LR_network.R | 48 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) create mode 100644 R/Omnipath_LR_network.R diff --git a/R/Omnipath_LR_network.R b/R/Omnipath_LR_network.R new file mode 100644 index 0000000..f9368f5 --- /dev/null +++ b/R/Omnipath_LR_network.R @@ -0,0 +1,48 @@ +#' Create protein-protein interactions. +#' +#' Create a protein-protein interaction network for the reconstruction of cell +#' signaling pathways. +#' +#' @param filename character; protein interaction data. +#' @param expr dgCMatrix; gene expression in each cell, found in +#' CellRouter@@ndata. +#' +#' @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 +createPPI <- function(expr, species){ + 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)) + 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) +} From 3ebced24a90c17ecf7fe2bb049ad9f14fb47e666 Mon Sep 17 00:00:00 2001 From: Goultard59 Date: Thu, 5 May 2022 16:29:04 +0200 Subject: [PATCH 10/17] Update Omnipath_LR_network.R --- R/Omnipath_LR_network.R | 35 ++++------------------------------- 1 file changed, 4 insertions(+), 31 deletions(-) diff --git a/R/Omnipath_LR_network.R b/R/Omnipath_LR_network.R index f9368f5..f7160ae 100644 --- a/R/Omnipath_LR_network.R +++ b/R/Omnipath_LR_network.R @@ -14,35 +14,8 @@ #' (character vector) with the genes not used in the interactions. #' #' @export -createPPI <- function(expr, species){ - 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)) - 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) +Omnipath_LR_network <- function(species){ + pairs <- nichenet_lr_network(only_omnipath = TRUE) + pairs$Pair.Name <- paste(pairs$from, pairs$to, sep = "_") + return(pairs) } From d748c968c8efe1610c63bb8f25653d8ba0232812 Mon Sep 17 00:00:00 2001 From: Goultard59 Date: Thu, 5 May 2022 16:38:23 +0200 Subject: [PATCH 11/17] Create Path_finding.R --- R/Path_finding.R | 46 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 46 insertions(+) create mode 100644 R/Path_finding.R diff --git a/R/Path_finding.R b/R/Path_finding.R new file mode 100644 index 0000000..c607e48 --- /dev/null +++ b/R/Path_finding.R @@ -0,0 +1,46 @@ +#' Compute cell-cell interaction network. +#' +#' Analyze the mean expression data generated by the computeValue method of +#' CellRouter to determine co-expression patterns of ligands and receptors. +#' +#' @param mean.expr list; the statistics of the expressed genes as calculated by +#' the computeValue method of CellRouter. +#' @param ligands vector; list of genes annotated as ligands. +#' @param receptors vector; list of genes annotated as receptors. +#' @param threshold numeric; threshold to select the genes in the population +#' according to their mean expression. +#' @param pairs the pairs; +#' +#' @return data frame; interactions with celltype1, celltype2, pair, ligand, +#' and receptor. +#' +#' @export +population.pairing <- function(mean.expr, ligands, receptors, threshold, pairs){ + interactions <- list() + for(ct1 in names(mean.expr)){ + for(ct2 in names(mean.expr)){ + ct1.expressed <- mean.expr[[ct1]][c(ligands, receptors),] + ct1.expressed <- ct1.expressed[which(ct1.expressed$p > threshold),] + + ct2.expressed <- mean.expr[[ct2]][c(ligands, receptors),] + ct2.expressed <- ct2.expressed[which(ct2.expressed$p > threshold),] + + combs <- expand.grid(rownames(ct1.expressed), rownames(ct2.expressed)) + combs <- paste(combs$Var1, combs$Var2, sep='_') + combs <- intersect(combs, pairs$Pair.Name) + if(length(combs) > 0){ + ligs <- sapply(strsplit(combs, split='_', fixed=TRUE), function(x) (x[1])) + recs <- sapply(strsplit(combs, split='_', fixed=TRUE), function(x) (x[2])) + df <- data.frame(celltype1=rep(ct1, times=length(combs)), + celltype2=rep(ct2, times=length(combs)), + pair=combs, ligand=ligs, receptor=recs, + stringsAsFactors = FALSE) + interactions[[paste(ct1,ct2,sep='_')]] <- df + } + + } + } + interactions <- dplyr::bind_rows(interactions)#rbind.fill(interactions) + + return(interactions) +} From 30cd82f58178a376b9e60b04e1d3603d46e39a74 Mon Sep 17 00:00:00 2001 From: Goultard59 Date: Thu, 5 May 2022 16:50:48 +0200 Subject: [PATCH 12/17] Update Path_finding.R --- R/Path_finding.R | 62 ++++++++++++++++++++++++++++++------------------ 1 file changed, 39 insertions(+), 23 deletions(-) diff --git a/R/Path_finding.R b/R/Path_finding.R index c607e48..70fda92 100644 --- a/R/Path_finding.R +++ b/R/Path_finding.R @@ -15,32 +15,48 @@ #' and receptor. #' #' @export -population.pairing <- function(mean.expr, ligands, receptors, threshold, pairs){ - interactions <- list() - for(ct1 in names(mean.expr)){ - for(ct2 in names(mean.expr)){ - ct1.expressed <- mean.expr[[ct1]][c(ligands, receptors),] - ct1.expressed <- ct1.expressed[which(ct1.expressed$p > threshold),] +pathfinding <- function(cellrouter, assay.type='ST', grn, gene.pairs_01, gene.pairs_02, markers, maindir=NULL, verbose=TRUE){ + tfs <- intersect(grn.data$tfs, rownames(slot(object, 'assays')[[assay.type]]@ndata)) + x <- rbind(gene.pairs_01$markers, gene.pairs_02$markers) - ct2.expressed <- mean.expr[[ct2]][c(ligands, receptors),] - ct2.expressed <- ct2.expressed[which(ct2.expressed$p > threshold),] + sts2 <- list() + for(i in unique(x$celltypes)){ + tmp <- x[which(x$celltypes == i),] + rs <- unique(tmp$receptor) + sts2[["receptor"]][[i]]$sources <- rs # receptors + sts2[["receptor"]][[i]]$targets <- tfs # TFs (all) + } + if (maindir == NULL) { + wd_package <- getwd() + setwd(wd_package) + maindir <- paste0(wd_package, '') + } + + interactionList <- expand.grid(x$celltypes) - combs <- expand.grid(rownames(ct1.expressed), rownames(ct2.expressed)) - combs <- paste(combs$Var1, combs$Var2, sep='_') - combs <- intersect(combs, pairs$Pair.Name) - if(length(combs) > 0){ - ligs <- sapply(strsplit(combs, split='_', fixed=TRUE), function(x) (x[1])) - recs <- sapply(strsplit(combs, split='_', fixed=TRUE), function(x) (x[2])) - df <- data.frame(celltype1=rep(ct1, times=length(combs)), - celltype2=rep(ct2, times=length(combs)), - pair=combs, ligand=ligs, receptor=recs, - stringsAsFactors = FALSE) - interactions[[paste(ct1,ct2,sep='_')]] <- df - } + cellcomm <- CreateCellComm() + for(interaction in interactionList){ + interaction <- gsub(" ", ".", interaction) + print(interaction) + celltype <- strsplit(interaction, split = "[_]")[[1]][2] + + cat("Using ", celltype, " coexpression network\n") + filename <- paste(celltype, ".txt", sep='') + dirname <- gsub(("[.]"), ".", interaction) + cellcomm <- findpaths.simpleRJava(cellcomm, sources.targets = sts2$receptor[[interaction]], file = filename, dir = dirname, maindir = maindir) + } - } + files <- list() + for(i in interactionList){ + f <- paste(strsplit(i, split = "_")[[1]][2], ".txt", sep = "") + files[[i]] <- f } - interactions <- dplyr::bind_rows(interactions)#rbind.fill(interactions) - return(interactions) + summary <- summarize.flow(files, prefix = "") + colnames(grn_table) <- c("target", "reg", "zscore", "corr") + sampTab <- slot(object, 'assays')[[assay.type]]@sampTab + apathways <- activeSignaling(data = slot(object, 'assays')[[assay.type]]@ndata, sampTab = slot(object, 'assays')[[assay.type]]@sampTab, + grn = grn$GRN_table, markers = markers, summary = summary, fc = 0.2) + npaths <- rankpaths(summary, apathways$pathways) + return(npaths) } From 6fbc81274b7898f2cb092591491c31a4bec2f4f6 Mon Sep 17 00:00:00 2001 From: Goultard59 Date: Fri, 6 May 2022 10:48:43 +0200 Subject: [PATCH 13/17] Update Omnipath_LR_network.R --- R/Omnipath_LR_network.R | 19 ++++--------------- 1 file changed, 4 insertions(+), 15 deletions(-) diff --git a/R/Omnipath_LR_network.R b/R/Omnipath_LR_network.R index f7160ae..fe68129 100644 --- a/R/Omnipath_LR_network.R +++ b/R/Omnipath_LR_network.R @@ -1,21 +1,10 @@ -#' Create protein-protein interactions. +#' Builds ligand-receptor network prior knowledge for fusca using Omnipath. #' -#' Create a protein-protein interaction network for the reconstruction of cell -#' signaling pathways. -#' -#' @param filename character; protein interaction data. -#' @param expr dgCMatrix; gene expression in each cell, found in -#' CellRouter@@ndata. -#' -#' @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. +#' @return data frame containing ligand and receptor. #' #' @export -Omnipath_LR_network <- function(species){ - pairs <- nichenet_lr_network(only_omnipath = TRUE) +Omnipath_LR_network <- function(){ + pairs <- nichenet_lr_network(ramilowski = NULL) pairs$Pair.Name <- paste(pairs$from, pairs$to, sep = "_") return(pairs) } From f265e1c948ce48afc5dc4a344f1232347edfdee4 Mon Sep 17 00:00:00 2001 From: Goultard59 Date: Fri, 6 May 2022 10:56:38 +0200 Subject: [PATCH 14/17] Update clusterPermutation_CellComm.R --- R/clusterPermutation_CellComm.R | 2 ++ 1 file changed, 2 insertions(+) diff --git a/R/clusterPermutation_CellComm.R b/R/clusterPermutation_CellComm.R index a15d0fd..f17d19b 100755 --- a/R/clusterPermutation_CellComm.R +++ b/R/clusterPermutation_CellComm.R @@ -11,6 +11,7 @@ #' 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. @@ -70,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. From 07159a58ca00988af6fbf0f11a15b1062fa1d517 Mon Sep 17 00:00:00 2001 From: Goultard59 Date: Fri, 6 May 2022 10:57:59 +0200 Subject: [PATCH 15/17] Update computeValue.R --- R/computeValue.R | 2 ++ 1 file changed, 2 insertions(+) diff --git a/R/computeValue.R b/R/computeValue.R index 2ffd2e8..182166c 100755 --- a/R/computeValue.R +++ b/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). @@ -81,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). From eb0f23f5aef64e9fb43ef1706872d648f461fc09 Mon Sep 17 00:00:00 2001 From: Goultard59 Date: Fri, 6 May 2022 12:18:07 +0200 Subject: [PATCH 16/17] Update createPPI_OmniPath.R --- R/createPPI_OmniPath.R | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/R/createPPI_OmniPath.R b/R/createPPI_OmniPath.R index f9368f5..fdb1d31 100644 --- a/R/createPPI_OmniPath.R +++ b/R/createPPI_OmniPath.R @@ -1,11 +1,12 @@ #' Create protein-protein interactions. #' #' Create a protein-protein interaction network for the reconstruction of cell -#' signaling pathways. +#' signaling pathways from OmniPath. #' -#' @param filename character; protein interaction data. #' @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 @@ -14,7 +15,7 @@ #' (character vector) with the genes not used in the interactions. #' #' @export -createPPI <- function(expr, species){ +createPPIOmniPath <- function(expr, species = "Hs", verbose = TRUE){ a <- expr if (species = "Hs") { organism <- 9606 @@ -26,8 +27,10 @@ createPPI <- function(expr, species){ organism <- 10116 } data <- as.data.frame(import_omnipath_interactions(datasets = 'omnipath', entity_types = 'protein', organism = organism)) - cat(dim(ppi)) - cat('\nProtein interaction data loaded!') + 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)) From 156242fe1afff5e849a4d8de6b86baaa38222a9d Mon Sep 17 00:00:00 2001 From: Goultard59 Date: Fri, 6 May 2022 12:27:44 +0200 Subject: [PATCH 17/17] Delete Path_finding.R --- R/Path_finding.R | 62 ------------------------------------------------ 1 file changed, 62 deletions(-) delete mode 100644 R/Path_finding.R diff --git a/R/Path_finding.R b/R/Path_finding.R deleted file mode 100644 index 70fda92..0000000 --- a/R/Path_finding.R +++ /dev/null @@ -1,62 +0,0 @@ -#' Compute cell-cell interaction network. -#' -#' Analyze the mean expression data generated by the computeValue method of -#' CellRouter to determine co-expression patterns of ligands and receptors. -#' -#' @param mean.expr list; the statistics of the expressed genes as calculated by -#' the computeValue method of CellRouter. -#' @param ligands vector; list of genes annotated as ligands. -#' @param receptors vector; list of genes annotated as receptors. -#' @param threshold numeric; threshold to select the genes in the population -#' according to their mean expression. -#' @param pairs the pairs; -#' -#' @return data frame; interactions with celltype1, celltype2, pair, ligand, -#' and receptor. -#' -#' @export -pathfinding <- function(cellrouter, assay.type='ST', grn, gene.pairs_01, gene.pairs_02, markers, maindir=NULL, verbose=TRUE){ - tfs <- intersect(grn.data$tfs, rownames(slot(object, 'assays')[[assay.type]]@ndata)) - x <- rbind(gene.pairs_01$markers, gene.pairs_02$markers) - - sts2 <- list() - for(i in unique(x$celltypes)){ - tmp <- x[which(x$celltypes == i),] - rs <- unique(tmp$receptor) - sts2[["receptor"]][[i]]$sources <- rs # receptors - sts2[["receptor"]][[i]]$targets <- tfs # TFs (all) - } - if (maindir == NULL) { - wd_package <- getwd() - setwd(wd_package) - maindir <- paste0(wd_package, '') - } - - interactionList <- expand.grid(x$celltypes) - - cellcomm <- CreateCellComm() - for(interaction in interactionList){ - interaction <- gsub(" ", ".", interaction) - print(interaction) - celltype <- strsplit(interaction, split = "[_]")[[1]][2] - - cat("Using ", celltype, " coexpression network\n") - filename <- paste(celltype, ".txt", sep='') - dirname <- gsub(("[.]"), ".", interaction) - cellcomm <- findpaths.simpleRJava(cellcomm, sources.targets = sts2$receptor[[interaction]], file = filename, dir = dirname, maindir = maindir) - } - - files <- list() - for(i in interactionList){ - f <- paste(strsplit(i, split = "_")[[1]][2], ".txt", sep = "") - files[[i]] <- f - } - - summary <- summarize.flow(files, prefix = "") - colnames(grn_table) <- c("target", "reg", "zscore", "corr") - sampTab <- slot(object, 'assays')[[assay.type]]@sampTab - apathways <- activeSignaling(data = slot(object, 'assays')[[assay.type]]@ndata, sampTab = slot(object, 'assays')[[assay.type]]@sampTab, - grn = grn$GRN_table, markers = markers, summary = summary, fc = 0.2) - npaths <- rankpaths(summary, apathways$pathways) - return(npaths) -}