decomposing community done

snaketron · snaketron · commit 6d8fbbfa3814 · 2025-04-07T14:24:01.000+02:00
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: ClustIRR
 Type: Package
 Title: Clustering of immune receptor repertoires
-Version: 1.5.48
+Version: 1.5.50
 Authors@R: c(
     person("Simo", "Kitanovski", email = "simokitanovski@gmail.com", 
     role = c("aut", "cre"), comment=c(ORCID="0000-0003-2909-5376")),
diff --git a/R/community.R b/R/community.R
@@ -438,7 +438,6 @@ decode_communities <- function(community_id,
                 V(sg)$components <- components(graph = sg)$membership
                 V(sg)$component_id <- paste0(V(sg)$key, '|', 
                                              V(sg)$components)
-                V(sg)$component_id <- as.numeric(as.factor(V(sg)$component_id))
                 return(disjoint_union(lapply(
                     X = unique(V(sg)$component_id), g = sg, 
                     FUN = function(x, g) {
@@ -447,7 +446,15 @@ decode_communities <- function(community_id,
                     })))
             })
         sgs <- disjoint_union(sgs)
-        return(sgs)
+        V(sgs)$component_id <- as.numeric(as.factor(V(sgs)$component_id))
+        
+        sgs_stat <- get_component_stats(x = sgs)
+        
+        vs <- as_data_frame(x = sgs, what = "vertices")
+        
+        return(list(community_graph = sgs, 
+                    component_stats = sgs_stat,
+                    node_summary = vs))
     } 
     else {
         sgs <- lapply(
@@ -469,8 +476,67 @@ decode_communities <- function(community_id,
                     })))
             })
         sgs <- disjoint_union(sgs)
-        return(sgs)
         
+        sgs_stat <- get_component_stats(x = sgs)
+        
+        vs <- as_data_frame(x = sgs, what = "vertices")
+        
+        return(list(community_graph = sgs, 
+                    component_stats = sgs_stat,
+                    node_summary = vs))
     }
 }
 
+
+get_component_stats <- function(x) {
+    
+    # what is he
+    community_id <- V(x)$community[1]
+    
+    # number of nodes per component
+    count_n <- table(V(x)$component_id)
+    
+    # components has more than one node
+    i <- which(count_n>1)
+    stats_components <- c()
+    if(length(i)>0) {
+        stats_components <- lapply(X = names(i), g = x, FUN = function(x, g) {
+            o <- subgraph(g, vids = which(V(g)$component_id==x))
+            
+            ncweight <- E(o)$ncweight
+            if(length(ncweight)==0) {
+                ncweight <- NA
+            }
+            nweight <- E(o)$nweight
+            if(length(nweight)==0) {
+                nweight <- NA
+            }
+            
+            return(data.frame(component_id = x,
+                              community = V(o)$community[1],
+                              mean_ncweight = mean(ncweight),
+                              mean_nweight = mean(nweight),
+                              n_nodes = length(o),
+                              n_edges = length(E(o)),
+                              n_clique_edges = length(o)*(length(o)-1)/2,
+                              diameter = diameter(o)))
+        })
+        stats_components <- do.call(rbind, stats_components)
+    }
+    
+    i <- which(count_n==1)
+    stats_singletons <- c()
+    if(length(i)>0) {
+        stats_singletons <- data.frame(component_id = names(i),
+                                       community = community_id,
+                                       mean_ncweight = NA,
+                                       mean_nweight = NA,
+                                       n_nodes = 1,
+                                       n_edges = 0,
+                                       n_clique_edges = 0,
+                                       diameter = 0)
+    }
+    
+    return(rbind(stats_components, stats_singletons))
+}
+
diff --git a/man/decode_communities.Rd b/man/decode_communities.Rd
@@ -42,7 +42,12 @@ have the same attribute values among \bold{ALL} provided attributes will
 be treated as a subcomponent.}
 }
 \value{
-The output is a "filtered" igraph object.
+The output is a list with the following elements
+\itemite{
+\item community_graph, "filtered" igraph object
+\item component_stats, data.frame with summary about each connected component
+\item node_summary, data.frame with summary about each node
+}
 }
 \examples{
 # load package input data
@@ -76,9 +81,9 @@ library(igraph)
 edge_attr_names(graph = gcd$graph)
 
 # For instance, the following edge-filter will instruct ClustIRR to keep 
-# edges with: edge attributes: nweight>=3 \bold{AND} ncweight>=3
-edge_filter <- rbind(data.frame(name = "nweight", value = 3, operation = ">="),
-                     data.frame(name = "ncweight", value = 3, operation = ">="))
+# edges with: edge attributes: nweight>=8 \bold{AND} ncweight>=8
+edge_filter <- rbind(data.frame(name = "nweight", value = 8, operation = ">="),
+                     data.frame(name = "ncweight", value = 8, operation = ">="))
 
 # In addition, we can construct filters based on the following node attributes:
 vertex_attr_names(graph = gcd$graph)
@@ -96,10 +101,14 @@ c1 <- decode_communities(community_id = 1,
 
 # Plot resulting igraph
 par(mar = c(0, 0, 0, 0))
-plot(c1, vertex.size = 10)
+plot(c1$community_graph, vertex.size = 10)
 
 # Now look at node attributes
-as_data_frame(x = c1, what = "vertices")[,c("name", "component_id", 
-                                            "CDR3b", "CDR3a", "Ag_gene")]
+as_data_frame(x = c1$community_graph, what = "vertices")
+    [,c("name", "component_id", "CDR3b", "CDR3a", "Ag_gene")]
+
+str(c1$component_stats)
+
+str(c1$node_summary)
 }
 
diff --git a/vignettes/User_manual.Rmd b/vignettes/User_manual.Rmd
@@ -201,9 +201,9 @@ cl <- clustirr(s = tcr_reps, meta = meta, control = list(gmi = 0.7))
 
 The output **complex** list with three elements:
 
-1.   `graph` = the joint graph $J$
+1.  `graph` = the joint graph $J$
 
-2.   `clust_irrs` = list with S4 `clust_irr` objects one for each TCR
+2.  `clust_irrs` = list with S4 `clust_irr` objects one for each TCR
     repertoire
 
 3.  `multigraph` = logical value which tell us whether $J$ contains
@@ -389,7 +389,8 @@ wrap_plots(honeycomb, nrow = 1)+
 ```
 
 Also see `community_summary`. In the data.frame `wide` we provide
-community summaries in each community across all TCR repertoires, including:
+community summaries in each community across all TCR repertoires,
+including:
 
 -   `clones_a`, `clone_b`, `clone_c`, `clones_n`: the frequency of
     clones in the community coming from repertoire `a`, `b`, `c` and in
@@ -406,8 +407,8 @@ community summaries in each community across all TCR repertoires, including:
 kable(head(gcd$community_summary$wide), digits = 1, row.names = FALSE)
 ```
 
-In the data.frame `tall` we provide community and repertoire
-summaries in each row.
+In the data.frame `tall` we provide community and repertoire summaries
+in each row.
 
 ```{r}
 kable(head(gcd$community_summary$tall), digits = 1, row.names = FALSE)
@@ -448,11 +449,11 @@ edge_attr_names(graph = gcd$graph)
 ```
 
 The following edge-filter will instruct ClustIRR to keep edges with with
-edge attributes: nweight $>=$ 4 **AND** ncweight $>=$ 4
+edge attributes: nweight $>=$ 8 **AND** ncweight $>=$ 8
 
 ```{r}
-edge_filter <- rbind(data.frame(name = "nweight", value = 4, operation = ">="),
-                     data.frame(name = "ncweight", value = 4, operation = ">="))
+edge_filter <- rbind(data.frame(name = "nweight", value = 8, operation = ">="),
+                     data.frame(name = "ncweight", value = 8, operation = ">="))
 ```
 
 ```{r}
@@ -477,17 +478,30 @@ c8 <- decode_communities(community_id = 8,
 
 ```{r, fig.width=6, fig.height=4}
 par(mar = c(0, 0, 0, 0))
-plot(c8, vertex.size = 10)
+plot(c8$community, vertex.size = 10)
 ```
 
 ```{r}
-kable(as_data_frame(x = c8, what = "vertices")[, c("name", "key", 
-                                                   "CDR3b", "TRBV", 
-                                                   "TRBJ", "CDR3a", 
-                                                   "TRAV", "TRAJ")],
+kable(as_data_frame(x = c8$community, what = "vertices")
+      [, c("name", "component_id", "CDR3b", "TRBV", 
+           "TRBJ", "CDR3a", "TRAV", "TRAJ")],
       row.names = FALSE)
 ```
 
+... or we can summarize the properties of each component in the next
+table as rows with:
+
+-   component id, community id (=8 in this example because this is what
+    we selected)
+-   mean edge weights (for all core and complete CDR3 pairs)
+-   number of nodes, number of edges and expected number edges if the
+    component were a clique
+-   diameter of the component
+
+```{r}
+kable(c8$component_stats, row.names = FALSE)
+```
+
 ## **Step 5.** differential community occupancy (DCO) with `dco`
 
 Do we see **expanding** or **contracting** communities in a given
@@ -604,7 +618,6 @@ community occupancy (DCO):
 \delta^{a-b}_i = \beta^a_i - \beta^b_i
 \end{align}
 ```
-
 Importantly, `r Biocpkg("ClustIRR")` always computes both contrasts ($a$
 vs. $b$ and $b$ vs. $a$): $\delta^{a-b}_i$ and $\delta^{b-a}_i$.
 
@@ -629,7 +642,6 @@ quantifies absolute differences in community probabilities:
 \mathrm{softmax}(\alpha + \beta^b)
 \end{align}
 ```
-
 Again, both contrasts are computed: $\epsilon^{a-b}_i$ and
 $\epsilon^{b-a}_i$.
 
