@@ -2,7 +2,7 @@ Type: Package

 Package: CellaRepertorium

 Title: Data structures, clustering and testing for single 

     cell immune receptor repertoires (scRNAseq RepSeq/AIRR-seq)

-Version: 1.1.2

+Version: 1.1.3

 Authors@R: 

     c(person(given = "Andrew",

              family = "McDavid",

@@ -60,7 +60,9 @@ Suggests:

     SingleCellExperiment,

     scater,

     broom.mixed,

-    cowplot

+    cowplot,

+    igraph,

+    ggraph

 LinkingTo: 

     Rcpp

 VignetteBuilder: 

@@ -31,6 +31,7 @@ export(modal_category)

 export(mutate_cdb)

 export(np)

 export(pairing_tables)

+export(plot_cluster_factors)

 export(purity)

 export(rank_chain_ccdb)

 export(rank_prevalence_ccdb)

@@ -1,67 +1,115 @@

-#' Visualization of clustering characteristics for a ContigCellDB object with inherited fine_clustering object

-#' 

+require_igraph = function(){

+  if(!requireNamespace('igraph') || !requireNamespace('ggraph')) stop('Install `igraph` and `ggraph`.')

+  TRUE

+}

+

+#' Visualization of pairs of cluster factor

+#'

+#' With `factors`, a pair of variables present in the `contig_tbl` and the `cluster_tbl`,

+#' generate and plot cross-tabs of the number of contigs, or its pearson residual.

 #' @param ccdb A ContigCellDB object.

-#' @param cluster_factors A character vector of variables in ccdb$contig table to visualize clusters by

+#' @param factors `character` length 2 of fields present

 #' @param type Type of visualization, a heatmap or a node-edge network plot

 #' @param statistic Cluster characteristics visualized by pearson residuals or raw contig counts

-#' @param ncluster Optional integer. Only plot characteristics that have at least nclusters (for clarity of visualization)

+#' @param ncluster `integer`.  Omit factors that occur less than `nclusters`. For clarity of visualization.

 #' @param chaintype Character in ccdb$contig_tbl$chain. If passed will subset contigs belonging to specified chain (IGH,IGK,IGL,TRA,TRB)

-#'

+#' @seealso canonicalize_cluster to "roll-up" additional contig variables into the `cluster_tbl``

 #' @return A ggraph object if type == 'network', and a ggplot object if type == 'heatmap'

+#' @export

 #' @examples

+#' library(ggraph)

 #' data(ccdb_ex)

-#' ccdb_germline_ex = cluster_germline(ccdb_ex, segment_keys = c('v_gene', 'j_gene', 'chain'), cluster_pk = 'segment_idx')

-#' ccdb_germline_ex = fine_clustering(germline_cluster, sequence_key = 'cdr3_nt', type = 'DNA')

-#' plot_clustering(ccdb_germline_ex,cluster_factors = c('v_gene','j_gene'), statistic = 'pearson', type = 'network' ,ncluster = 50, chaintype = 'IGH')

+#' ccdb_germline_ex = cluster_germline(ccdb_ex, segment_keys = c('v_gene', 'j_gene', 'chain'),

+#' cluster_pk = 'segment_idx')

+#' ccdb_germline_ex = fine_clustering(ccdb_germline_ex, sequence_key = 'cdr3_nt', type = 'DNA')

+#' plot_cluster_factors(ccdb_germline_ex,factors = c('v_gene','j_gene'),

+#' statistic = 'pearson', type = 'network' ,ncluster = 10, chaintype = 'TRB')

+#' plot_cluster_factors(ccdb_germline_ex,factors = c('v_gene','j_gene'),

+#' statistic = 'contigs', type = 'heatmap')

+#' plot_cluster_factors(ccdb_germline_ex,factors = c('v_gene','j_gene'),

+#' statistic = 'contigs', type = 'network', ncluster = 10)

+plot_cluster_factors = function(ccdb,

+                           factors,

+                           type = c('heatmap', 'network'),

+                           statistic = c('pearson', 'contigs'),

+                           ncluster = 0,

+                           chaintype) {

+  if (!inherits(ccdb,  "ContigCellDB") || !has_fineclustering(ccdb))

+    stop('ccdb must be  `ContigCellDB` with `fine_clustering`.')

+  check_plot_infra()

+  statistic = match.arg(statistic, c('pearson', 'contigs'))

+  type = match.arg(type, c('heatmap', 'network'))

+  if (length(factors) != 2 ||

+      !is.character(factors))

+    stop('`factors` must be character vector that specifies exactly two columns from ccdb')

-  

-plot_clustering = function(ccdb,cluster_factors,type = c('heatmap','network'),statistic = c('pearson','contigs'),ncluster,chaintype) {

+  cluster_tbl =  dplyr::filter(ccdb$cluster_tbl, n_cluster > ncluster)

-if(!inherits(ccdb,  "ContigCellDB")) stop('ccdb must have class ContigCellDb')

-if(length(cluster_factors) != 2 | !is.character(cluster_factors)) stop('cluster_factors must be character vector that specifies exactly two columns from ccdb')

-  

-if(!missing(ncluster)) {tmp = ccdb$cluster_tbl %>% filter(n_cluster > ncluster)} else {tmp = ccdb$cluster_tbl}

+  if (!missing(chaintype)) {

+    cluster_tbl = cluster_tbl %>% filter(chain %in% chaintype)

+  }

-if(!missing(chaintype)) {tmp = tmp %>% filter(chain == chaintype)} else {tmp = tmp}

-  

-if(statistic == 'pearson') {

+  if(!nrow(cluster_tbl) > 0) stop("No clusters remain after filtering.")

-  if(!missing(chaintype)) {ccdb$contig_tbl = ccdb$contig_tbl %>% filter(chain == chaintype)}

-  

-  resids = chisq.test(table(ccdb$contig_tbl %>% pull(cluster_factors[1]), ccdb$contig_tbl %>% pull(cluster_factors[2])))

-  

-  resid_tab = as.data.frame(resids$residuals)

-  

-  colnames(resid_tab) = c(cluster_factors,'residuals')

-  

-  tmp = left_join(tmp,resid_tab,by = cluster_factors)

-  

-  fillby = 'residuals'

-  

-}

-  

-else {fillby = 'n_cluster'}

-  

-tmp = tmp %>% relocate(cluster_factors)

+  if (statistic == 'pearson') {

+    if (!missing(chaintype)) {

+      ccdb$contig_tbl = ccdb$contig_tbl %>% filter(chain %in% chaintype)

+    }

+    if(!nrow(ccdb$contig_tbl)>0) stop("No contigs left after filtering.")

+    resids = stats::chisq.test(table(

+      ccdb$contig_tbl[[factors[1]]],

+      ccdb$contig_tbl[[factors[2]]]

+    ))

+    resid_tab = as.data.frame(resids$residuals)

+    colnames(resid_tab) = c(factors, 'residuals')

+    cluster_tbl = left_join_warn(cluster_tbl, resid_tab, by = factors)

+    fillby = 'residuals'

+  } else {

+    fillby = 'n_cluster'

+  }

+  cluster_tbl = cluster_tbl %>% dplyr::relocate(factors)

-if (type == 'network') {

-  

-  nodes = data.frame(names = c(unique(tmp %>% pull(cluster_factors[1])),unique(tmp %>% pull(cluster_factors[2]))))

-  

-  graph = graph_from_data_frame(tmp,vertices = nodes)

+  if (type == 'network') {

+    require_igraph()

+    nodes = data.frame(names = c(unique(cluster_tbl[[factors[1]]]),unique(cluster_tbl[[factors[2]]])))

+    graph = igraph::graph_from_data_frame(cluster_tbl, vertices = nodes)

-  plot = ggraph(graph,layout = 'circle') + geom_node_point() + geom_edge_fan(aes_string(colour = fillby)) + geom_node_label(aes(label = V(graph)$name), size = 2) + ggtitle(paste('Network Plot for factors', cluster_factors[1],'and',cluster_factors[2], 'using',statistic))

-  

-  if(statistic == 'pearson') {plot = plot + scale_edge_color_gradient2()}

-  

-  return(plot) }

-  

-else {

-  plot = ggplot(tmp) + aes_string(x = cluster_factors[1],y = cluster_factors[2],fill = fillby) + geom_tile() + theme(axis.text.x = element_text(angle = 90)) +  ggtitle(paste('Heatmap for factors', cluster_factors[1],'and',cluster_factors[2], 'using',statistic))

-  

-  if (statistic == 'pearson') {plot = plot + scale_fill_gradient2()}

+    plot = ggraph::ggraph(graph, layout = 'circle') + ggraph::geom_node_point() + ggraph::geom_edge_fan(ggplot2::aes_string(colour = fillby)) + ggraph::geom_node_label(ggplot2::aes(label = igraph::V(graph)$name), size = 2) + ggplot2::ggtitle(

+      paste(

+        'Network Plot for factors',

+        factors[1],

+        'and',

+        factors[2],

+        'using',

+        statistic

+      )

+    )

-return(plot) }

+    if (statistic == 'pearson') {

+      plot = plot + ggraph::scale_edge_color_gradient2()

+    }

+

+    return(plot)

+  }

+

+  else {

+    plot = ggplot2::ggplot(cluster_tbl) + ggplot2::aes_string(x = factors[1], y = factors[2], fill = fillby) + ggplot2::geom_tile() + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 90)) +  ggplot2::ggtitle(

+      paste(

+        'Heatmap for factors',

+        factors[1],

+        'and',

+        factors[2],

+        'using',

+        statistic

+      )

+    )

+

+    if (statistic == 'pearson') {

+      plot = plot + ggplot2::scale_fill_gradient2()

+    }

+

+    return(plot)

+  }

 }

-  

+

new file mode 100644

@@ -0,0 +1,50 @@

+% Generated by roxygen2: do not edit by hand

+% Please edit documentation in R/plot_clustering.R

+\name{plot_cluster_factors}

+\alias{plot_cluster_factors}

+\title{Visualization of pairs of cluster factor}

+\usage{

+plot_cluster_factors(

+  ccdb,

+  factors,

+  type = c("heatmap", "network"),

+  statistic = c("pearson", "contigs"),

+  ncluster = 0,

+  chaintype

+)

+}

+\arguments{

+\item{ccdb}{A ContigCellDB object.}

+

+\item{factors}{\code{character} length 2 of fields present}

+

+\item{type}{Type of visualization, a heatmap or a node-edge network plot}

+

+\item{statistic}{Cluster characteristics visualized by pearson residuals or raw contig counts}

+

+\item{ncluster}{\code{integer}.  Omit factors that occur less than \code{nclusters}. For clarity of visualization.}

+

+\item{chaintype}{Character in ccdb$contig_tbl$chain. If passed will subset contigs belonging to specified chain (IGH,IGK,IGL,TRA,TRB)}

+}

+\value{

+A ggraph object if type == 'network', and a ggplot object if type == 'heatmap'

+}

+\description{

+With \code{factors}, a pair of variables present in the \code{contig_tbl} and the \code{cluster_tbl},

+generate and plot cross-tabs of the number of contigs, or its pearson residual.

+}

+\examples{

+library(ggraph)

+data(ccdb_ex)

+ccdb_germline_ex = cluster_germline(ccdb_ex, segment_keys = c('v_gene', 'j_gene', 'chain'), cluster_pk = 'segment_idx')

+ccdb_germline_ex = fine_clustering(ccdb_germline_ex, sequence_key = 'cdr3_nt', type = 'DNA')

+plot_cluster_factors(ccdb_germline_ex,factors = c('v_gene','j_gene'),

+statistic = 'pearson', type = 'network' ,ncluster = 10, chaintype = 'TRB')

+plot_cluster_factors(ccdb_germline_ex,factors = c('v_gene','j_gene'),

+statistic = 'contigs', type = 'heatmap')

+plot_cluster_factors(ccdb_germline_ex,factors = c('v_gene','j_gene'),

+statistic = 'contigs', type = 'network', ncluster = 10)

+}

+\seealso{

+canonicalize_cluster to "roll-up" additional contig variables into the `cluster_tbl``

+}

deleted file mode 100644

@@ -1,40 +0,0 @@

-% Generated by roxygen2: do not edit by hand

-% Please edit documentation in R/plot_clustering.R

-\name{plot_clustering}

-\alias{plot_clustering}

-\title{Visualization of clustering characteristics for a ContigCellDB object with inherited fine_clustering object}

-\usage{

-plot_clustering(

-  ccdb,

-  cluster_factors,

-  type = c("heatmap", "network"),

-  statistic = c("pearson", "contigs"),

-  ncluster,

-  chaintype

-)

-}

-\arguments{

-\item{ccdb}{A ContigCellDB object.}

-

-\item{cluster_factors}{A character vector of variables in ccdb$contig table to visualize clusters by}

-

-\item{type}{Type of visualization, a heatmap or a node-edge network plot}

-

-\item{statistic}{Cluster characteristics visualized by pearson residuals or raw contig counts}

-

-\item{ncluster}{Optional integer. Only plot characteristics that have at least nclusters (for clarity of visualization)}

-

-\item{chaintype}{Character in ccdb$contig_tbl$chain. If passed will subset contigs belonging to specified chain (IGH,IGK,IGL,TRA,TRB)}

-}

-\value{

-A ggraph object if type == 'network', and a ggplot object if type == 'heatmap'

-}

-\description{

-Visualization of clustering characteristics for a ContigCellDB object with inherited fine_clustering object

-}

-\examples{

-data(ccdb_ex)

-ccdb_germline_ex = cluster_germline(ccdb_ex, segment_keys = c('v_gene', 'j_gene', 'chain'), cluster_pk = 'segment_idx')

-ccdb_germline_ex = fine_clustering(germline_cluster, sequence_key = 'cdr3_nt', type = 'DNA')

-plot_clustering(ccdb_germline_ex,cluster_factors = c('v_gene','j_gene'), statistic = 'pearson', type = 'network' ,ncluster = 50, chaintype = 'IGH')

-}

@@ -111,16 +111,16 @@ We can cluster by any other feature of the contigs. Here we cluster each contig

 ```{r}

 germline_cluster = fine_clustering(germline_cluster, sequence_key = 'cdr3_nt', type = 'DNA')

-ggplot(germline_cluster$cluster_tbl %>% filter(chain == 'TRB'), aes(x = v_gene, y = j_gene, fill = n_cluster)) + geom_tile() + theme(axis.text.x = element_text(angle = 90))

+filter_cdb(germline_cluster, chain == 'TRB') %>% plot_cluster_factors(factors = c('v_gene','j_gene'), statistic = 'contigs', type = 'heatmap')

 ```

-Number of pairs

+Number of pairs.  The pearson residual (showing the difference from expected counts given marginals) is probably more informative, set `statistic = 'residual'` for this.

 ```{r}

 ggplot(germline_cluster$cluster_tbl %>% filter(chain == 'TRB'), aes(x = v_gene, y = j_gene, fill = avg_distance)) + geom_tile() + theme(axis.text.x = element_text(angle = 90))

 ```

-Average Levenshtein distance of CDR3 within each pair.  This might be turned into a z-score by fitting a weighted linear model with sum-to-zero contrasts and returning the pearson residuals.  This could determine if a pairing has an unexpected small, or large, within cluster distance.

+Average Levenshtein distance of CDR3 within each pair.  This might be turned into a z-score by fitting a weighted linear model with sum-to-zero contrasts and returning the studentized residuals.  This could determine if a pairing has an unexpected small, or large, within cluster distance.