new file mode 100644

@@ -0,0 +1,3 @@

+^.*\.Rproj$

+^\.Rproj\.user$

+extdata/.*\.json$

@@ -1,70 +1,20 @@

-#' For each cell, return a single, canonical chain-cluster

-#'

-#' In single cell data, multiple chains (heavy-light or alpha-beta) are expected.  In some cases, there could be more than two (eg multiple alpha alleles for T cells).

-#' This picks a cluster id for each cell based on the overall prevalence of cluster ids over all cells in `tbl`.

-#' If order = 1 then the canonical chain-cluster will be the most prevalent, and if order = 2, it will be the 2nd most prevalent, and so on.  Ties are broken arbitrarily (possibly by lexicographic order of `cluster_idx`).

-#' @param tbl `data.frame` containing columns specified in `cell_identifiers` and `cluster_idx`

-#' @param cell_identifiers `character` vector specifying columns in `tbl` that identify a cell

-#' @param cluster_idx `character` specifying the column in `tbl` that identifies a

-#' @param order return the 1st, 2nd, 3rd, etc, most common chain-cluster

-#'

-#'

-#' @return `data.frame` with columns from `cell_identifiers` and a single `cluster_idx` for each cell

-#' @export

-get_canonical_chain = function(tbl, cell_identifiers = 'barcode', cluster_idx = 'cluster_idx', order = 1){

-    prevalence = tbl %>% group_by(!!sym(cluster_idx)) %>% summarize(prev = n())

-    tbl_prevalence = left_join(tbl, prevalence, by = cluster_idx)

-    tbl_order = tbl_prevalence %>% group_by(!!!syms(cell_identifiers)) %>% summarize(!!cluster_idx := dplyr::nth(!!sym(cluster_idx), -order, prev))

-    ungroup(tbl_order)

+cluster_germline = function(tbl, segment_identifiers = c('v_gene', 'j_gene', 'chain')){

+    seg_types = tbl %>% group_by(!!!syms(segment_identifiers)) %>% summarize() %>% ungroup() %>% mutate(germline_idx = seq_len(nrow(.)))

+    tbl %>% left_join(seg_types, by = segment_identifiers)

 }

-# Define canonical chain types per cell

-

-#' Given a family of similar sequences, return a "representative"

-#'

-#' @param seqs character vector

-#' @param medoid_idx optional index into seqs

-#' @param warn_if_distinct Should a warning be emitted if there are distinct elements in seqs?

-#'

-#' If medoid_idx is supplied, the medoid sequence is returned, otherwise the longest

-#' sequence is returned

-#'

-#' @return character vector

-#' @export

-#'

-#' @examples

-#' get_canonical_representative(c('apple', 'manzana', 'pomme'))

-get_canonical_representative = function(seqs, medoid_idx, warn_if_distinct = FALSE){

-    if(!missing(medoid_idx)){

-        if(!is.integer(medoid_idx) || medoid_idx < 1 || medoid_idx > length(seqs)) stop("Illegal `medoid_idx`")

-        ii = medoid_idx

-    } else{

-        len = str_length(seqs)

-        ii = which.max(len)

-    }

-    rep = seqs[ii]

-

-    if(warn_if_distinct && length(nomatch <- which(seqs != rep)) > 0) warning("At indices ", nomatch, " sequences did not match the representative")

-

-    return(rep)

-}

-

-

-

-

-

 fine_cluster_by = function(seqs, by, max_dist = NULL, ...){

     seq_list = split(seqs, by)

     fc_list = lapply(seq_list, fine_cluster, ...)

     if(is.null(max_dist)){

-        max_max = max(map_dbl(fc_list, 'max_dist'))

+        max_max = max(purrr::map_dbl(fc_list, 'max_dist'))

     } else {

         max_max = max_dist

     }

-    affinities = map(fc_list, ~ max_max - .[['distance']])

+    affinities = purrr::map(fc_list, ~ max_max - .[['distance']])

     aff_matrix = Matrix::.bdiag(affinities)

-    homologies = unlist(map(fc_list, 'homology'), use.names = FALSE)

-    medoid = map_dbl(fc_list, 'medoid')

+    homologies = unlist(purrr::map(fc_list, 'homology'), use.names = FALSE)

+    medoid = purrr::map_dbl(fc_list, 'medoid')

     list(affinities, aff_matrix, homologies, medoid)

 }

@@ -72,12 +22,15 @@ fine_cluster_by = function(seqs, by, max_dist = NULL, ...){

 #'

 #' The distances between AA sequences is defined to be 1-score/max(score) times the median length of the input sequences.

 #' The distances between nucleotide sequences is defined to be edit_distance/max(edit_distance) times the median length of input sequences.

+#'

 #' @param seqs character vector, DNAStringSet or AAStringSet

 #' @param type character either `AA` or `DNA` specifying type of `seqs`

 #' @param big_memory_brute attempt to cluster more than 4000 sequences?  Clustering is quadratic, so this will take a long time and might exhaust memory

 #' @param method one of 'substitutionMatrix' or 'levenshtein'

 #' @param substitution_matrix a character vector naming a substition matrix available in Biostrings, or a substitution matrix itself

+#' @param cluster_fun `character`, one of "hclust" or "none", determining if distance matrices should also be clustered with `hclust`

 #' @param cluster_method character passed to `hclust`

+#'

 #' @seealso hclust, stringDist

 #' @return dendrogram of class `hclust`

 #' @export

@@ -87,10 +40,10 @@ fine_cluster_by = function(seqs, by, max_dist = NULL, ...){

 #' aaseq = Biostrings::readAAStringSet(fasta_path)[1:100]

 #' cls = fine_cluster(aaseq)

 #' plot(cls$cluster)

-fine_cluster = function(seqs, type = 'AA', big_memory_brute = FALSE, method = 'levenshtein', substitution_matrix = 'BLOSUM100', cluster = 'hclust', cluster_method = 'complete'){

+fine_cluster = function(seqs, type = 'AA', big_memory_brute = FALSE, method = 'levenshtein', substitution_matrix = 'BLOSUM100', cluster_fun = 'hclust', cluster_method = 'complete'){

     if(length(seqs) > 4000 & !big_memory_brute) stop("Won't try to cluster ", length(seqs), " sequences unless `big_memory_brute` = TRUE.  (Make sure you actually have lots of memory!)")

     type = match.arg(type, choices = c('AA', 'DNA'))

-    cluster = match.arg(cluster, c('hclust', 'none'))

+    cluster_fun = match.arg(cluster_fun, c('hclust', 'none'))

     method = match.arg(method, c('substitutionMatrix', 'levenshtein'))

     if(is.character(seqs)){

@@ -122,7 +75,7 @@ fine_cluster = function(seqs, type = 'AA', big_memory_brute = FALSE, method = 'l

         #if(length(ss)>1) sd = sd/max(sd)*median(Biostrings::nchar(ss))

     }

-    if(length(seqs)>1 && cluster != 'none'){

+    if(length(seqs)>1 && cluster_fun != 'none'){

         hc = stats::hclust(as.dist(sd), method = cluster_method)

         hc$labels = names(ss)

     } else{

@@ -169,6 +122,7 @@ entropy = function(v, pseudo_count = length(v)/1000, na.action = na.fail){

 #' @export

 #' @describeIn entropy The number of categories exceeding `p` proportion of the total

+#' @param p proportion threshold

 np = function(v, p = .05, pseudo_count = p/5, na.action = na.fail){

     v = na.action(v)

     if(length(v) == 0) return(0)

new file mode 100644

@@ -0,0 +1,209 @@

+#' For each cell, return a single, canonical chain-cluster

+#'

+#' In single cell data, multiple chains (heavy-light or alpha-beta) are expected.  In some cases, there could be more than two (eg multiple alpha alleles for T cells).

+#' This picks a cluster id for each cell based on the overall prevalence of cluster ids over all cells in `tbl`.

+#' If order = 1 then the canonical chain-cluster will be the most prevalent, and if order = 2, it will be the 2nd most prevalent, and so on.  Ties are broken arbitrarily (possibly by lexicographic order of `cluster_idx`).

+#' @param tbl `data.frame` containing columns specified in `cell_identifiers`, `cluster_idx` and optionally `chain_identifiers`

+#' @param cell_identifiers `character` vector specifying columns in `tbl` that identify a cell

+#' @param cluster_idx `character` specifying the column in `tbl` that identifies a clsuter

+#' @param order return the 1st, 2nd, 3rd, etc, most common chain-cluster

+#'

+#' @return `data.frame` with columns from `cell_identifiers` and a single `cluster_idx` for each cell

+#' @export

+canonicalize_by_prevalence = function(tbl, cell_identifiers = 'barcode', cluster_idx = 'cluster_idx', order = 1){

+    prevalence = tbl %>% group_by(!!sym(cluster_idx)) %>% summarize(prev = dplyr::n())

+    tbl_prevalence = left_join(tbl, prevalence, by = cluster_idx)

+    tbl_order = tbl_prevalence %>% group_by(!!!syms(cell_identifiers)) %>% summarize(!!cluster_idx := dplyr::nth(!!sym(cluster_idx), -order, prev))

+    ungroup(tbl_order)

+}

+

+

+

+

+#' @param sort_factors `character` vector naming columns in `tbl` to sorted on, within  `cell_identifier`. Sorted by first element first, then ties broken by subsequent elements.  Sorted in decreasing order for each element.

+#' @param chain_levels an optional `character` vector providing the sort order of the `chain` column in `tbl`.  Set to length zero to disable.

+#' @export

+#' @describeIn canonicalize_by_prevalence return a canonical contig by chain type, with TRB/IGH returned first. By default, ties are broken by umis and reads.

+canonicalize_by_chain = function(tbl,  cell_identifiers = 'barcode', sort_factors = c('chain', 'umis', 'reads'), cluster_idx = 'cluster_idx', order = 1, chain_levels = c('IGL', 'IGK', 'TRA', 'TRB', 'IGH')){

+    #ochain = tbl[[sort_factors[1]]]

+    if(length(chain_levels) > 0){

+        if(!('chain' %in% names(tbl))) stop('`tbl` must contain a `chain` column to set `chain_levels`')

+        tbl = tbl %>% mutate(chain = factor(chain, levels = chain_levels, ordered = TRUE))

+    }

+    arranging = purrr::map(sort_factors, ~ rlang::quo(desc(!!sym(.x))))

+    tbl %>% group_by(!!!syms(cell_identifiers)) %>% dplyr::arrange(!!!arranging) %>% mutate(rank = seq_along(!!sym(sort_factors[[1]]))) %>% dplyr::filter(rank == order) %>% dplyr::select(-rank)

+

+}

+

+# Define canonical chain types per cell

+

+#' Given a family of similar sequences, return a "representative"

+#'

+#' @param seqs character vector

+#' @param medoid_idx optional index into seqs

+#' @param warn_if_distinct Should a warning be emitted if there are distinct elements in seqs?

+#'

+#' If medoid_idx is supplied, the medoid sequence is returned, otherwise the longest

+#' sequence is returned

+#'

+#' @return character vector

+#' @export

+#'

+#' @examples

+#' get_canonical_representative(c('apple', 'manzana', 'pomme'))

+get_canonical_representative = function(seqs, medoid_idx, warn_if_distinct = FALSE){

+    if(!missing(medoid_idx)){

+        if(!is.integer(medoid_idx) || medoid_idx < 1 || medoid_idx > length(seqs)) stop("Illegal `medoid_idx`")

+        ii = medoid_idx

+    } else{

+        len = str_length(seqs)

+        ii = which.max(len)

+    }

+    rep = seqs[ii]

+

+    if(warn_if_distinct && length(nomatch <- which(seqs != rep)) > 0) warning("At indices ", nomatch, " sequences did not match the representative")

+

+    return(rep)

+}

+

+#' Generate a list of tables representing clusters paired in cells

+#'

+#' A contingency table of every combination of `cluster_idx` up to `table_order` is generated.

+#' Combinations that are found in at least `min_expansion` number of cells are reported.  All cells that have these combinations are returned, as well as cells that only have `orphan_level` of matching `cluster_idx`.

+#'

+#' For example, if `table_order=2` and `min_expansion=2` then heavy/light or alpha/beta pairs found two or more times will be returned (as well as alpha-alpha pairs, etc, if those are present).

+#' If `orphan_level=1` then all cells that share just a single chain with an expanded clone will be returned.

+#'

+#' The `cluster_idx.1_fct` and `cluster_idx.2_fct` fields in `cell_tbl`, `idx1_tbl`, `idx2_tbl` are cast to factors and ordered such that pairings will tend to occur along the diagonal when they are cross-tabulated.

+#' This facilitates plotting.

+#'

+#' @section Caveats and warnings:

+#'  The cell_idx -> cluster_idx map is generally one-to-many, and is resolved by `canonicalize_fun`.  For `table_order>1`, few collisions are expected as most cells will contain no more than 2 clusters.  Any collisions are resolved by returning the most prevalent cluster, across samples.  When two clusters are tied for most prevalent within a cell, the `cluster_idx` returned is arbitrary. Therefore, when `table_order=1`, it is strongly recommended to subset the `cluster_tbl` to just a single chain.

+#'

+#' @param cluster_tbl a table with all combinations of clusters in all cells

+#' @param cell_identifiers character vector naming fields that key a cell

+#' @param cluster_idx character naming a single field IDing the clusters

+#' @param min_expansion the minimal number of times a pairing needs to occur for it to be reported

+#' @param cell_tbl optional, ancillary table with additional cell features.  Must also be keyed by `cell_identifiers`

+#' @param feature_tbl optional, ancillary table with additional cluster features.  Must also be keyed by `cluster_idx`

+#' @param table_order Integer larger than 1. What order of cluster_idx will be paired, eg, order = 2 means that the most common and second most common cluster_idx will be sought for each cell

+#' @param orphan_level Integer larger than 0 and less than or equal to `table_order`.  Given that at least `min_expansion` cells are found that have `table_order` chains identical, how many `cluster_idx` pairs will we match on to select other cells.  Example: `ophan_level=1` means that cells that share just a single chain with the

+#' @param cluster_whitelist a table of "cluster_idx" that should always be reported.  In contrast to the `cluster_tbl`, here the clusters must be named "cluster_idx.1", "cluster_idx.2" (if order-2 pairs are being selected).

+#' @param cluster_blacklist a table of "cluster_idx" that will never be reported.  Must be named as per `cluster_whitelist`.

+#' @param canonicalize_fun a function with signature `canonicalize_fun(cluster_tbl, cell_identifiers, cluster_idx, order = i)` that for each `cell_identifier` returns a single contig that depends on the `order`.  For instance \link{canonicalize_by_prevalence} or \link{canonicalize_by_chain}.

+#'

+#' @return list of tables.  The `cell_tbl` is keyed by the `cell_identifiers`, with fields "cluster_idx.1", "cluster_idx.2", etc, IDing the contigs present in each cell. "cluster_idx.1_fct" and "cluster_idx.2_fct" cast these fields to factors and are reordered to maximize the number of pairs along the diagonal. The `idx1_tbl` and `idx2_tbl` report information (passed in about the `cluster_idx` by `feature_tbl`.)  The `cluster_pair_tbl` reports all pairings found of contigs, and the number of times observed.

+#' @export

+#'

+#' @seealso canonicalize_by_prevalence, canonicalize_by_chain

+#' @importFrom tibble as_data_frame

+#' @importFrom dplyr bind_rows left_join ungroup summarize anti_join

+#' @importFrom stringr str_length str_c

+#' @importFrom rlang sym syms :=

+#' @examples

+#' library(dplyr)

+#' cluster_tbl = data_frame(clust_idx = gl(3, 2), cell_idx = rep(1:3, times = 2))

+#' # no pairs found twice

+#' pt1 = pairing_tables(cluster_tbl, 'cell_idx', 'clust_idx', canonicalize_by_prevalence)

+#' # all pairs found, found once.

+#' pt2 = pairing_tables(cluster_tbl, 'cell_idx', 'clust_idx',

+#'     canonicalize_by_prevalence, min_expansion = 1)

+#' pt2$cell_tbl

+#' cluster_tbl2 = bind_rows(cluster_tbl, cluster_tbl %>% mutate(cell_idx = rep(4:6, times = 2)))

+#' #all pairs found twice

+#' pt3 = pairing_tables(cluster_tbl2, 'cell_idx', 'clust_idx', canonicalize_by_prevalence, min_expansion = 1)

+#' pt3$cell_tbl

+#' # `canonicalize_by_chain` expects fields `umis`, `reads`

+#' # to break ties,  wrap the function to change this

+#' cluster_tbl3 = cluster_tbl2 %>%

+#'     mutate(umis = 1, reads = 1, chain = rep(c('TRA', 'TRB'), times = 6))

+#' pt4 = pairing_tables(cluster_tbl3, 'cell_idx', 'clust_idx',

+#'     canonicalize_by_chain, min_expansion = 1, table_order = 2)

+pairing_tables = function(cluster_tbl, cell_identifiers = 'barcode', cluster_idx = 'cluster_idx', canonicalize_fun = canonicalize_by_chain, table_order = 2, min_expansion = 2,  orphan_level = 1, cluster_whitelist = NULL, cluster_blacklist = NULL, cell_tbl = NULL, feature_tbl = NULL ){

+

+    if(orphan_level > table_order) stop('`ophan_level` must be less than or equal to `table_order`')

+    if(table_order < 1) stop('Table order must be at least 1')

+

+    # get `table_order` most common clusters for each cell

+    # forcibly rename cluster_idx -> "cluster_idx"

+    bar_chain_tbls = purrr::map(seq_len(table_order), function(i) canonicalize_fun(cluster_tbl, cell_identifiers = cell_identifiers, cluster_idx = cluster_idx, order = i) %>% dplyr::select(!!!c(syms(cell_identifiers), rlang::quo(cluster_idx))) %>% dplyr::rename( !!paste0('cluster_idx.', i) := !!cluster_idx))

+    # for each cell, what clusters are present

+    oligo_cluster_pairs = purrr::reduce(bar_chain_tbls, left_join, by = cell_identifiers)

+

+    # Set up cluster_ids for indexing into the pairing tables

+    cluster_ids = str_c('cluster_idx.', seq_len(table_order))

+    cluster_ids_to_select = cluster_ids[seq_len(orphan_level)]

+

+    # In how many cells do each cluster pairing appear?

+    cluster_pair_tbl = oligo_cluster_pairs %>% group_by(!!!syms(cluster_ids)) %>% summarize(n_clone_pairs = dplyr::n())

+    # which clusters are expanded

+    expanded_cluster = cluster_pair_tbl %>% dplyr::filter(n_clone_pairs >= min_expansion) %>% dplyr::filter_at(.vars = cluster_ids, .vars_predicate = all_vars(!is.na(.)))

+    expanded_cluster = ungroup(expanded_cluster) %>% dplyr::select(!!!syms(cluster_ids_to_select), max_pairs = n_clone_pairs)

+    if(!is.null(cluster_whitelist)){

+        expanded_cluster = bind_rows(expanded_cluster, cluster_whitelist)

+    }

+    if(!is.null(cluster_blacklist)) expanded_cluster = anti_join(expanded_cluster, cluster_blacklist)

+

+    # Could have duplicated cluster_ids after binding to the whitelist or from considering orphans

+    expanded_cluster = expanded_cluster[!duplicated(expanded_cluster %>% dplyr::select(-max_pairs)),]

+    expanded_c1 = oligo_cluster_pairs %>% dplyr::inner_join(expanded_cluster, by = cluster_ids_to_select)

+    if(anyDuplicated(expanded_c1 %>% dplyr::select(!!!syms(cell_identifiers)))) stop("Ruhoh, duplicated cell identifiers, this is a bug!")

+

+    # cross-tab pairings in order to figure out how to order the cluster_idx to put most common pairing on diagonal

+    if(nrow(expanded_c1)>0 && table_order > 1){

+        expanded_counts = reshape2::dcast(expanded_c1, cluster_idx.1 ~ cluster_idx.2, fun.aggregate = length, value.var = 'max_pairs')

+        mat = expanded_counts[,-1]

+        #rownames(mat) = rowid[['cluster_idx']]

+        ro = hclust(dist(mat))$order

+        co = hclust(dist(t(mat)))$order

+    } else if (nrow(expanded_c1) == 0) {

+        warning('No pairs found')

+        return()

+    } else{

+        expanded_counts = reshape2::dcast(expanded_c1, cluster_idx.1 ~ ., fun.aggregate = length, value.var = 'max_pairs')

+        ro = order(expanded_counts[[2]])

+        co = NA_integer_

+    }

+

+    ci_class = class(cluster_tbl[[cluster_idx]])

+    as_method = if(ci_class == 'factor') as.factor else function(x) as(x, ci_class)

+    rowid = data_frame(cluster_idx = expanded_counts$cluster_idx.1 %>% as_method, plot_order = ro)

+    colid = suppressWarnings(data_frame(cluster_idx = colnames(expanded_counts)[-1] %>% as_method, plot_order = co))

+    rowid[['cluster_idx.1_fct']] = factor(rowid[['cluster_idx']], levels = rowid[['cluster_idx']][ro])

+    colid[['cluster_idx.2_fct']] = factor(colid[['cluster_idx']], levels = colid[['cluster_idx']][co])

+

+

+    # also fix levels of this tbl

+    expanded_c1 = expanded_c1 %>% mutate(cluster_idx.1_fct = factor(cluster_idx.1, levels = levels(rowid[['cluster_idx.1_fct']])))

+    if(table_order>1) expanded_c1 = expanded_c1 %>% mutate(cluster_idx.2_fct = factor(cluster_idx.2, levels = levels(colid[['cluster_idx.2_fct']])))

+

+    if(!is.null(cell_tbl)){

+        if(anyDuplicated(cell_tbl %>% select(!!!syms(cell_identifiers)))) stop('`cell_tbl` must not have duplicate `cell_identifiers`.')

+        expanded_c1 = left_join(expanded_c1, cell_tbl, by = cell_identifiers)

+    }

+

+    idx1_tbl = rowid %>% dplyr::rename(!!cluster_idx := cluster_idx)

+    idx2_tbl = colid %>% dplyr::rename(!!cluster_idx := cluster_idx)

+    if(!is.null(feature_tbl)){

+        if(anyDuplicated(feature_tbl %>% select(!!sym(cluster_idx)))) stop('`feature_tbl` must not have duplicate `cluster_idx`.')

+        idx1_tbl = left_join(idx1_tbl, feature_tbl, by = cluster_idx)

+        idx2_tbl = left_join(idx2_tbl, feature_tbl, by = cluster_idx)

+    }

+

+    list(cell_tbl = expanded_c1, idx1_tbl = idx1_tbl, idx2_tbl = idx2_tbl, cluster_pair_tbl = cluster_pair_tbl)

+

+}

+

+plot_pairing = function(pairing_list, color_labels_by){

+    pl = pairing_list

+    pairs_plt = ggplot(pairing_list$cell_tbl, aes(x = cluster_idx.1_fct, y = cluster_idx.2_fct, color = sample, shape = pop)) + geom_jitter(width = .3, height = .3)

+

+    ylab = data_frame(!!color_labels_by :=  ggplot_build(pairs_plt)$layout$panel_params[[1]]$y.label) %>% left_join(feature_tbl) %>% mutate(class_color = ifelse(is.na(class_color), '#E41A1C', class_color))

+

+    xlab = data_frame(!!color_labels_by :=  ggplot_build(pairs_plt)$layout$panel_params[[1]]$x.label) %>% left_join(feature_tbl) %>% mutate(class_color = ifelse(is.na(class_color), '#E41A1C', class_color))

+

+    pairs_plt = pairs_plt + theme(axis.text.x = element_text(angle = 90, color = xlab$class_color, size = 8), axis.text.y = element_text(color = ylab$class_color, size = 8))

+

+

+}

+

@@ -12,9 +12,11 @@

 #'

 #' @examples

 #' library(dplyr)

-#' contig_anno_path = system.file('extdata', 'cellranger_contig_annotation.csv', package = 'CellaRepertorium')

+#' contig_anno_path = system.file('extdata', 'cellranger_contig_annotation.csv',

+#'     package = 'CellaRepertorium')

 #' contig_anno = readr::read_csv(contig_anno_path)

-#' contig_anno = contig_anno %>% mutate(fancy_name = fancy_name_contigs(., prefix = paste(sample, pop, sep = '_')))

+#' contig_anno = contig_anno %>% mutate(fancy_name =

+#'     fancy_name_contigs(., prefix = paste(sample, pop, sep = '_')))

 #' stopifnot(!any(duplicated(contig_anno$fancy_name)))

 fancy_name_contigs = function(contig_frame, prefix){

     notin = setdiff(c('chain', 'v_gene', 'd_gene', 'j_gene'), names(contig_frame))

@@ -31,138 +33,6 @@ variable_genes = function(ref = '10X'){

 }

-#' Generate a list of tables representing clusters paired in cells

-#'

-#' A contingency table of every combination of `cluster_idx` up to `table_order` is generated.

-#' Combinations that are found in at least `min_expansion` number of cells are reported.  All cells that have these combinations are returned, as well as cells that only have `orphan_level` of matching `cluster_idx`.

-#'

-#' For example, if `table_order=2` and `min_expansion=2` then heavy/light or alpha/beta pairs found two or more times will be returned (as well as alpha-alpha pairs, etc, if those are present).

-#' If `orphan_level=1` then all cells that share just a single chain with an expanded clone will be returned.

-#'

-#' The `cluster_idx.1_fct` and `cluster_idx.2_fct` fields in `cell_tbl`, `idx1_tbl`, `idx2_tbl` are cast to factors and ordered such that pairings will tend to occur along the diagonal when they are cross-tabulated.

-#' This facilitates plotting.

-#'

-#' @section Caveats and warnings:

-#'  The cell_idx -> cluster_idx map is generally one-to-many.  For `table_order>1`, few collisions are expected as most cells will contain no more than 2 clusters.  Any collisions are resolved by returning the most prevalent cluster, across samples.  When two clusters are tied for most prevalent within a cell, the `cluster_idx` returned is arbitrary. Therefore, when `table_order=1`, it is strongly recommended to subset the `cluster_tbl` to just a single chain.

-#' @param cluster_tbl a table with all combinations of clusters in all cells

-#' @param cell_identifiers character vector naming fields that key a cell

-#' @param cluster_idx character naming a single field IDing the clusters

-#' @param min_expansion the minimal number of times a pairing needs to occur for it to be reported

-#' @param cell_tbl optional, ancillary table with additional cell features.  Must also be keyed by `cell_identifiers`

-#' @param feature_tbl optional, ancillary table with additional cluster features.  Must also be keyed by `cluster_idx`

-#' @param table_order Integer larger than 1. What order of cluster_idx will be paired, eg, order = 2 means that the most common and second most common cluster_idx will be sought for each cell

-#' @param orphan_level Integer larger than 0 and less than or equal to `table_order`.  Given that at least `min_expansion` cells are found that have `table_order` chains identical, how many `cluster_idx` pairs will we match on to select other cells.  Example: `ophan_level=1` means that cells that share just a single chain with the

-#' @param cluster_whitelist a table of "cluster_idx" that should always be reported.  In contrast to the `cluster_tbl`, here the clusters must be named "cluster_idx.1", "cluster_idx.2" (if order-2 pairs are being selected).

-#' @param cluster_blacklist a table of "cluster_idx" that will never be reported.  Must be named as per `cluster_whitelist`.

-#'

-#' @return list of tables.  The `cell_tbl` is keyed by the `cell_identifiers`, with fields "cluster_idx.1", "cluster_idx.2", etc, IDing the contigs present in each cell. "cluster_idx.1_fct" and "cluster_idx.2_fct" cast these fields to factors and are reordered to maximize the number of pairs along the diagonal. The `idx1_tbl` and `idx2_tbl` report information (passed in about the `cluster_idx` by `feature_tbl`.)  The `cluster_pair_tbl` reports all pairings found of contigs, and the number of times observed.

-#' @export

-#'

-#' @seealso get_canonical_chain

-#' @importFrom tibble as_data_frame

-#' @importFrom dplyr bind_rows left_join ungroup summarize anti_join

-#' @importFrom stringr str_length str_c

-#' @importFrom rlang sym syms :=

-#' @examples

-#' library(dplyr)

-#' cluster_tbl = data_frame(clust_idx = gl(3, 2), cell_idx = rep(1:3, times = 2))

-#' # no pairs found twice

-#' pairing_tables(cluster_tbl, 'cell_idx', 'clust_idx')

-#' # all pairs found, found once.

-#' pairing_tables(cluster_tbl, 'cell_idx', 'clust_idx', min_expansion = 1)

-#' cluster_tbl2 = bind_rows(cluster_tbl, cluster_tbl %>% mutate(cell_idx = rep(4:6, times = 2)))

-#' #all pairs found twice

-#' pairing_tables(cluster_tbl2, 'cell_idx', 'clust_idx', min_expansion = 1)

-#' # A warning if `table_order = 1` -- best to stratify by chain

-#' # (could be accomplished by setting the `cell_identifiers` to include a field that identifies the chain).

-#' pairing_tables(cluster_tbl2, 'cell_idx', 'clust_idx', min_expansion = 1, table_order = 1)

-pairing_tables = function(cluster_tbl, cell_identifiers = 'barcode', cluster_idx = 'cluster_idx', table_order = 2, min_expansion = 2,  orphan_level = 1, cluster_whitelist = NULL, cluster_blacklist = NULL, cell_tbl = NULL, feature_tbl = NULL ){

-

-    if(orphan_level > table_order) stop('`ophan_level` must be less than or equal to `table_order`')

-    if(table_order < 1) stop('Table order must be at least 1')

-

-    # get `table_order` most common clusters for each cell

-    # forcibly rename cluster_idx -> "cluster_idx"

-    bar_chain_tbls = purrr::map(seq_len(table_order), function(i) get_canonical_chain(cluster_tbl, cell_identifiers = cell_identifiers, cluster_idx = cluster_idx, order = i) %>% dplyr::rename( !!paste0('cluster_idx.', i) := !!cluster_idx))

-    # for each cell, what clusters are present

-    oligo_cluster_pairs = purrr::reduce(bar_chain_tbls, left_join, by = cell_identifiers)

-

-    # Set up cluster_ids for indexing into the pairing tables

-    cluster_ids = str_c('cluster_idx.', seq_len(table_order))

-    cluster_ids_to_select = cluster_ids[seq_len(orphan_level)]

-

-    # In how many cells do each cluster pairing appear?

-    cluster_pair_tbl = oligo_cluster_pairs %>% group_by(!!!syms(cluster_ids)) %>% summarize(n_clone_pairs = n())

-    # which clusters are expanded

-    expanded_cluster = cluster_pair_tbl %>% dplyr::filter(n_clone_pairs >= min_expansion) %>% dplyr::filter_at(.vars = cluster_ids_to_select, .vars_predicate = all_vars(!is.na(.)))

-    expanded_cluster = ungroup(expanded_cluster) %>% dplyr::select(!!!syms(cluster_ids_to_select), max_pairs = n_clone_pairs)

-    if(!is.null(cluster_whitelist)){

-        expanded_cluster = bind_rows(expanded_cluster, cluster_whitelist)

-    }

-    if(!is.null(cluster_blacklist)) expanded_cluster = anti_join(expanded_cluster, cluster_blacklist)

-

-    # Could have duplicated cluster_ids after binding to the whitelist or from considering orphans

-    expanded_cluster = expanded_cluster[!duplicated(expanded_cluster %>% dplyr::select(-max_pairs)),]

-    expanded_c1 = oligo_cluster_pairs %>% right_join(expanded_cluster, by = cluster_ids_to_select)

-    if(anyDuplicated(expanded_c1 %>% dplyr::select(!!!syms(cell_identifiers)))) stop("Ruhoh, duplicated cell identifiers, this is a bug!")

-

-    # cross-tab pairings in order to figure out how to order the cluster_idx to put most common pairing on diagonal

-    if(nrow(expanded_c1)>0 && table_order > 1){

-        expanded_counts = reshape2::dcast(expanded_c1, cluster_idx.1 ~ cluster_idx.2, fun.aggregate = length, value.var = 'max_pairs')

-        mat = expanded_counts[,-1]

-        #rownames(mat) = rowid[['cluster_idx']]

-        ro = hclust(dist(mat))$order

-        co = hclust(dist(t(mat)))$order

-    } else if (nrow(expanded_c1) == 0) {

-       warning('No pairs found')

-        return()

-    } else{

-            expanded_counts = reshape2::dcast(expanded_c1, cluster_idx.1 ~ ., fun.aggregate = length, value.var = 'max_pairs')

-            ro = order(expanded_counts[[2]])

-            co = NA_integer_

-    }

-

-    ci_class = class(cluster_tbl[[cluster_idx]])

-    as_method = if(ci_class == 'factor') as.factor else function(x) as(x, ci_class)

-    rowid = data_frame(cluster_idx = expanded_counts$cluster_idx.1 %>% as_method, plot_order = ro)

-    colid = suppressWarnings(data_frame(cluster_idx = colnames(expanded_counts)[-1] %>% as_method, plot_order = co))

-    rowid[['cluster_idx.1_fct']] = factor(rowid[['cluster_idx']], levels = rowid[['cluster_idx']][ro])

-    colid[['cluster_idx.2_fct']] = factor(colid[['cluster_idx']], levels = colid[['cluster_idx']][co])

-

-

-    # also fix levels of this tbl

-    expanded_c1 = expanded_c1 %>% mutate(cluster_idx.1_fct = factor(cluster_idx.1, levels = levels(rowid[['cluster_idx.1_fct']])))

-    if(table_order>1) expanded_c1 = expanded_c1 %>% mutate(cluster_idx.2_fct = factor(cluster_idx.2, levels = levels(colid[['cluster_idx.2_fct']])))

-

-    if(!is.null(cell_tbl)){

-        if(anyDuplicated(cell_tbl %>% select(!!!syms(cell_identifiers)))) stop('`cell_tbl` must not have duplicate `cell_identifiers`.')

-        expanded_c1 = left_join(expanded_c1, cell_tbl, by = cell_identifiers)

-    }

-

-    idx1_tbl = rowid %>% dplyr::rename(!!cluster_idx := cluster_idx)

-    idx2_tbl = colid %>% dplyr::rename(!!cluster_idx := cluster_idx)

-    if(!is.null(feature_tbl)){

-        if(anyDuplicated(feature_tbl %>% select(!!sym(cluster_idx)))) stop('`feature_tbl` must not have duplicate `cluster_idx`.')

-        idx1_tbl = left_join(idx1_tbl, feature_tbl, by = cluster_idx)

-        idx2_tbl = left_join(idx2_tbl, feature_tbl, by = cluster_idx)

-    }

-

-    list(cell_tbl = expanded_c1, idx1_tbl = idx1_tbl, idx2_tbl = idx2_tbl, cluster_pair_tbl = cluster_pair_tbl)

-

-}

-

-plot_pairing = function(pairing_list, color_labels_by){

-    pl = pairing_list

-    pairs_plt = ggplot(pairing_list$cell_tbl, aes(x = cluster_idx.1_fct, y = cluster_idx.2_fct, color = sample, shape = pop)) + geom_jitter(width = .3, height = .3)

-

-    ylab = data_frame(!!color_labels_by :=  ggplot_build(pairs_plt)$layout$panel_params[[1]]$y.label) %>% left_join(feature_tbl) %>% mutate(class_color = ifelse(is.na(class_color), '#E41A1C', class_color))

-

-    xlab = data_frame(!!color_labels_by :=  ggplot_build(pairs_plt)$layout$panel_params[[1]]$x.label) %>% left_join(feature_tbl) %>% mutate(class_color = ifelse(is.na(class_color), '#E41A1C', class_color))

-

-    pairs_plt = pairs_plt + theme(axis.text.x = element_text(angle = 90, color = xlab$class_color, size = 8), axis.text.y = element_text(color = ylab$class_color, size = 8))

-

-

-}

 cleanup_annotations = function(json){

     anno = json[['annotations']]

new file mode 100644

@@ -0,0 +1,4 @@

+# External data

+T and B cells from healthy PBMC donor, V3 chemistry

+https://support.10xgenomics.com/single-cell-vdj/datasets/3.0.0/vdj_v1_hs_pbmc2_t

+

new file mode 100644

@@ -0,0 +1,41 @@

+% Generated by roxygen2: do not edit by hand

+% Please edit documentation in R/pairing-methods.R

+\name{canonicalize_by_prevalence}

+\alias{canonicalize_by_prevalence}

+\alias{canonicalize_by_chain}

+\title{For each cell, return a single, canonical chain-cluster}

+\usage{

+canonicalize_by_prevalence(tbl, cell_identifiers = "barcode",

+  cluster_idx = "cluster_idx", order = 1)

+

+canonicalize_by_chain(tbl, cell_identifiers = "barcode",

+  sort_factors = c("chain", "umis", "reads"),

+  cluster_idx = "cluster_idx", order = 1, chain_levels = c("IGL",

+  "IGK", "TRA", "TRB", "IGH"))

+}

+\arguments{

+\item{tbl}{`data.frame` containing columns specified in `cell_identifiers`, `cluster_idx` and optionally `chain_identifiers`}

+

+\item{cell_identifiers}{`character` vector specifying columns in `tbl` that identify a cell}

+

+\item{cluster_idx}{`character` specifying the column in `tbl` that identifies a clsuter}

+

+\item{order}{return the 1st, 2nd, 3rd, etc, most common chain-cluster}

+

+\item{sort_factors}{`character` vector naming columns in `tbl` to sorted on, within  `cell_identifier`. Sorted by first element first, then ties broken by subsequent elements.  Sorted in decreasing order for each element.}

+

+\item{chain_levels}{an optional `character` vector providing the sort order of the `chain` column in `tbl`.  Set to length zero to disable.}

+}

+\value{

+`data.frame` with columns from `cell_identifiers` and a single `cluster_idx` for each cell

+}

+\description{

+In single cell data, multiple chains (heavy-light or alpha-beta) are expected.  In some cases, there could be more than two (eg multiple alpha alleles for T cells).

+This picks a cluster id for each cell based on the overall prevalence of cluster ids over all cells in `tbl`.

+If order = 1 then the canonical chain-cluster will be the most prevalent, and if order = 2, it will be the 2nd most prevalent, and so on.  Ties are broken arbitrarily (possibly by lexicographic order of `cluster_idx`).

+}

+\section{Functions}{

+\itemize{

+\item \code{canonicalize_by_chain}: return a canonical contig by chain type, with TRB/IGH returned first. By default, ties are broken by umis and reads.

+}}

+

@@ -34,7 +34,8 @@ For example, \code{clusters} could be a clonal ID, while \code{labels} is a subj

 \examples{

 library(dplyr)

 purity = function(clusters, labels, covariates){

-n_label_cluster = data_frame(labels = labels, clusters = clusters) \%>\% group_by(clusters, labels) \%>\% summarize(n = n()) \%>\% ungroup()

+n_label_cluster = data_frame(labels = labels, clusters = clusters) \%>\%

+    group_by(clusters, labels) \%>\% summarize(n = n()) \%>\% ungroup()

 singletons = mean(n_label_cluster$n == 1)

 singletons

 }

@@ -18,6 +18,8 @@ modal_category(v, na.action = na.fail)

 \item{pseudo_count}{number of pseudo counts to add on, to stablize empty categories}

 \item{na.action}{how to handle NA values}

+

+\item{p}{proportion threshold}

 }

 \value{

 the sample entropy

@@ -21,6 +21,7 @@ Generate a legible name for a series of contigs

 library(dplyr)

 contig_anno_path = system.file('extdata', 'cellranger_contig_annotation.csv', package = 'CellaRepertorium')

 contig_anno = readr::read_csv(contig_anno_path)

-contig_anno = contig_anno \%>\% mutate(fancy_name = fancy_name_contigs(., prefix = paste(sample, pop, sep = '_')))

+contig_anno = contig_anno \%>\% mutate(fancy_name =

+    fancy_name_contigs(., prefix = paste(sample, pop, sep = '_')))

 stopifnot(!any(duplicated(contig_anno$fancy_name)))

 }

new file mode 100644

@@ -0,0 +1,27 @@

+% Generated by roxygen2: do not edit by hand

+% Please edit documentation in R/pairing-methods.R

+\name{get_canonical_representative}

+\alias{get_canonical_representative}

+\title{Given a family of similar sequences, return a "representative"}

+\usage{

+get_canonical_representative(seqs, medoid_idx, warn_if_distinct = FALSE)

+}

+\arguments{

+\item{seqs}{character vector}

+

+\item{medoid_idx}{optional index into seqs}

+

+\item{warn_if_distinct}{Should a warning be emitted if there are distinct elements in seqs?

+

+If medoid_idx is supplied, the medoid sequence is returned, otherwise the longest

+sequence is returned}

+}

+\value{

+character vector

+}

+\description{

+Given a family of similar sequences, return a "representative"

+}

+\examples{

+get_canonical_representative(c('apple', 'manzana', 'pomme'))

+}

new file mode 100644

@@ -0,0 +1,76 @@

+% Generated by roxygen2: do not edit by hand

+% Please edit documentation in R/pairing-methods.R

+\name{pairing_tables}

+\alias{pairing_tables}

+\title{Generate a list of tables representing clusters paired in cells}

+\usage{

+pairing_tables(cluster_tbl, cell_identifiers = "barcode",

+  cluster_idx = "cluster_idx",

+  canonicalize_fun = canonicalize_by_chain, table_order = 2,

+  min_expansion = 2, orphan_level = 1, cluster_whitelist = NULL,

+  cluster_blacklist = NULL, cell_tbl = NULL, feature_tbl = NULL)

+}

+\arguments{

+\item{cluster_tbl}{a table with all combinations of clusters in all cells}

+

+\item{cell_identifiers}{character vector naming fields that key a cell}

+

+\item{cluster_idx}{character naming a single field IDing the clusters}

+

+\item{canonicalize_fun}{a function with signature `canonicalize_fun(cluster_tbl, cell_identifiers, cluster_idx, order = i)` that for each `cell_identifier` returns a single contig that depends on the `order`.  For instance \link{canonicalize_by_prevalence} or \link{canonicalize_by_chain}.}

+

+\item{table_order}{Integer larger than 1. What order of cluster_idx will be paired, eg, order = 2 means that the most common and second most common cluster_idx will be sought for each cell}

+

+\item{min_expansion}{the minimal number of times a pairing needs to occur for it to be reported}

+

+\item{orphan_level}{Integer larger than 0 and less than or equal to `table_order`.  Given that at least `min_expansion` cells are found that have `table_order` chains identical, how many `cluster_idx` pairs will we match on to select other cells.  Example: `ophan_level=1` means that cells that share just a single chain with the}

+

+\item{cluster_whitelist}{a table of "cluster_idx" that should always be reported.  In contrast to the `cluster_tbl`, here the clusters must be named "cluster_idx.1", "cluster_idx.2" (if order-2 pairs are being selected).}

+

+\item{cluster_blacklist}{a table of "cluster_idx" that will never be reported.  Must be named as per `cluster_whitelist`.}

+

+\item{cell_tbl}{optional, ancillary table with additional cell features.  Must also be keyed by `cell_identifiers`}

+

+\item{feature_tbl}{optional, ancillary table with additional cluster features.  Must also be keyed by `cluster_idx`}

+}

+\value{

+list of tables.  The `cell_tbl` is keyed by the `cell_identifiers`, with fields "cluster_idx.1", "cluster_idx.2", etc, IDing the contigs present in each cell. "cluster_idx.1_fct" and "cluster_idx.2_fct" cast these fields to factors and are reordered to maximize the number of pairs along the diagonal. The `idx1_tbl` and `idx2_tbl` report information (passed in about the `cluster_idx` by `feature_tbl`.)  The `cluster_pair_tbl` reports all pairings found of contigs, and the number of times observed.

+}

+\description{

+A contingency table of every combination of `cluster_idx` up to `table_order` is generated.

+Combinations that are found in at least `min_expansion` number of cells are reported.  All cells that have these combinations are returned, as well as cells that only have `orphan_level` of matching `cluster_idx`.

+}

+\details{

+For example, if `table_order=2` and `min_expansion=2` then heavy/light or alpha/beta pairs found two or more times will be returned (as well as alpha-alpha pairs, etc, if those are present).

+If `orphan_level=1` then all cells that share just a single chain with an expanded clone will be returned.

+

+The `cluster_idx.1_fct` and `cluster_idx.2_fct` fields in `cell_tbl`, `idx1_tbl`, `idx2_tbl` are cast to factors and ordered such that pairings will tend to occur along the diagonal when they are cross-tabulated.

+This facilitates plotting.

+}

+\section{Caveats and warnings}{

+

+ The cell_idx -> cluster_idx map is generally one-to-many, and is resolved by `canonicalize_fun`.  For `table_order>1`, few collisions are expected as most cells will contain no more than 2 clusters.  Any collisions are resolved by returning the most prevalent cluster, across samples.  When two clusters are tied for most prevalent within a cell, the `cluster_idx` returned is arbitrary. Therefore, when `table_order=1`, it is strongly recommended to subset the `cluster_tbl` to just a single chain.

+}

+

+\examples{

+library(dplyr)

+cluster_tbl = data_frame(clust_idx = gl(3, 2), cell_idx = rep(1:3, times = 2))

+# no pairs found twice

+pt1 = pairing_tables(cluster_tbl, 'cell_idx', 'clust_idx', canonicalize_by_prevalence)

+# all pairs found, found once.

+pt2 = pairing_tables(cluster_tbl, 'cell_idx', 'clust_idx',

+    canonicalize_by_prevalence, min_expansion = 1)

+pt2$cell_tbl

+cluster_tbl2 = bind_rows(cluster_tbl, cluster_tbl \%>\% mutate(cell_idx = rep(4:6, times = 2)))

+#all pairs found twice

+pt3 = pairing_tables(cluster_tbl2, 'cell_idx', 'clust_idx', canonicalize_by_prevalence, min_expansion = 1)

+pt3$cell_tbl

+# canonicalize_by_chain by default expects fields umis, reads to break ties, could wrap the function to change this

+cluster_tbl3 = cluster_tbl2 \%>\%

+    mutate(umis = 1, reads = 1, chain = rep(c('TRA', 'TRB'), times = 6))

+pt4 = pairing_tables(cluster_tbl3, 'cell_idx', 'clust_idx',

+    canonicalize_by_chain, min_expansion = 1, table_order = 2)

+}

+\seealso{

+canonicalize_by_prevalence, canonicalize_by_chain

+}