Fix "inconsistent counts" issue

With this fix, the following counts are now correct:

    db_name <- install_IMGT_germline_db("202614-2", "Homo_sapiens",
                                        tcr.db=TRUE)
    all_counts <- list_germline_dbs(long.listing=TRUE)  # long listing
    colSums(all_counts[[db_name]])
    #   V   D   J
    # 355   6 101

More generally speaking, this fix makes sure that the counts reported
by list_germline_dbs()'s and list_c_region_dbs()'s long listings are
consistent with the counts reported by their short listings.

Author: hpages

R/clean_allele_set.R

@@ -142,7 +142,7 @@ list_loci_in_germline_fasta_dir <-
 ### Create the three "region dbs": one V-, one D-, and one J-region db.
 .do_create_germline_db <- function(destdir, fasta_dir, loci,
                                    gapped=FALSE, with.intdata=FALSE,
-                                   excluded_J_alleles=NULL,
+                                   excluded_J_alleles=character(0),
                                    with.auxdata=FALSE, imgt.fasta=FALSE,
                                    disambiguate.allele.names=FALSE,
                                    verbose=FALSE)
@@ -158,7 +158,7 @@ list_loci_in_germline_fasta_dir <-
                        verbose=verbose)
     J_fasta_files <- .collect_fasta_files(fasta_dir, "J", loci)
     create_J_region_db(J_fasta_files, destdir,
-                       excluded_J_alleles=excluded_J_alleles,
+                       excluded_alleles=excluded_J_alleles,
                        with.auxdata=with.auxdata, imgt.fasta=imgt.fasta,
                        disambiguate.allele.names=disambiguate.allele.names,
                        verbose=verbose)
@@ -176,7 +176,7 @@ list_loci_in_germline_fasta_dir <-
 ### and 'disambiguate.allele.names' arguments.
 create_germline_db <- function(destdir, fasta_dir, loci,
                                gapped=FALSE, with.intdata=FALSE,
-                               excluded_J_alleles=NULL,
+                               excluded_J_alleles=character(0),
                                with.auxdata=FALSE, imgt.fasta=FALSE,
                                disambiguate.allele.names=FALSE,
                                overwrite=FALSE, verbose=FALSE)

R/create_germline_db.R

@@ -309,6 +309,26 @@ display_data_frame_in_browser <- function(df)
 }
 
 
+### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+### read_version_file()
+###
+
+read_version_file <- function(dirpath)
+{
+    stopifnot(isSingleNonWhiteString(dirpath), dir.exists(dirpath))
+    version_path <- file.path(dirpath, "version")
+    if (!file.exists(version_path))
+        stop(wmsg("missing 'version' file in ", dirpath, "/"))
+    version <- readLines(version_path)
+    if (length(version) != 1L)
+        stop(wmsg("file '", version_path, "' should contain exactly one line"))
+    version <- trimws2(version)
+    if (version == "")
+        stop(wmsg("file '", version_path, "' contains only white spaces"))
+    version
+}
+
+
 ### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 ### list_fasta_files()
 ###
@@ -326,55 +346,6 @@ list_fasta_files <- function(path, pattern="\\.fasta$", full.names=TRUE,
 }
 
 
-### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-### copy_fasta_file()
-###
-
-### By "repeated" alleles we mean alleles with identical sequences **and**
-### names. Note that:
-### - we keep alleles with identical sequences but different names;
-### - we also keep alleles with identical names but different sequences.
-### Note that the .drop_repeated_alleles() function in R/clean_allele_set.R
-### implements the same semantic.
-copy_fasta_file <- function(infasta, outfasta, excluded_alleles=NULL,
-                            drop.repeated.alleles=FALSE)
-{
-    stopifnot(isSingleNonWhiteString(infasta),
-              !dir.exists(infasta), file.exists(infasta),
-              isSingleNonWhiteString(outfasta),
-              !dir.exists(outfasta), !file.exists(outfasta),
-              isTRUEorFALSE(drop.repeated.alleles))
-    if (!is.null(excluded_alleles) || drop.repeated.alleles) {
-        headers <- names(fasta.seqlengths(infasta))
-        allele_names <- clean_imgt_fasta_headers(headers)
-        drop_me <- logical(length(allele_names))
-        if (!is.null(excluded_alleles)) {
-            stopifnot(is.character(excluded_alleles))
-            drop_me <- drop_me | (allele_names %in% excluded_alleles)
-        }
-        if (drop.repeated.alleles) {
-            alleles <- readBStringSet(infasta)
-            stopifnot(identical(names(alleles), headers))
-            m1 <- match(alleles, alleles)
-            m2 <- match(allele_names, allele_names)
-            is_repeated <- selfmatchIntegerPairs(m1, m2) != seq_along(alleles)
-            drop_me <- drop_me | is_repeated
-        }
-        drop_idx <- which(drop_me)
-        if (length(drop_idx) != 0L) {
-            lines <- readLines(infasta)
-            record_starts <- grep("^>", lines)
-            stopifnot(identical(lines[record_starts], paste0(">", headers)))
-            record_ends <- c(tail(record_starts - 1L, n=-1L), length(lines))
-            record_bounds <- IRanges(record_starts, record_ends)
-            lines <- extractROWS(lines, record_bounds[-drop_idx])
-            return(writeLines(lines, outfasta))
-        }
-    }
-    stopifnot(file.copy(infasta, outfasta, overwrite=FALSE))
-}
-
-
 ### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 ### fasta_files_as_one_uncompressed_file()
 ###

R/create_region_db.R

@@ -109,7 +109,7 @@ install_IMGT_germline_db <- function(release, organism="Homo sapiens",
     if (organism == "Homo_sapiens" && loci_prefix == "IG") {
         excluded_J_alleles <- .EXCLUDED_IMGT_HUMAN_J_ALLELES
     } else {
-        excluded_J_alleles <- NULL
+        excluded_J_alleles <- character(0)
     }
 
     ## Create and install germline db.

R/file-utils.R

@@ -54,7 +54,7 @@ stop_on_existing_cached_germline_db <- function(db_name)
 ### Not exported!
 install_germline_db <- function(install_dir, db_name, fasta_dir, loci,
                                 gapped=FALSE, with.intdata=FALSE,
-                                excluded_J_alleles=NULL,
+                                excluded_J_alleles=character(0),
                                 with.auxdata=FALSE, imgt.fasta=FALSE,
                                 disambiguate.allele.names=FALSE,
                                 if.exists=c("error", "overwrite", "no-op"),

R/install_IMGT_germline_db.R

@@ -130,6 +130,65 @@ check_seqlens <- function(seqlens, varname)
 }
 
 
+### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
+### make_unique_allele_names()
+###
+
+### Generates a pool of unique suffixes made of lowercase letters.
+.make_pool_of_suffixes <- function(min_pool_size)
+{
+    max_pool_size <- (length(letters)**8 - 1) / (length(letters) - 1) - 1
+    if (min_pool_size > max_pool_size)
+        stop(wmsg("too many duplicated allele names"))
+    ans <- character(0)
+    for (i in 1:7) {
+        ans <- c(ans, mkAllStrings(letters, i))
+        if (length(ans) >= min_pool_size)
+            return(ans)
+    }
+    ## Should never happen because we checked for this condition earlier (see
+    ## above).
+    stop(wmsg("too many duplicated allele names"))
+}
+
+### Returns the suffixes in a named character vector parallel to 'x'.
+.make_disambiguation_suffixes <- function(x)
+{
+    stopifnot(is.character(x))
+    if (length(x) <= 1L) {
+        suffixes <- character(length(x))
+    } else {
+        oo <- order(x)
+        ordered_x <- x[oo]
+        ir <- IRanges(1L, runLength(Rle(ordered_x)))
+        pool_of_suffixes <- .make_pool_of_suffixes(max(width(ir)))
+        suffixes <- extractList(pool_of_suffixes, ir)  # CharacterList
+        suffixes[lengths(suffixes) == 1L] <- ""
+        rev_oo <- S4Vectors:::reverseIntegerInjection(oo, length(oo))
+        suffixes <- unlist(suffixes, use.names=FALSE)[rev_oo]
+    }
+    setNames(suffixes, names(x))
+}
+
+### Has its own tests in tests/testthat/test-utils.R!
+### We make the allele names unique by adding a disambiguation suffix to
+### them. Note that this is not guaranteed to make the allele names unique
+### because adding the suffixes can create clashes with some other allele
+### names that already looks like they have a disambiguation suffix.
+### So we check for that and raise an error if it happens.
+make_unique_allele_names <- function(allele_names, suffixes.only=FALSE)
+{
+    stopifnot(isTRUEorFALSE(suffixes.only))
+    suffixes <- .make_disambiguation_suffixes(allele_names)
+    ans <- add_suffix(allele_names, suffixes)
+    if (anyDuplicated(ans))
+        stop(wmsg("Our allele name disambiguation scheme doesn't work for ",
+                  "your set of alleles, sorry. Please file an issue at ",
+                  "https://github.com/HyrienLab/igblastr/issues"))
+    if (suffixes.only) suffixes else ans
+}
+
+
 ### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 ### interweave()
 ###
@@ -168,26 +227,6 @@ align_vectors_by_names <- function(vectors)
 }
 
 
-### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-### read_version_file()
-###
-
-read_version_file <- function(dirpath)
-{
-    stopifnot(isSingleNonWhiteString(dirpath), dir.exists(dirpath))
-    version_path <- file.path(dirpath, "version")
-    if (!file.exists(version_path))
-        stop(wmsg("missing 'version' file in ", dirpath, "/"))
-    version <- readLines(version_path)
-    if (length(version) != 1L)
-        stop(wmsg("file '", version_path, "' should contain exactly one line"))
-    version <- trimws2(version)
-    if (version == "")
-        stop(wmsg("file '", version_path, "' contains only white spaces"))
-    version
-}
-
-
 ### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
 ### rows_with_same_key_are_identical()
 ###

R/install_custom_germline_db.R

@@ -7,18 +7,6 @@ o More testing of combine_germline_dbs(). Expand the examples in its man
 o In addition to combine_germline_dbs(), we need combine_c_region_dbs()
   to combine C-region dbs.
 
-o Fix the "inconsistent counts" issue:
-
-    db_name <- install_IMGT_germline_db("202614-2", "Homo_sapiens", tcr.db=TRUE)
-    all_counts <- list_germline_dbs(long.listing=TRUE)  # long listing
-    colSums(all_counts[[db_name]])
-    #   V   D   J 
-    # 370   6 101 
-
-  These counts are wrong and inconsistent with the counts reported
-  by list_germline_dbs()'s short listing (V:355 D:6 J:101).
-  See .finish_cleanup() in R/clean_allele_set.R why.
-
 o Maybe replace 'without.intdata' and 'without.auxdata' arguments with
   'with.intdata' and 'with.auxdata'. With default of TRUE for 'with.intdata'
   and "auto" for 'with.auxdata' ("auto" meaning "include the computed auxdata

R/utils.R

@@ -20,6 +20,47 @@ test_that(".drop_repeated_alleles()", {
     expect_identical(as.character(object), toupper(y[-c(5:7, 10L)]))
 })
 
+test_that(".make_DORsummary()", {
+    make_DORsummary <- igblastr:::.make_DORsummary
+
+    expected <- data.frame(allele_name="A", suffix="")
+    expect_identical(make_DORsummary(c(A="ct")), expected)
+
+    x <- c(A="ct", A="ct")
+    expected <- data.frame(allele_name=names(x), suffix=c("", NA))
+    expect_identical(make_DORsummary(x), expected)
+    expect_identical(make_DORsummary(BStringSet(x)), expected)
+    expect_identical(make_DORsummary(DNAStringSet(x)), expected)
+
+    x <- c(A="ct", A="ct", A="ct")
+    expected <- data.frame(allele_name=names(x), suffix=c("", NA, NA))
+    expect_identical(make_DORsummary(x), expected)
+    expect_identical(make_DORsummary(BStringSet(x)), expected)
+    expect_identical(make_DORsummary(DNAStringSet(x)), expected)
+
+    x <- c(A="ct", B="", A="ct", A="ct", C="", D="ggg", C="")
+    expected <- data.frame(allele_name=names(x),
+                           suffix=c("", "", NA, NA, "", "", NA))
+    expect_identical(make_DORsummary(x), expected)
+    expect_identical(make_DORsummary(BStringSet(x)), expected)
+    expect_identical(make_DORsummary(DNAStringSet(x)), expected)
+
+    x <- c(A="tt", B="g", C="gg", B="", C="g", A="tt", D="t",
+           B="ggg", B="", C="", B="g", C="gg", E="tt", B="tt")
+    expected <- data.frame(allele_name=names(x),
+                           suffix=c("", "a", "a", "b", "b", NA, "",
+                                    "c", NA, "c", NA, NA, "", "d"))
+    expect_identical(
+        make_DORsummary(x, disambiguate.allele.names=TRUE),
+        expected)
+    expect_identical(
+        make_DORsummary(BStringSet(x), disambiguate.allele.names=TRUE),
+        expected)
+    expect_identical(
+        make_DORsummary(DNAStringSet(x), disambiguate.allele.names=TRUE),
+        expected)
+})
+
 .check_object_mcols <- function(object)
 {
     object_mcols <- mcols(object)

TODO

@@ -41,12 +41,12 @@ test_that("combine_germline_dbs()", {
     db_name2 <- "_OGRDB.mouse.PWD_PhJ.IGH+IGK+IGL.202410"
 
     ## Combine 'db_name1' and 'db_name2'.
-    db_name12  <- "comb.OGRDB.human+mouse.IGH+IGK+IGL"
+    db_name12 <- "comb.OGRDB.human+mouse.IGH+IGK+IGL"
     combine_germline_dbs(db_name12, db_name1, db_name2, "_Hs", "_Mm")
     .check_combined_germline_db(db_name12, db_name1, db_name2, "_Hs", "_Mm")
 
     ## Combine 'db_name2' and 'db_name1'.
-    db_name21  <- "comb.OGRDB.mouse+human.IGH+IGK+IGL"
+    db_name21 <- "comb.OGRDB.mouse+human.IGH+IGK+IGL"
     combine_germline_dbs(db_name21, db_name2, db_name1, "_Mm", "_Hs")
     .check_combined_germline_db(db_name21, db_name2, db_name1, "_Mm", "_Hs")
 
@@ -93,5 +93,8 @@ test_that("combine_germline_dbs()", {
     expect_true(all(AIRR_df12$locus %in% c("IGK", "IGL")))
     expect_true(all(igblastr:::has_suffix(AIRR_df21$v_call, "_Hs")))
     expect_true(all(igblastr:::has_suffix(AIRR_df21$j_call, "_Hs")))
+
+    rm_germline_db(db_name12)
+    rm_germline_db(db_name21)
 })