diff --git a/NAMESPACE b/NAMESPACE
index 5b41425e..d7b39638 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -1,8 +1,11 @@
 # Generated by roxygen2: do not edit by hand
 
 export(add_rl_column)
+export(avg_exp_for_complexes)
+export(avg_reclig_expr)
 export(build_domino)
 export(circos_ligand_receptor)
+export(circos_ligand_receptor_general)
 export(cor_heatmap)
 export(cor_scatter)
 export(count_linkage)
@@ -23,9 +26,14 @@ export(dom_zscores)
 export(feat_heatmap)
 export(gene_network)
 export(incoming_signaling_heatmap)
+export(invert_rec_lig_expr)
+export(invert_rec_lig_linkages)
 export(mean_ligand_expression)
+export(outgoing_network)
 export(plot_differential_linkages)
 export(rename_clusters)
+export(resolve_complexes)
+export(resolve_names)
 export(signaling_heatmap)
 export(signaling_network)
 export(summarize_linkages)
diff --git a/R/class_definitions.R b/R/class_definitions.R
index 41a19b75..d5ec84b7 100644
--- a/R/class_definitions.R
+++ b/R/class_definitions.R
@@ -37,7 +37,9 @@ domino <- methods::setClass(
     clust_de = "matrix",
     misc = "list",
     cl_signaling_matrices = "list",
-    signaling = "matrix"
+    cr_signaling_matrices = "list",
+    signaling = "matrix",
+    rec_signaling = "matrix"
   ),
   prototype = list(
     misc = list("build" = FALSE)
diff --git a/R/convenience_fxns.R b/R/convenience_fxns.R
index 841d4bbc..7be69337 100644
--- a/R/convenience_fxns.R
+++ b/R/convenience_fxns.R
@@ -97,3 +97,105 @@ table_convert_genes <- function(genes, from, to, conversion_table) {
   genesV2 <- cbind(col1[which(col1 %in% genes)], col2[which(col1 %in% genes)])
   return(genesV2)
 }
+
+#' Resolve Names
+#' @description
+#' Return gene names of ligands with non-standard names
+#' @param dom domino object
+#' @param genes vector of gene names to resolve
+#' @details Ligand names (which are stored in dom@linkages$rec_lig) do not always match the gene name
+#' Search all provided names and output the gene name if the ligand name is non-standard
+#' @return vector of length(genes) with applicable values replaced using dom@misc$rl_map information
+#' @export
+#' 
+resolve_names <- function(dom, genes, rec_lig = "lig") {
+  rl_map = dom@misc[["rl_map"]]
+  if (rec_lig == "lig") {
+    gene <- "L.gene"
+    name <- "L.name"
+  } else if (rec_lig == "rec") {
+    gene <- "R.gene"
+    name <- "R.name"
+  } else {
+    stop("rec_lig must be one of 'lig' or 'rec'.\n")
+  }
+  
+  genes_resolved <- sapply(genes, function(l){
+    int <- rl_map[rl_map[[name]] == l, ][1,] 
+    if((int[[name]] != int[[gene]]) & !grepl("\\,", int[[gene]])){
+      int[[gene]]
+    } else { 
+      int[[name]]
+    }
+  })
+  return(genes_resolved)
+} # resolve_names
+
+#' Resolve Complexes
+#' @description
+#' Expand any complex names into their component gene names
+#' @param dom domino object
+#' @param genes vector of complex names to resolve
+#' @details Ligand names (which are stored in dom@linkages$rec_lig) can refer to complexes that are
+#' detailed in dom@linkages$complexes. Search all provided names and if the name is a complex, output component genes, otherwise output the name
+#' @return list of length(genes), with names(list) == genes. List values are the same as the names if not in a complex and are the complex genes if they are in a complex.
+#' @export
+#' 
+resolve_complexes <- function(dom, genes) {
+  genes_list <- lapply(genes, function(l){
+    if(l %in% names(dom@linkages$complexes)){
+      return(dom@linkages$complexes[[l]])
+    } else {
+      return(l)
+    }
+  })
+  names(genes_list) <- genes
+  return(genes_list)
+}
+
+#' Get All Ligands
+#' @description
+#' Get all unique ligands present in dom@linkages$rec_lig
+#' @param dom domino object
+#' @param expressed_only logical indicating whether to subset ligands based on expression in dom@z_scores
+#' @details Get all unique ligands in a domino object, expanding all ligand complexes as well. Optionally subset by expression.
+#' @return vector of ligands if expressed_only = T, list if F
+get_all_reclig <- function(dom, expressed_only = T, rec_lig = "rec") {
+  
+  ### Vector of all ligands expressed
+  if (rec_lig == "rec") {
+    all <- unlist(dom@linkages$rec_lig)
+    resolve_rec_lig <- "lig"
+  } else if (rec_lig == "lig") {
+    if (!"lig_rec" %in% names(dom@linkages)) stop("Must run invert_rec_lig_linkages if rec_lig is set to 'lig'.\n")
+    all <- unlist(dom@linkages$lig_rec)
+    resolve_rec_lig <- "rec"
+  } else {
+    stop("lig_rec must be one of 'lig' or 'rec")
+  }
+  
+  all <- unique(all)
+  all <- all[!all == ""]
+  
+  ### Resolve non-standard ligand names
+  all_names_resolved <- resolve_names(dom, all, rec_lig = resolve_rec_lig)
+  all_names_resolved <- unique(all_names_resolved)
+  
+  ### Resolve complexes
+  if(length(dom@linkages$complexes) > 0){
+    all_complexes_resolved_list <- resolve_complexes(dom, all_names_resolved)
+    all_names_resolved <- unlist(all_complexes_resolved_list)
+  }
+  
+  if (expressed_only) {
+    ### Subset for ligands expressed in the data
+    genes <- intersect(all_names_resolved, rownames(dom@z_scores))
+  } else {
+    genes <- all_names_resolved
+    #genes <- all_complexes_resolved_list
+  } 
+  
+  out_ls <- list("genes" = genes, "complex" = all_complexes_resolved_list)
+  return(out_ls)
+  
+} # get_all_reclig
diff --git a/R/import_fxns.R b/R/import_fxns.R
index 22357ceb..136c21e6 100644
--- a/R/import_fxns.R
+++ b/R/import_fxns.R
@@ -448,6 +448,7 @@ create_domino <- function(
   }
   dom@linkages[["rec_lig"]] <- rec_lig_linkage
   dom@misc[["rl_map"]] <- rl_reading
+  dom <- invert_rec_lig_linkages(dom) # untested
   # Get z-score and cluster info
   if (verbose) {
     message("Getting z_scores, clusters, and counts")
diff --git a/R/plot_fxns.R b/R/plot_fxns.R
index b6183fb7..7b1e00ff 100644
--- a/R/plot_fxns.R
+++ b/R/plot_fxns.R
@@ -79,6 +79,8 @@ signaling_heatmap <- function(
 #' @param clusts Vector of clusters to be included. If NULL then all clusters are used.
 #' @param min_thresh Minimum signaling threshold for plotting. Defaults to -Inf for no threshold.
 #' @param max_thresh Maximum signaling threshold for plotting. Defaults to Inf for no threshold.
+#' @param display_top Number of top receptor-ligand pairs to display in output (rather than all results)
+#' @param display_method How to determine top pairs. Options are 'mean' or 'median' of pair's expression across all incoming populations, or 'max' to take the top N invidivual expression values with no repeats
 #' @param scale How to scale the values (after thresholding). Options are 'none', 'sqrt' for square root, or 'log' for log10.
 #' @param normalize Options to normalize the matrix. Accepted inputs are 'none' for no normalization, 'rec_norm' to normalize to the maximum value with each receptor cluster, or 'lig_norm' to normalize to the maximum value within each ligand cluster 
 #' @param title Either a string to use as the title or a boolean describing whether to include a title. In order to pass the 'main' parameter to  [ComplexHeatmap::Heatmap()]  you must set title to FALSE.
@@ -91,6 +93,7 @@ signaling_heatmap <- function(
 #'
 incoming_signaling_heatmap <- function(
     dom, rec_clust, clusts = NULL, min_thresh = -Inf, max_thresh = Inf,
+    display_top = NULL, display_method = "mean",
     scale = "none", normalize = "none", title = TRUE, ...) {
   if (!dom@misc[["build"]]) {
     stop("Please run domino_build prior to generate signaling network.")
@@ -126,12 +129,34 @@ incoming_signaling_heatmap <- function(
   } else if (normalize != "none") {
     stop("Do not recognize normalize input")
   }
+  # Attempt to get top genes and plot heatmap for only that
+  # Need to figure out color scale first in order for it to work I think.
+  # if (!is.null(display_top)) {
+  #   if (display_method == "mean") {
+  #     top_v <- unique(names(rev(sort(rowMeans(mat, na.rm = T)))[1:display_top]))
+  #   } else if (display_method == "median") {
+  #     top_v <- unique(names(rev(sort(apply(mat, 1, function(x) median(x, na.rm = T))))[1:display_top]))
+  #   } else if (display_method == "max") {
+  #     top_v <- unique(names(rev(sort(apply(mat, 1, function(x) max(x, na.rm = T))))[1:display_top]))
+  #   } else {
+  #     stop("display_method must be one of 'mean', 'median', or 'max'.\n")
+  #   }
+  #   if (length(top_v) < display_top*.5) warning("At least half of the top hits are the same. Consider adjusting min_thresh if you intend to subset ligands")
+  #   mat <- mat[rownames(mat) %in% top_v,]
+  #   #cols <- c("#0000FFFF", "white", "#FF0000FF")
+  #   cols <- colorRampPalette(brewer.pal(11, "RdBu"))(100)
+  #   cols <- determineBreaks(mat, cut_v = 3, verbose_v = F)$colors
+  # } else {
+  #   cols <- NULL
+  # }
+  cols <- NULL
   if (title == TRUE) {
     return(
       Heatmap(
         mat, 
         name = "expression",
         column_title = paste0("Expression of ligands targeting cluster ", rec_clust), 
+        col = cols,
         ...
       )
     )
@@ -140,6 +165,7 @@ incoming_signaling_heatmap <- function(
       Heatmap(
         mat, 
         name = "expression",
+        col = cols,
         ...
       )
     )
@@ -149,6 +175,7 @@ incoming_signaling_heatmap <- function(
         mat,
         name = "expression",
         column_title = title, 
+        col = cols,
         ...
       )
     )
@@ -175,6 +202,7 @@ incoming_signaling_heatmap <- function(
 #' @param scale_by How to size vertices. Options are 'lig_sig' for summed outgoing signaling, 'rec_sig' for summed incoming signaling, and 'none'. In the former two cases the values are scaled with asinh after summing all incoming or outgoing signaling.
 #' @param vert_scale Integer used to scale size of vertices with our without variable scaling from size_verts_by.
 #' @param plot_title Text for the plot's title.
+#' @param offset_sides optional parameter to offset the side labels wihen layout is type 'circle'. Options are 'above' or 'below'
 #' @param ... Other parameters to be passed to plot when used with an `{igraph}` object.
 #' @return an igraph rendered to the active graphics device
 #' @export signaling_network
@@ -190,7 +218,7 @@ incoming_signaling_heatmap <- function(
 signaling_network <- function(
     dom, cols = NULL, edge_weight = 0.3, clusts = NULL, showOutgoingSignalingClusts = NULL,
     showIncomingSignalingClusts = NULL, min_thresh = -Inf, max_thresh = Inf, normalize = "none", scale = "sq",
-    layout = "circle", scale_by = "rec_sig", vert_scale = 3, plot_title = NULL, ...) {
+    layout = "circle", scale_by = "rec_sig", vert_scale = 3, plot_title = NULL, offset_sides = NULL, ...) {
   if (!length(dom@clusters)) {
     stop("This domino object was not built with clusters so there is no intercluster signaling.")
   }
@@ -253,7 +281,7 @@ signaling_network <- function(
   }
   graph <- igraph::graph(links)
   # Get vert colors and scale size if desired.
-  igraph::V(graph)$label.dist <- 1.5
+  igraph::V(graph)$label.dist <- 3.5
   igraph::V(graph)$label.color <- "black"
   v_cols <- cols[names(igraph::V(graph))]
   if (scale_by == "lig_sig" & all(gsub("L_", "", colnames(mat)) %in% names(igraph::V(graph)))) {
@@ -271,6 +299,19 @@ signaling_network <- function(
   if (layout == "circle") {
     v_angles <- 1:length(igraph::V(graph))
     v_angles <- -2 * pi * (v_angles - 1) / length(v_angles)
+    if (!is.null(offset_sides)) {
+      if (offset_sides == "above") {
+        first <- -0.3; middle <- 0.85
+      } else if (offset_sides == "below") {
+        first <- 0.3; middle <- 1.15
+      } else {
+        warning(sprintf("offset_sides argument not NULL, but not 'above' or 'below'. Will be ignored.\n"))
+        first <- v_angles[1]
+        middle <- 1
+      }
+      if (v_angles[1] == 0) v_angles[1] <- first
+      if (length(v_angles) %% 2 == 0) v_angles[(length(v_angles)/2)+1] <- v_angles[(length(v_angles)/2)+1] * middle
+    }
     igraph::V(graph)$label.degree <- v_angles
   }
   names(v_cols) <- c()
@@ -328,7 +369,8 @@ signaling_network <- function(
 #'
 gene_network <- function(dom, clust = NULL, OutgoingSignalingClust = NULL, 
     class_cols = c(lig = "#FF685F",rec = "#47a7ff", feat = "#39C740"),
-    cols = NULL, lig_scale = 1, layout = "grid", ...) {
+    cols = NULL, lig_scale = 1, layout = "grid", 
+    out_name = NULL, in_name = NULL, ...) {
   if (!dom@misc[["build"]]) {
     warning("Please build a signaling network with domino_build prior to plotting.")
   }
@@ -383,7 +425,13 @@ gene_network <- function(dom, clust = NULL, OutgoingSignalingClust = NULL,
     allowed_ligs <- c()
     for (cl in cl_with_signaling) {
       if (!is.null(OutgoingSignalingClust)) {
-        OutgoingSignalingClust <- paste0("L_", OutgoingSignalingClust)
+        OutgoingSignalingClust <- paste0("L_", gsub("^L_", "", OutgoingSignalingClust))
+        temp <- dom@cl_signaling_matrices[[cl]]
+        missing_clust <- setdiff(OutgoingSignalingClust, colnames(temp))
+        if (length(missing_clust)) {
+          warning(sprintf("No signaling for %s", paste0(missing_clust, collapse = "; ")))
+          OutgoingSignalingClust <- intersect(OutgoingSignalingClust, colnames(mat))
+        }
         mat <- dom@cl_signaling_matrices[[cl]][, OutgoingSignalingClust]
         if (is.null(dim(mat))) {
           allowed_ligs <- names(mat[mat > 0])
@@ -430,7 +478,7 @@ gene_network <- function(dom, clust = NULL, OutgoingSignalingClust = NULL,
     all_sums <- all_sums[names(all_sums) %in% names(v_size)]
     v_size[names(all_sums)] <- 0.5 * all_sums * lig_scale
   }
-  names(v_size) <- c()
+  names(v_size) <- c() # ? I guess need them for the lig_scale mapping, but can't have them in graph?
   igraph::V(graph)$size <- v_size
   igraph::V(graph)$label.degree <- pi
   igraph::V(graph)$label.offset <- 2
@@ -460,7 +508,17 @@ gene_network <- function(dom, clust = NULL, OutgoingSignalingClust = NULL,
   } else if (layout == "kk") {
     l <- igraph::layout_with_kk(graph)
   }
-  plot(graph, layout = l, main = paste0("Signaling ", OutgoingSignalingClust, " to ", clust), ...)
+  if (!is.null(out_name)) {
+    main <- paste0("Signaling from ", out_name)
+  } else {
+    main <- paste0("Signaling")
+  }
+  if (!is.null(in_name)) {
+    main <- paste0(main, " to ", in_name)
+  } else {
+    main <- paste0(main, " to ", clust)
+  }
+  plot(graph, layout = l, main = main, ...)
   return(invisible(list(graph = graph, layout = l)))
 }
 
@@ -475,6 +533,7 @@ gene_network <- function(dom, clust = NULL, OutgoingSignalingClust = NULL,
 #' @param title Either a string to use as the title or a boolean describing whether to include a title. In order to pass the 'main' parameter to  [ComplexHeatmap::Heatmap()]  you must set title to FALSE.
 #' @param norm Boolean indicating whether or not to normalize the transcrption factors to their max value.
 #' @param feats Either a vector of features to include in the heatmap or 'all' for all features. If left NULL then the features selected for the signaling network will be shown.
+#' @param clust A vector of clusters in the domino object that should be included in the output. NULL (default) will use all clusters.
 #' @param ann_cols Boolean indicating whether to include cell cluster as a column annotation. Colors can be defined with cols. If FALSE then custom annotations can be passed to NMF.
 #' @param cols Named vector of colors to annotate cells by cluster color. Values are taken as colors and names as cluster. If left as NULL then default ggplot colors will be generated.
 #' @param min_thresh Minimum threshold for color scaling if not a boolean heatmap
@@ -491,7 +550,7 @@ gene_network <- function(dom, clust = NULL, OutgoingSignalingClust = NULL,
 #'   max_thresh = 0.6, norm = TRUE, bool = FALSE)
 #' 
 feat_heatmap <- function(
-    dom, feats = NULL, bool = FALSE, bool_thresh = 0.2, title = TRUE, norm = FALSE,
+    dom, feats = NULL, clust = NULL, bool = FALSE, bool_thresh = 0.2, title = TRUE, norm = FALSE,
     cols = NULL, ann_cols = TRUE, min_thresh = NULL, max_thresh = NULL, ...) {
   if (!length(dom@clusters)) {
     warning("This domino object wasn't built with clusters. Cells will not be ordered.")
@@ -499,7 +558,17 @@ feat_heatmap <- function(
   }
   mat <- dom@features
   cl <- dom@clusters
-  cl <- sort(cl)
+  #cl <- sort(cl)
+  if (!is.null(clust)) {
+    lvls <- levels(cl)
+    cells <- which(as.character(cl) %in% clust)
+    mat <- mat[,cells]
+    cl <- cl[cells]
+    cl <- factor(as.character(cl), levels = intersect(lvls, clust))
+  }
+  names(cl) <- colnames(mat); cl <- sort(cl)
+  #mat <- mat[,match(colnames(mat), names(cl))]
+  mat <- mat[,match(names(cl), colnames(mat))]
   if (norm & (!is.null(min_thresh) | !is.null(max_thresh))) {
     warning("You are using norm with min_thresh and max_thresh. Note that values will be thresholded AFTER normalization.")
   }
@@ -524,6 +593,7 @@ feat_heatmap <- function(
   if (is.null(feats)) {
     feats <- c()
     links <- dom@linkages$clust_tf
+    if (!is.null(clust)) links <- links[clust]
     for (i in links) {
       feats <- c(feats, i)
     }
@@ -539,9 +609,9 @@ feat_heatmap <- function(
   } else if (feats == "all") {
     feats <- rownames(mat)
   }
-  if (length(cl)) {
-    mat <- mat[feats, names(cl)]
-  }
+  # if (length(cl)) {
+  #   mat <- mat[feats, names(cl)]
+  # }
   if (ann_cols) {
     ac <- list(Cluster = cl)
     names(ac[[1]]) <- c()
@@ -550,6 +620,13 @@ feat_heatmap <- function(
       names(cols) <- levels(cl)
     }
     # cols <- list(Cluster = cols)
+    missing_color <- setdiff(levels(cl), names(cols))
+    if (length(missing_color)) {
+      warning(sprintf("%s cluster level(s) missing in provided colors. Using randomly-generated colors for: \n%s\n",
+                      length(missing_color), paste0(missing_color, collapse = "\n ")))
+      new_color <- ggplot_col_gen(length(missing_color)); names(new_color) <- missing_color
+      cols <- c(cols, new_color)
+    }
     feat_anno <- columnAnnotation(
       Cluster = cl,
       col = list(Cluster = cols)
@@ -741,37 +818,124 @@ cor_scatter <- function(dom, tf, rec, remove_rec_dropout = TRUE, ...) {
 #' 
 circos_ligand_receptor <- function(
     dom, receptor, ligand_expression_threshold = 0.01, cell_idents = NULL,
-    cell_colors = NULL) {
-  ligands <- dom@linkages$rec_lig[[receptor]]
+    cell_colors = NULL, title = paste0(receptor, " Signaling")) {
+ 
+  warning("circos_ligand_receptor is deprecated in favor of circos_ligand_receptor_general.\n")
+  
+  circos_ligand_receptor_general(dom = dom, receptor = receptor, ligand = NULL,
+                                 expression_threshold = ligand_expression_threshold, cell_idents = cell_idents, 
+                                 cell_colors = cell_colors, title = title)
+}
+
+#' Plot expression of a receptor's ligands by other cell types as a chord plot
+#'
+#' Creates a chord plot of expression of ligands that can activate a specified
+#' receptor where chord widths correspond to mean ligand expression by the cluster.
+#'
+#' @param dom Domino object that has undergone network building with [build_domino()]
+#' @param receptor Name of a receptor active in at least one cell type in the domino object
+#' @param ligand Name of a ligand active in at least one cell type in the domino object
+#' @param which_v If both receptor and ligand are set, this must be set to either "receptor" or "ligand" to determine which to use.
+#' @param expression_threshold Minimum mean expression value of a ligand by a cell type for a chord to be rendered between the cell type and the receptor
+#' @param cell_idents Vector of cell types from cluster assignments in the domino object to be included in the plot.
+#' @param cell_colors Named vector of color names or hex codes where names correspond to the plotted cell types and the color values
+#' @param title if provided, final title will be paste0(receptor/ligand, " Signaling")
+#' @param label_arcs logical indicating if arc populations should be labeled
+#' @param multi_plot numeric vector that, if present indicates which plot in a list of plots this one is (for legend output)
+#' @return renders a circos plot to the active graphics device
+#' @export circos_ligand_receptor_general
+#' 
+circos_ligand_receptor_general <- function(
+    dom, receptor = NULL, ligand = NULL, which_v = NULL, expression_threshold = 0.01, cell_idents = NULL,
+    cell_colors = NULL, title = "Signaling", label_arcs = T, multi_plot = NULL) {
+  
+  if (!is.null(receptor) & !is.null(ligand)) {
+    
+    if (is.null(which_v)) {
+      
+      stop("Both receptor and ligand are set and neither is specified with which_v. Must choose one or the other.\n")
+      
+    } else if (which_v == "receptor") {
+      
+      arc_genes <- dom@linkages$rec_lig[[receptor]]
+      arc_genes <- intersect(arc_genes, ligand)
+      sig_mat <- "cl"
+      dest <- receptor
+      
+    } else if (which_v == "ligand") {
+      
+      if (!"lig_rec" %in% names(dom@linkages)) {
+        warning("ligand argument set, but dom@linkages$lig_rec doesn't exist. Running internally.\n")
+        dom <- invert_rec_lig_linkages(dom)
+      } # fi
+      arc_genes <- dom@linkages$lig_rec[[ligand]]
+      arc_genes <- intersect(arc_genes, receptor)
+      sig_mat <- "cr"
+      dest <- ligand
+    } # fi
+    
+  } else if (!is.null(receptor)) {
+    
+    arc_genes <- dom@linkages$rec_lig[[receptor]]
+    sig_mat <- "cl"
+    dest <- receptor
+    
+  } else if (!is.null(ligand)) {
+    
+    if (!"lig_rec" %in% names(dom@linkages)) {
+      warning("ligand argument set, but dom@linkages$lig_rec doesn't exist. Running internally.\n")
+      dom <- invert_rec_lig_linkages(dom)
+    } # fi
+    arc_genes <- dom@linkages$lig_rec[[ligand]]
+    sig_mat <- "cr"
+    dest <- ligand
+  } # fi
+  
   signaling_df <- NULL
   if (is.null(cell_idents)) {
     # default to all cluster labels in domino object in alphabetical order
     cell_idents <- sort(unique(dom@clusters))
   }
   # obtain expression values from cl_signaling matrices
-  active_chk <- sapply(dom@linkages$clust_rec, function(x) {
-    receptor %in% x
-  })
+  if (sig_mat == "cl") {
+    active_chk <- sapply(dom@linkages$clust_rec, function(x) {
+      dest %in% x
+    })
+  } else if (sig_mat == "cr") {
+    active_chk <- sapply(dom@linkages$clust_incoming_lig, function(x) {
+      dest %in% x
+    })
+  } 
+  
   if (sum(active_chk)) {
     # obtain a signaling matrix where receptor is active
     active_cell <- names(active_chk[active_chk == TRUE])
-    sig <- dom@cl_signaling_matrices[active_cell][[1]]
-    cell_names <- gsub("^L_", "", colnames(sig))
-    for (l in ligands) {
-      df <- data.frame(origin = paste0(cell_names, "-", l), destination = receptor, mean.expression = unname(sig[rownames(sig) ==
-        l, ]))
+    if (sig_mat == "cl") {
+      sig <- dom@cl_signaling_matrices[active_cell][[1]]
+    } else if (sig_mat == "cr") {
+      sig <- dom@cr_signaling_matrices[active_cell][[1]]
+    }
+    
+    cell_names <- gsub("^L_|^R_", "", colnames(sig))
+    for (a in arc_genes) {
+      if (length(sig[rownames(sig) == a,]) == 0) next # have to add this because dom@linkages$clust_incoming_ligand isn't filtered for expression
+      df <- data.frame(origin = paste0(cell_names, "-", a), destination = dest, mean.expression = unname(sig[rownames(sig) == a, ]))
+      df <- df[grepl(paste0(cell_idents, collapse = "|"), df$origin),]
       signaling_df <- rbind(signaling_df, df)
     }
   } else {
-    stop(paste0("No clusters have active ", receptor, " signaling"))
+    warning(paste0("No clusters have active ", dest, " signaling"))
+    return()
   }
+  arc_genes <- intersect(arc_genes, unique(gsub("^.*-", "", signaling_df$origin)))
   signaling_df$mean.expression[is.na(signaling_df$mean.expression)] <- 0
   # create a scaled mean expression plot for coord widths greater than 1 by dividing by the max
   # expression [range (0-1)] scaled.mean will only be used when the max expression is > 1
   signaling_df$scaled.mean.expression <- signaling_df$mean.expression / max(signaling_df$mean.expression)
   # exit function if no ligands are expressed above ligand expression threshold
-  if (sum(signaling_df[["mean.expression"]] > ligand_expression_threshold) == 0) {
-    stop(paste0("No ligands of ", receptor, " exceed ligand expression threshold."))
+  if (sum(signaling_df[["mean.expression"]] > expression_threshold) == 0) {
+    warning(paste0("No pairs of ", dest, " exceed expression threshold."))
+    return()
   }
   # initialize chord diagram with even ligand arcs
   arc_df <- signaling_df[, c("origin", "destination")]
@@ -779,33 +943,38 @@ circos_ligand_receptor <- function(
   # receptor arc will always sum to 4 no matter how many ligands and cell idents are plotted
   arc_df["receptor.arc"] <- 4 / (nrow(signaling_df))
   # name grouping based on [cell_ident]
-  nm <- c(receptor, arc_df$origin)
-  group <- structure(c(nm[1], gsub("-.*", "", nm[-1])), names = nm)
+  nm <- c(dest, arc_df$origin)
+  sub_nm <- sapply(nm[-1], function(x) {
+    y <- strsplit(x, split = "-")[[1]]
+    y <- y[1:length(y)-1]
+    y <- paste0(y, collapse = "-")}, USE.NAMES = F)
+  group <- structure(c(nm[1], sub_nm), names = nm)
   # order group as a factor with the receptor coming first
   group <- factor(group, levels = c(
-    receptor, sort(unique(gsub("-.*", "", nm))[-1]) # alphabetical order of the other cell idents
+    dest, sort(unique(sub_nm)) # alphabetical order of the other cell idents
   ))
   # colors for ligand chords
-  lig_colors <- ggplot_col_gen(length(ligands))
-  names(lig_colors) <- ligands
+  arc_colors <- ggplot_col_gen(length(arc_genes))
+  names(arc_colors) <- arc_genes
   # colors for [cell_ident] arcs
   if (is.null(cell_colors)) {
     cell_colors <- ggplot_col_gen(length(cell_idents))
     names(cell_colors) <- cell_idents
   }
   grid_col <- c("#FFFFFF") # hide the arc corresponding to the receptor by coloring white
-  for (i in 1:length(ligands)) {
-    grid_col <- c(grid_col, rep(lig_colors[i], length(cell_idents)))
+  for (i in 1:length(arc_colors)) {
+    grid_col <- c(grid_col, rep(arc_colors[i], length(cell_idents)))
   }
-  names(grid_col) <- c(receptor, signaling_df$origin)
+  names(grid_col) <- c(dest, signaling_df$origin)
   circlize::circos.clear()
-  circlize::circos.par(start.degree = 0)
+  circlize::circos.par(start.degree = 0, circle.margin = 0.5)
   circlize::chordDiagram(arc_df,
-    group = group, grid.col = grid_col, link.visible = FALSE, annotationTrack = c("grid"),
-    preAllocateTracks = list(track.height = circlize::mm_h(4), track.margin = c(circlize::mm_h(2), 0)), big.gap = 2
+                         group = group, grid.col = grid_col, link.visible = FALSE, annotationTrack = c("grid"),
+                         preAllocateTracks = list(track.height = circlize::mm_h(4), track.margin = c(circlize::mm_h(2), 0)), big.gap = 2
   )
-  for (send in signaling_df$origin) {
-    if (signaling_df[signaling_df$origin == send, ][["mean.expression"]] > ligand_expression_threshold) {
+  if (!is.null(title)) title(paste0(dest, " ", title))
+  for (send in signaling_df$origin[1:length(grid_col)-1]) {
+    if (signaling_df[signaling_df$origin == send, ][["mean.expression"]] > expression_threshold) {
       if (max(signaling_df[["mean.expression"]]) > 1) {
         expr <- signaling_df[signaling_df$origin == send, ][["scaled.mean.expression"]]
         max_width <- signif(max(signaling_df[["mean.expression"]]), 2)
@@ -813,29 +982,34 @@ circos_ligand_receptor <- function(
         expr <- signaling_df[signaling_df$origin == send, ][["mean.expression"]]
         max_width <- 1
       }
-      circlize::circos.link(send, c(0.5 - (expr / 2), 0.5 + (expr / 2)), receptor, 2, col = paste0(
+      circlize::circos.link(send, c(0.5 - (expr / 2), 0.5 + (expr / 2)), dest, 2, col = paste0(
         grid_col[[send]],
         "88"
       ))
     }
   }
   sector_names <- circlize::get.all.sector.index()
-  cell_sectors <- cell_idents[cell_idents %in% gsub("-.*", "", sector_names)]
+  sub_sector_names <- c(sector_names[1], sapply(sector_names[-1], function(x) {
+    y <- strsplit(x, split = "-")[[1]]
+    y <- y[1:length(y)-1]
+    y <- paste0(y, collapse = "-")}, USE.NAMES = F))
+  cell_sectors <- cell_idents[cell_idents %in% sub_sector_names]
   for (cell in cell_sectors) {
     row_pick <- sector_names[grepl(paste0("^", cell), sector_names)]
     if (length(row_pick)) {
+      if (label_arcs) { cell_text <- cell } else { cell_text <- NULL }
       circlize::highlight.sector(sector_names[grepl(paste0("^", cell, "-"), sector_names)],
-        track.index = 1,
-        col = cell_colors[cell], text = cell, cex = 1, facing = "inside", text.col = "black",
-        niceFacing = FALSE, text.vjust = -1.5
+                                 track.index = 1,
+                                 col = cell_colors[cell], text = cell_text, cex = 1, facing = "inside", text.col = "black",
+                                 niceFacing = FALSE, text.vjust = -1.5
       )
     }
   }
-  # highlight receptor sector
-  circlize::highlight.sector(sector_names[grepl(paste0("^", receptor, "$"), sector_names)],
-    track.index = 1,
-    col = "#FFFFFF", text = receptor, cex = 1.5, facing = "clockwise", text.col = "black", niceFacing = TRUE,
-    pos = 4
+  # highlight sector
+  circlize::highlight.sector(sector_names[grepl(paste0("^", dest, "$"), sector_names)],
+                             track.index = 1,
+                             col = "#FFFFFF", text = dest, cex = 1.5, facing = "clockwise", text.col = "black", niceFacing = TRUE,
+                             pos = 4
   )
   # create legends
   lgd_cells <- ComplexHeatmap::Legend(
@@ -843,20 +1017,36 @@ circos_ligand_receptor <- function(
     title_position = "topleft", title = "cell identity"
   )
   lgd_ligands <- ComplexHeatmap::Legend(
-    at = ligands, type = "grid", legend_gp = grid::gpar(fill = lig_colors), title_position = "topleft",
+    at = arc_genes, type = "grid", legend_gp = grid::gpar(fill = arc_colors), title_position = "topleft",
     title = "ligand"
   )
-  chord_width <- 10 / (4 + length(cell_idents) * length(ligands))
+  chord_width <- 10 / (4 + length(cell_idents) * length(arc_genes))
   lgd_chord <- ComplexHeatmap::Legend(
-    at = c(ligand_expression_threshold, max_width), col_fun = circlize::colorRamp2(c(
-      ligand_expression_threshold,
+    at = c(expression_threshold, max_width), col_fun = circlize::colorRamp2(c(
+      expression_threshold,
       max_width
     ), c("#DDDDDD", "#DDDDDD")), legend_height = grid::unit(chord_width, "in"), title_position = "topleft",
     title = "ligand expression"
   )
-  lgd_list_vertical <- ComplexHeatmap::packLegend(lgd_cells, lgd_ligands, lgd_chord)
-  ComplexHeatmap::draw(lgd_list_vertical, x = grid::unit(0.02, "npc"), y = grid::unit(0.98, "npc"), just = c("left", "top"))
-}
+  if (is.null(multi_plot)) {
+    lgd_list_vertical <- ComplexHeatmap::packLegend(lgd_cells, lgd_ligands, lgd_chord)
+    ComplexHeatmap::draw(lgd_list_vertical, x = grid::unit(0.02, "npc"), y = grid::unit(0.98, "npc"), just = c("left", "top"))
+  } else {
+    lgd_list_vertical <- ComplexHeatmap::packLegend(lgd_cells, lgd_ligands, lgd_chord)
+    #lgd_list_vertical <- lgd_chord
+    if (multi_plot %% 2 == 0) { 
+      ComplexHeatmap::draw(lgd_list_vertical, x = grid::unit(0.9, "npc"), y = grid::unit(0.98-(0.13*(multi_plot-2)), "npc"), just = c("left", "top"))
+    } else { 
+      ComplexHeatmap::draw(lgd_list_vertical, x = grid::unit(0.02, "npc"), y = grid::unit(0.98-(0.13*(multi_plot-1)), "npc"), just = c("left", "top"))
+      # if (multi_plot == 3) {
+      #   ComplexHeatmap::draw(lgd_list_vertical, x = grid::unit(0.02, "npc"), y = grid::unit(0.7, "npc"), just = c("left", "top"))
+      # } else {
+      #   ComplexHeatmap::draw(lgd_list_vertical, x = grid::unit(0.02, "npc"), y = grid::unit(0.7-(0.13*(multi_plot-2)), "npc"), just = c("left", "top"))
+      # }
+    }
+    
+  }
+} # circos_ligand_receptor_general
 
 #' Plot differential linkages among domino results ranked by a comparative statistic
 #'
diff --git a/R/processing_fxns.R b/R/processing_fxns.R
index 3ed45315..4bab99b9 100644
--- a/R/processing_fxns.R
+++ b/R/processing_fxns.R
@@ -107,6 +107,9 @@ build_domino <- function(
       clust_ligs[[clust]] <- vec
     }
     dom@linkages[["clust_incoming_lig"]] <- clust_ligs
+    # In order to build: dom@linkages[["clust_lig"]]
+    # I have to get dom@linkages$clust_tf_lig to know which are expressed
+    # Doesn't matter for now, because I can subset the ligands in my plot.
     # Build signaling matrices for each cluster
     cl_signaling_matrices <- list()
     signaling <- matrix(0, ncol = length(levels(dom@clusters)), nrow = length(levels(dom@clusters)))
@@ -219,3 +222,305 @@ lc <- function(list, list_names) {
   }
   return(vec)
 }
+
+#' Invert Receptor Ligand Data
+#' @description
+#' Reformat exisiting domino object to also record expression data/receptor mapping
+#' from the ligand-receptor direction rather than receptor-ligand
+#' @param dom domino object built by domino build
+#' @return domino object with updated `rec_signaling` slot and `cr_signaling_matrices` slot
+#' @export
+#' 
+invert_rec_lig_expr <- function(dom) {
+  
+  ### Check
+  if (!dom@misc[["build"]]) {
+    stop("Must build a signaling network with domino_build prior to reformatting")
+  }
+  
+  ### Add linkages, if missing
+  if (!"lig_rec" %in% names(dom@linkages)) dom <- invert_rec_lig_linkages(dom)
+
+  ### Grab object stuff
+  clust_ligs <- dom@linkages$clust_incoming_lig
+  clust_recs <- dom@linkages$clust_rec
+  
+  # Build signaling matrices for each cluster 
+  #(This is identical to lines 111-185 in processing_fxns.R/build_domino, just switch all receptor-ligand references)
+  cr_signaling_matrices <- list()
+  signaling <- matrix(0, ncol = length(levels(dom@clusters)), nrow = length(levels(dom@clusters)))
+  rownames(signaling) <- paste0("L_", levels(dom@clusters))
+  colnames(signaling) <- paste0("R_", levels(dom@clusters))
+  for (clust in levels(dom@clusters)) {
+    inc_recs <- clust_recs[[clust]]
+    rl_map <- dom@misc[["rl_map"]]
+    inc_recs <- sapply(inc_recs, function(r) {
+      int <- rl_map[rl_map$R.name == r, ][1, ]
+      if ((int$R.name != int$R.gene) & !grepl("\\,", int$R.gene)) {
+        int$R.gene
+      } else {
+        int$R.name
+      }
+    })
+    if (length(dom@linkages$complexes) > 0) {
+      # if complexes were used
+      inc_recs_list <- lapply(inc_recs, function(r) {
+        if (r %in% names(dom@linkages$complexes)) {
+          return(dom@linkages$complexes[[r]])
+        } else {
+          return(r)
+        }
+      })
+      names(inc_recs_list) <- inc_recs
+      inc_recs <- unlist(inc_recs_list)
+    }
+    rec_genes <- intersect(inc_recs, rownames(dom@z_scores))
+    if (length(rec_genes) %in% c(0, 1)) {
+      rec_genes <- numeric(0)
+    }
+    cr_sig_mat <- matrix(0, ncol = length(levels(dom@clusters)), nrow = length(rec_genes))
+    colnames(cr_sig_mat) <- colnames(signaling)
+    rownames(cr_sig_mat) <- rec_genes
+    for (c2 in levels(dom@clusters)) {
+      n_cell <- length(which(dom@clusters == c2))
+      if (n_cell > 1) {
+        expr <- matrix(dom@z_scores[rec_genes, which(dom@clusters == c2)], nrow = length(rec_genes))
+        sig <- rowMeans(expr)
+      } else if (n_cell == 1) {
+        sig <- dom@z_scores[rec_genes, which(dom@clusters == c2)]
+      } else {
+        sig <- rep(0, length(rec_genes))
+        names(sig) <- rec_genes
+      }
+      # mean scaled expression less than 0 is brought up to 0 as a floor
+      sig[which(sig < 0)] <- 0
+      cr_sig_mat[, paste0("R_", c2)] <- sig
+    }
+    if (length(dom@linkages$complexes) > 0) {
+      # if complexes were used
+      cr_sig_list <- lapply(seq_along(inc_recs_list), function(x) {
+        if (all(inc_recs_list[[x]] %in% rec_genes)) {
+          # Some of the ligands in the list object may not be present in the data
+          if (length(inc_recs_list[[x]]) > 1) {
+            return(colMeans(cr_sig_mat[inc_recs_list[[x]], ]))
+          } else {
+            return(cr_sig_mat[inc_recs_list[[x]], ])
+          }
+        }
+      })
+      names(cr_sig_list) <- names(inc_recs_list)
+      if (length(cr_sig_list) > 1) {
+        cr_sig_mat <- do.call(rbind, cr_sig_list)
+      }
+    }
+    cr_signaling_matrices[[clust]] <- cr_sig_mat
+    if (nrow(cr_sig_mat) > 1) {
+      signaling[paste0("L_", clust), ] <- colSums(cr_sig_mat)
+    } else {
+      signaling[paste0("L_", clust), ] <- 0
+    }
+  }
+  dom@cr_signaling_matrices <- cr_signaling_matrices
+  dom@rec_signaling <- signaling
+  
+  return(dom)
+  
+} # invert_rec_lig_expr
+
+#' Invert Receptor Ligand Linkages
+#' @description
+#' Use dom@misc$rl_map to create inverted version of dom@linkages$rec_lig
+#' @param dom domino object built by domino build
+#' @return domino object with new entry in dom@linkages called `lig_rec`, which is a list
+#' with one element per ligand and the list values are all receptors mapped to that ligand by rl_map
+#' @export
+#' 
+invert_rec_lig_linkages <- function(dom) {
+  
+  ### Check
+  if (!dom@misc[["build"]]) {
+    stop("Must build a signaling network with domino_build prior to reformatting")
+  }
+  
+  rl_map <- dom@misc$rl_map
+  
+  ### Add linkages, if missing
+  if (!"lig_rec" %in% names(dom@linkages)) {
+    
+    lig_genes <- unique(unlist(strsplit(rl_map[["L.gene"]], split = "\\,")))
+    lig_names <- rl_map[["L.name"]]
+    
+    lig_rec_linkage <- list()
+    for (lig in lig_names) {
+      inter <- rl_map[rl_map[["L.name"]] == lig, ]
+      recs <- inter[["R.name"]]
+      lig_rec_linkage[[lig]] <- recs
+    } # for lig
+    
+    dom@linkages[["lig_rec"]] <- lig_rec_linkage
+    
+  } else {
+    message("dom@linkages$lig_rec already exists.\n")
+  } # fi
+  
+  return(dom)
+  
+} # invertRecLigLinkages
+
+#' Average Complex Expression
+#' @description
+#' Average complex expression
+#' @param exp_mat expression matrix of ligands that may or may not be in a complex
+#' @param complexes_list named list. names = ligands, elements = ligands and/or complex compononents
+#' @details Search each element of complexes_list and determine if actually part of a complex. 
+#' If it is, calculate colMeans for all ligands in the complex. If not, return the ligand expression.
+#' rewrite with purrr?
+#' @return list
+#' @export
+#' 
+avg_exp_for_complexes <- function(exp_mat, complexes_list) {
+  trim_list <- complexes_list %>% purrr::keep(~{all(.x %in% rownames(exp_mat))})
+  gene_exp_list <- lapply(seq_along(trim_list), function(x) {
+    if (length(trim_list[[x]]) > 1) {
+      return(colMeans(exp_mat[trim_list[[x]], ,drop=F]))
+    } else {
+      return(exp_mat[trim_list[[x]], ])
+    }
+  })
+  names(gene_exp_list) <- names(trim_list)
+  #if(length(gene_exp_list)>1){mat <- do.call(rbind, gene_exp_list)}
+  return(gene_exp_list)
+}
+
+avg_reclig_expr <- function(dom, cluster = NULL, genes, complexes) {
+  #' Average Receptor/Ligand Expression
+  #' @description
+  #' Create a matrix of average expression of cells in provided clsuters
+  #' @param dom domino object
+  #' @param cluster vector of clusters to include (NULL uses all)
+  #' @param genes genes to get averages of
+  #' @param complexes list of resolved complexes
+  #' @details
+  #' foo
+  #' @return foo
+  #' @export
+  #'
+  
+  if (is.null(cluster)) cluster <- levels(dom@clusters)
+  
+  cl_expr <- matrix(0, ncol = length(cluster), nrow = length(genes))
+  colnames(cl_expr) <- cluster
+  rownames(cl_expr) <- genes
+  
+  for(c2 in cluster){
+    n_cell = length(which(dom@clusters == c2))
+    if(n_cell > 1){
+      sig = rowMeans(dom@z_scores[genes, which(dom@clusters == c2)])
+    } else if(n_cell == 1){
+      sig = dom@z_scores[genes, which(dom@clusters == c2)]
+    } else {
+      sig = rep(0, length(genes))
+      names(sig) = genes
+    }
+    sig[which(sig < 0)] = 0
+    cl_expr[,c2] = sig
+  }
+  
+  if(length(dom@linkages$complexes) > 0){
+    cl_expr_coll_list <- avg_exp_for_complexes(cl_expr, complexes)
+    if(length(cl_expr_coll_list)>1){mat <- do.call(rbind, cl_expr_coll_list)}
+    cl_expr <- mat
+  }
+  
+  return(cl_expr)
+  
+} # avg_reclig_expr
+
+#' Create outgoing signaling network
+#' @description
+#' Create outgoing signaling network
+#' @param dom domino object
+#' @param outgoing_cluster vector of cluster names
+#' @param rec_cluster vector of cluster names
+#' @param plot_ligands vector of ligands to include
+#' @details Make a table
+#' @return table
+#' @export
+#' 
+outgoing_network <- function(dom, outgoing_cluster = NULL, rec_clusters = NULL, plot_ligands = NULL, verbose = T) {
+  
+  # Get ligands and receptors
+  lig_ls <- get_all_reclig(dom, rec_lig = "rec", expressed_only = T)
+  rec_ls <- get_all_reclig(dom, rec_lig = "lig", expressed_only = T)
+  
+  if (is.null(outgoing_cluster)) outgoing_cluster <- levels(dom@clusters)
+  if (is.null(rec_clusters)) rec_clusters <- levels(dom@clusters)
+  
+  #Create a matrix of all ligands expressed by the outgoing clusters 
+  # values are mean expression of all cells in a cluster
+  cl_ligands <- avg_reclig_expr(dom, cluster = outgoing_cluster, genes = lig_ls$genes, complexes = lig_ls$complex)
+  cl_receptors <- avg_reclig_expr(dom, cluster = rec_clusters, genes = rec_ls$genes, complexes = rec_ls$complex)
+  
+  # Melt and subset
+  cl_ligands_sub <- reshape2::melt(cl_ligands)
+  colnames(cl_ligands_sub) <- c("ligand", "cluster", "mean_z_score")
+  cl_ligands_sub <- cl_ligands_sub[cl_ligands_sub$mean_z_score > 0, ]
+  
+  cl_receptors_sub <- reshape2::melt(cl_receptors)
+  colnames(cl_receptors_sub) <- c("receptor", "cluster", "mean_z_score")
+  cl_receptors_sub <- cl_receptors_sub[cl_receptors_sub$mean_z_score > 0, ]
+  
+  #If a list of interesting genes provided - ignore for now
+  if(! is.null(plot_ligands)) {
+    cl_ligands_sub <- cl_ligands_sub[cl_ligands_sub$ligand %in% plot_ligands, ]
+  }
+  
+  if(is.null(rec_clusters)) {
+    rec_clusters <- levels(dom@clusters)
+  }
+  
+  
+  #Create a df of this type for the shortlisted ligands and specified recieving clusters.
+  #  ligand transcription.factor receptor ligand_exp receptor_exp sending           recv
+  #1 CXCL12            POU2F2...    CXCR4  0.4049881    1.1570554 Pericytes Mono and Mac
+  #2  ITGB1            POU2F2...   TYROBP  0.7428652    0.5838909 Pericytes Mono and Mac
+  #3 SEMA6D            POU2F2...   TYROBP  0.2379150    0.5838909 Pericytes Mono and Mac
+  #4 CXCL12            POU2F2...   LILRB2  0.4049881            0 Pericytes Mono and Mac
+  #5   CSF1            POU2F2...    CSF1R  0.2704230   1.66534274 Pericytes Mono and Mac
+  #6   IL34            POU2F2...    CSF1R  0.4082755   1.66534274 Pericytes Mono and Mac
+  
+  df <- data.frame()
+  
+  ### Can't figure out how to add the receptor expression.....
+  ### May be that this is entirely incorrect and I have to do this twice, once for 
+  
+  for(cl in rec_clusters) {
+    if (verbose) cat(sprintf("Building network for: %s\n", cl))
+    for(tf in names(dom@linkages$clust_tf_rec[[cl]])) {
+      recs = dom@linkages$clust_tf_rec[[cl]][[tf]]
+      if(length(recs)){
+        for(rec in recs) {
+          ligs <- dom@linkages$rec_lig[[rec]]
+          ligs <- resolve_names(dom, ligs)
+          if(length(ligs)) {
+            df.temp <- cl_ligands_sub[cl_ligands_sub$ligand %in% ligs, ]
+            df.temp2 <- cl_receptors_sub[cl_receptors_sub$receptor %in% rec,]
+            if(nrow(df.temp)) { #The df is non-empty
+              if (nrow(df.temp2) > 0) { rec_exp <- df.temp2$mean_z_score } else { rec_exp <- 0}
+              new.row <- data.frame(ligand = df.temp$ligand, 
+                                    transcription.factor = tf, 
+                                    receptor = rec, 
+                                    ligand_exp = df.temp$mean_z_score, 
+                                    receptor_exp = rep(rec_exp, nrow(df.temp)),
+                                    sending = df.temp$cluster,
+                                    recv = cl)
+              df <- rbind(df, new.row)  
+            } # fi nrow(df.temp)
+          } # fi length(ligs)
+        } # for rec
+      } # if length(recs)
+    } # for tf
+  } # for cl
+  return(df)
+} # outgoingNetwork
+
diff --git a/R/utils.R b/R/utils.R
index 65fd7639..fb744fb9 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -114,7 +114,7 @@ dom_correlations <- function(dom) {
 #' tf_rec_by_cluster <- dom_linkages(domino2:::pbmc_dom_built_tiny, "tf-receptor", TRUE)
 #' 
 dom_linkages <- function(dom, link_type = c(
-                            "complexes", "receptor-ligand",
+                            "complexes", "receptor-ligand", "ligand-receptor",
                             "tf-target", "tf-receptor", "receptor", "incoming-ligand"
                         ), by_cluster = FALSE) {
     links <- slot(dom, "linkages")
@@ -137,6 +137,8 @@ dom_linkages <- function(dom, link_type = c(
             return(links$tf_targets)
         } else if (link_type == "tf-receptor") {
             return(links$tf_rec)
+	} else if (link_type == "ligand-receptor") {
+	    return(links$lig_rec)
         } else {
             stop("This linkage type is not available.")
         }
diff --git a/man/avg_exp_for_complexes.Rd b/man/avg_exp_for_complexes.Rd
new file mode 100644
index 00000000..32577004
--- /dev/null
+++ b/man/avg_exp_for_complexes.Rd
@@ -0,0 +1,24 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/processing_fxns.R
+\name{avg_exp_for_complexes}
+\alias{avg_exp_for_complexes}
+\title{Average Complex Expression}
+\usage{
+avg_exp_for_complexes(exp_mat, complexes_list)
+}
+\arguments{
+\item{exp_mat}{expression matrix of ligands that may or may not be in a complex}
+
+\item{complexes_list}{named list. names = ligands, elements = ligands and/or complex compononents}
+}
+\value{
+list
+}
+\description{
+Average complex expression
+}
+\details{
+Search each element of complexes_list and determine if actually part of a complex.
+If it is, calculate colMeans for all ligands in the complex. If not, return the ligand expression.
+rewrite with purrr?
+}
diff --git a/man/avg_reclig_expr.Rd b/man/avg_reclig_expr.Rd
new file mode 100644
index 00000000..13e12833
--- /dev/null
+++ b/man/avg_reclig_expr.Rd
@@ -0,0 +1,26 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/processing_fxns.R
+\name{avg_reclig_expr}
+\alias{avg_reclig_expr}
+\title{Average Receptor/Ligand Expression}
+\usage{
+avg_reclig_expr(dom, cluster = NULL, genes, complexes)
+}
+\arguments{
+\item{dom}{domino object}
+
+\item{cluster}{vector of clusters to include (NULL uses all)}
+
+\item{genes}{genes to get averages of}
+
+\item{complexes}{list of resolved complexes}
+}
+\value{
+foo
+}
+\description{
+Create a matrix of average expression of cells in provided clsuters
+}
+\details{
+foo
+}
diff --git a/man/circos_ligand_receptor.Rd b/man/circos_ligand_receptor.Rd
index 0cad9c2e..d794c7a4 100644
--- a/man/circos_ligand_receptor.Rd
+++ b/man/circos_ligand_receptor.Rd
@@ -9,7 +9,8 @@ circos_ligand_receptor(
   receptor,
   ligand_expression_threshold = 0.01,
   cell_idents = NULL,
-  cell_colors = NULL
+  cell_colors = NULL,
+  title = paste0(receptor, " Signaling")
 )
 }
 \arguments{
diff --git a/man/circos_ligand_receptor_general.Rd b/man/circos_ligand_receptor_general.Rd
new file mode 100644
index 00000000..45176ed2
--- /dev/null
+++ b/man/circos_ligand_receptor_general.Rd
@@ -0,0 +1,47 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plot_fxns.R
+\name{circos_ligand_receptor_general}
+\alias{circos_ligand_receptor_general}
+\title{Plot expression of a receptor's ligands by other cell types as a chord plot}
+\usage{
+circos_ligand_receptor_general(
+  dom,
+  receptor = NULL,
+  ligand = NULL,
+  which_v = NULL,
+  expression_threshold = 0.01,
+  cell_idents = NULL,
+  cell_colors = NULL,
+  title = "Signaling",
+  label_arcs = T,
+  multi_plot = NULL
+)
+}
+\arguments{
+\item{dom}{Domino object that has undergone network building with \code{\link[=build_domino]{build_domino()}}}
+
+\item{receptor}{Name of a receptor active in at least one cell type in the domino object}
+
+\item{ligand}{Name of a ligand active in at least one cell type in the domino object}
+
+\item{which_v}{If both receptor and ligand are set, this must be set to either "receptor" or "ligand" to determine which to use.}
+
+\item{expression_threshold}{Minimum mean expression value of a ligand by a cell type for a chord to be rendered between the cell type and the receptor}
+
+\item{cell_idents}{Vector of cell types from cluster assignments in the domino object to be included in the plot.}
+
+\item{cell_colors}{Named vector of color names or hex codes where names correspond to the plotted cell types and the color values}
+
+\item{title}{if provided, final title will be paste0(receptor/ligand, " Signaling")}
+
+\item{label_arcs}{logical indicating if arc populations should be labeled}
+
+\item{multi_plot}{numeric vector that, if present indicates which plot in a list of plots this one is (for legend output)}
+}
+\value{
+renders a circos plot to the active graphics device
+}
+\description{
+Creates a chord plot of expression of ligands that can activate a specified
+receptor where chord widths correspond to mean ligand expression by the cluster.
+}
diff --git a/man/dom_linkages.Rd b/man/dom_linkages.Rd
index 42238c50..cc170697 100644
--- a/man/dom_linkages.Rd
+++ b/man/dom_linkages.Rd
@@ -6,8 +6,8 @@
 \usage{
 dom_linkages(
   dom,
-  link_type = c("complexes", "receptor-ligand", "tf-target", "tf-receptor", "receptor",
-    "incoming-ligand"),
+  link_type = c("complexes", "receptor-ligand", "ligand-receptor", "tf-target",
+    "tf-receptor", "receptor", "incoming-ligand"),
   by_cluster = FALSE
 )
 }
diff --git a/man/feat_heatmap.Rd b/man/feat_heatmap.Rd
index 11cf3ae2..18ef4c1a 100644
--- a/man/feat_heatmap.Rd
+++ b/man/feat_heatmap.Rd
@@ -7,6 +7,7 @@
 feat_heatmap(
   dom,
   feats = NULL,
+  clust = NULL,
   bool = FALSE,
   bool_thresh = 0.2,
   title = TRUE,
@@ -23,6 +24,8 @@ feat_heatmap(
 
 \item{feats}{Either a vector of features to include in the heatmap or 'all' for all features. If left NULL then the features selected for the signaling network will be shown.}
 
+\item{clust}{A vector of clusters in the domino object that should be included in the output. NULL (default) will use all clusters.}
+
 \item{bool}{Boolean indicating whether the heatmap should be continuous or boolean. If boolean then bool_thresh will be used to determine how to define activity as positive or negative.}
 
 \item{bool_thresh}{Numeric indicating the threshold separating 'on' or 'off' for feature activity if making a boolean heatmap.}
diff --git a/man/gene_network.Rd b/man/gene_network.Rd
index b252f1cd..ba13ae7a 100644
--- a/man/gene_network.Rd
+++ b/man/gene_network.Rd
@@ -12,6 +12,8 @@ gene_network(
   cols = NULL,
   lig_scale = 1,
   layout = "grid",
+  out_name = NULL,
+  in_name = NULL,
   ...
 )
 }
diff --git a/man/get_all_reclig.Rd b/man/get_all_reclig.Rd
new file mode 100644
index 00000000..1f8f9232
--- /dev/null
+++ b/man/get_all_reclig.Rd
@@ -0,0 +1,22 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/convenience_fxns.R
+\name{get_all_reclig}
+\alias{get_all_reclig}
+\title{Get All Ligands}
+\usage{
+get_all_reclig(dom, expressed_only = T, rec_lig = "rec")
+}
+\arguments{
+\item{dom}{domino object}
+
+\item{expressed_only}{logical indicating whether to subset ligands based on expression in dom@z_scores}
+}
+\value{
+vector of ligands if expressed_only = T, list if F
+}
+\description{
+Get all unique ligands present in dom@linkages$rec_lig
+}
+\details{
+Get all unique ligands in a domino object, expanding all ligand complexes as well. Optionally subset by expression.
+}
diff --git a/man/incoming_signaling_heatmap.Rd b/man/incoming_signaling_heatmap.Rd
index 17102364..20b6f6c4 100644
--- a/man/incoming_signaling_heatmap.Rd
+++ b/man/incoming_signaling_heatmap.Rd
@@ -10,6 +10,8 @@ incoming_signaling_heatmap(
   clusts = NULL,
   min_thresh = -Inf,
   max_thresh = Inf,
+  display_top = NULL,
+  display_method = "mean",
   scale = "none",
   normalize = "none",
   title = TRUE,
@@ -27,6 +29,10 @@ incoming_signaling_heatmap(
 
 \item{max_thresh}{Maximum signaling threshold for plotting. Defaults to Inf for no threshold.}
 
+\item{display_top}{Number of top receptor-ligand pairs to display in output (rather than all results)}
+
+\item{display_method}{How to determine top pairs. Options are 'mean' or 'median' of pair's expression across all incoming populations, or 'max' to take the top N invidivual expression values with no repeats}
+
 \item{scale}{How to scale the values (after thresholding). Options are 'none', 'sqrt' for square root, or 'log' for log10.}
 
 \item{normalize}{Options to normalize the matrix. Accepted inputs are 'none' for no normalization, 'rec_norm' to normalize to the maximum value with each receptor cluster, or 'lig_norm' to normalize to the maximum value within each ligand cluster}
diff --git a/man/invert_rec_lig_expr.Rd b/man/invert_rec_lig_expr.Rd
new file mode 100644
index 00000000..929a5918
--- /dev/null
+++ b/man/invert_rec_lig_expr.Rd
@@ -0,0 +1,18 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/processing_fxns.R
+\name{invert_rec_lig_expr}
+\alias{invert_rec_lig_expr}
+\title{Invert Receptor Ligand Data}
+\usage{
+invert_rec_lig_expr(dom)
+}
+\arguments{
+\item{dom}{domino object built by domino build}
+}
+\value{
+domino object with updated \code{rec_signaling} slot and \code{cr_signaling_matrices} slot
+}
+\description{
+Reformat exisiting domino object to also record expression data/receptor mapping
+from the ligand-receptor direction rather than receptor-ligand
+}
diff --git a/man/invert_rec_lig_linkages.Rd b/man/invert_rec_lig_linkages.Rd
new file mode 100644
index 00000000..5155b75d
--- /dev/null
+++ b/man/invert_rec_lig_linkages.Rd
@@ -0,0 +1,18 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/processing_fxns.R
+\name{invert_rec_lig_linkages}
+\alias{invert_rec_lig_linkages}
+\title{Invert Receptor Ligand Linkages}
+\usage{
+invert_rec_lig_linkages(dom)
+}
+\arguments{
+\item{dom}{domino object built by domino build}
+}
+\value{
+domino object with new entry in dom@linkages called \code{lig_rec}, which is a list
+with one element per ligand and the list values are all receptors mapped to that ligand by rl_map
+}
+\description{
+Use dom@misc$rl_map to create inverted version of dom@linkages$rec_lig
+}
diff --git a/man/outgoing_network.Rd b/man/outgoing_network.Rd
new file mode 100644
index 00000000..c80e9725
--- /dev/null
+++ b/man/outgoing_network.Rd
@@ -0,0 +1,32 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/processing_fxns.R
+\name{outgoing_network}
+\alias{outgoing_network}
+\title{Create outgoing signaling network}
+\usage{
+outgoing_network(
+  dom,
+  outgoing_cluster = NULL,
+  rec_clusters = NULL,
+  plot_ligands = NULL,
+  verbose = T
+)
+}
+\arguments{
+\item{dom}{domino object}
+
+\item{outgoing_cluster}{vector of cluster names}
+
+\item{plot_ligands}{vector of ligands to include}
+
+\item{rec_cluster}{vector of cluster names}
+}
+\value{
+table
+}
+\description{
+Create outgoing signaling network
+}
+\details{
+Make a table
+}
diff --git a/man/resolve_complexes.Rd b/man/resolve_complexes.Rd
new file mode 100644
index 00000000..279d97ac
--- /dev/null
+++ b/man/resolve_complexes.Rd
@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/convenience_fxns.R
+\name{resolve_complexes}
+\alias{resolve_complexes}
+\title{Resolve Complexes}
+\usage{
+resolve_complexes(dom, genes)
+}
+\arguments{
+\item{dom}{domino object}
+
+\item{genes}{vector of complex names to resolve}
+}
+\value{
+list of length(genes), with names(list) == genes. List values are the same as the names if not in a complex and are the complex genes if they are in a complex.
+}
+\description{
+Expand any complex names into their component gene names
+}
+\details{
+Ligand names (which are stored in dom@linkages$rec_lig) can refer to complexes that are
+detailed in dom@linkages$complexes. Search all provided names and if the name is a complex, output component genes, otherwise output the name
+}
diff --git a/man/resolve_names.Rd b/man/resolve_names.Rd
new file mode 100644
index 00000000..37bdb97a
--- /dev/null
+++ b/man/resolve_names.Rd
@@ -0,0 +1,23 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/convenience_fxns.R
+\name{resolve_names}
+\alias{resolve_names}
+\title{Resolve Names}
+\usage{
+resolve_names(dom, genes, rec_lig = "lig")
+}
+\arguments{
+\item{dom}{domino object}
+
+\item{genes}{vector of gene names to resolve}
+}
+\value{
+vector of length(genes) with applicable values replaced using dom@misc$rl_map information
+}
+\description{
+Return gene names of ligands with non-standard names
+}
+\details{
+Ligand names (which are stored in dom@linkages$rec_lig) do not always match the gene name
+Search all provided names and output the gene name if the ligand name is non-standard
+}
diff --git a/man/signaling_network.Rd b/man/signaling_network.Rd
index 03db3a48..111129b0 100644
--- a/man/signaling_network.Rd
+++ b/man/signaling_network.Rd
@@ -19,6 +19,7 @@ signaling_network(
   scale_by = "rec_sig",
   vert_scale = 3,
   plot_title = NULL,
+  offset_sides = NULL,
   ...
 )
 }
@@ -51,6 +52,8 @@ signaling_network(
 
 \item{plot_title}{Text for the plot's title.}
 
+\item{offset_sides}{optional parameter to offset the side labels wihen layout is type 'circle'. Options are 'above' or 'below'}
+
 \item{...}{Other parameters to be passed to plot when used with an \code{{igraph}} object.}
 }
 \value{