igblastr 1.3.1: Fix representation of missing values in custom annotations

We cannot use the string NA to represent missing values in the custom "ndm"
or ".aux" files that we create, because this crashes IgBLAST!
So now we use the following special values instead:
  o "N/A" for missing character values;
  o -1 for missing integer values.
See long comment at the top of R/igdata-IO.R for more information.

Author: hpages

DESCRIPTION

@@ -17,7 +17,7 @@ Description: The igblastr package provides functions to conveniently install
 biocViews: Immunology, Immunogenetics, ImmunoOncology, CellBiology
 URL: https://bioconductor.org/packages/igblastr
 BugReports: https://github.com/HyrienLab/igblastr/issues
-Version: 1.3.0
+Version: 1.3.1
 License: Artistic-2.0
 Encoding: UTF-8
 Authors@R: c(
@@ -42,7 +42,7 @@ Collate:
 	utils.R
 	internet-utils.R
 	long_to_wide_airr.R
-	read_broken_table.R
+	igdata-IO.R
 	translate_codons.R
 	allele2gene.R
 	parse_imgt_fasta_headers.R

R/OGRDB-utils.R

@@ -15,7 +15,7 @@ OGRDB_URL <- "https://ogrdb.airr-community.org/"
 encode_OGRDB_URL_component <- function(component)
 {
     stopifnot(isSingleString(component))
-    gsub("/", "%25252f", URLencode(component))
+    gsub("/", "%25252f", URLencode(component), fixed=TRUE)
 }
 
 .OGRDB_organism_url <- function(organism, for.download=FALSE)
@@ -280,7 +280,7 @@ download_OGRDB_germline_set <-
 .encode_OGRDB_path_component <- function(path_component)
 {
     stopifnot(isSingleNonWhiteString(path_component))
-    gsub("/", ".slash.", chartr(" ", "_", path_component))
+    gsub("/", ".slash.", chartr(" ", "_", path_component), fixed=TRUE)
 }
 
 .get_path_to_stored_OGRDB_germline_set <-

R/auxdata-IO.R

@@ -48,15 +48,12 @@ check_auxdata_col2class <- function(auxdata, what="'auxdata'")
 
 read_auxdata <- function(filepath)
 {
-    read_broken_table(filepath, AUXDATA_COL2CLASS)
+    read_igdata(filepath, AUXDATA_COL2CLASS)
 }
 
 write_auxdata <- function(auxdata, file="")
 {
     check_auxdata_col2class(auxdata)
-    header <- paste0("#", paste(colnames(auxdata), collapse=", "))
-    cat(header, "\n", sep="", file=file)
-    write.table(auxdata, file, append=TRUE, quote=FALSE,
-                sep="\t", row.names=FALSE, col.names=FALSE)
+    write_igdata(auxdata, file=file)
 }
 

R/clean_allele_set.R

@@ -317,7 +317,7 @@ checkarg_with.intdata <- function(with.intdata, gapped)
 
     ## Annotate the J alleles.
     auxdata <- compute_auxdata(allele_set, codon_starts=codon_starts,
-                                           no.warning=TRUE)
+                                           no.warnings=TRUE)
 
     ## Sanity checks.
     stopifnot(identical(auxdata[ , "allele_name"], allele_names))

R/compute_auxdata.R

@@ -177,16 +177,33 @@ J_alleles_with_missing_coding_frame_start <- function(auxdata)
     bad_alleles <- J_alleles_with_missing_coding_frame_start(auxdata)
     if (length(bad_alleles) != 0L) {
         in1string <- paste(bad_alleles, collapse=", ")
-        warning(wmsg("the \"coding frame start\" could not be determined ",
-                     "for J allele(s): ", in1string),
+        warning(wmsg("The \"coding frame start\" could not be determined ",
+                     "for J allele(s): ", in1string, "."),
                 "\n  ",
                 wmsg("--> coding_frame_start, cdr3_end, and extra_bps ",
                      "were set to NA for these alleles"))
     }
 }
 
+### Very rare but can happen (we've only seen this for rhesus monkey
+### allele IGHJ0-ZXTW*01 from OGRDB so far).
+warn_if_negative_cdr3_end <- function(auxdata, what)
+{
+    bad_idx <- which(auxdata[ , "cdr3_end"] < 0L)
+    if (length(bad_idx) != 0L) {
+        in1string <- paste(auxdata[bad_idx , "allele_name"], collapse=", ")
+        warning(wmsg("The ", what, " is negative for J allele(s): ",
+                     in1string, "."),
+                "\n  ",
+                wmsg("Note that you won't be able to save the returned ",
+                     "data.frame with write_auxdata() unless you replace ",
+                     "the negative \"cdr3_end\" value(s) with NA(s) or drop ",
+                     "the row(s) with a negative \"cdr3_end\" value."))
+    }
+}
+
 ### Returns a data.frame with 1 row per sequence in 'J_alleles'.
-compute_auxdata <- function(J_alleles, codon_starts=NULL, no.warning=FALSE)
+compute_auxdata <- function(J_alleles, codon_starts=NULL, no.warnings=FALSE)
 {
     if (!is(J_alleles, "DNAStringSet"))
         stop(wmsg("'J_alleles' must be DNAStringSet object"))
@@ -202,8 +219,8 @@ compute_auxdata <- function(J_alleles, codon_starts=NULL, no.warning=FALSE)
     if (!is.null(codon_starts))
         codon_starts <- .normarg_codon_starts(codon_starts, allele_names)
 
-    if (!isTRUEorFALSE(no.warning))
-        stop(wmsg("'no.warning' must be TRUE or FALSE"))
+    if (!isTRUEorFALSE(no.warnings))
+        stop(wmsg("'no.warnings' must be TRUE or FALSE"))
 
     JH_alleles <- J_alleles[allele_groups == "IGHJ"]
     JK_alleles <- J_alleles[allele_groups == "IGKJ"]
@@ -223,8 +240,11 @@ compute_auxdata <- function(J_alleles, codon_starts=NULL, no.warning=FALSE)
     ## Add "extra_bps" column.
     ans$extra_bps <- (width(J_alleles) - ans[ , "coding_frame_start"]) %% 3L
 
-    if (!no.warning)
+    if (!no.warnings) {
         .warn_if_missing_coding_frame_starts(ans)
+        what <- "computed CDR3 end position"
+        warn_if_negative_cdr3_end(ans, what)
+    }
 
     ans
 }

R/create_region_db.R

@@ -280,12 +280,11 @@ create_V_region_db <- function(fasta_files, destdir,
 ### create_J_region_db()
 ###
 
-### Not exported!
 ### Adds the auxdata to the db ONLY if its "coding_frame_start" column
 ### contains no NAs.
 ### Returns TRUE if that's actually the case (and therefore the auxdata
 ### was added), and FALSE otherwise.
-add_computed_auxdata_to_db <- function(auxdata, destdir, verbose=FALSE)
+.add_computed_auxdata_to_db <- function(auxdata, destdir, verbose=FALSE)
 {
     bad_alleles <- J_alleles_with_missing_coding_frame_start(auxdata)
     if (length(bad_alleles) != 0L) {
@@ -313,6 +312,11 @@ add_computed_auxdata_to_db <- function(auxdata, destdir, verbose=FALSE)
     auxdata_dir <- dirname(auxdata_path)
     stopifnot(!dir.exists(auxdata_dir))
     dir.create(auxdata_dir)
+
+    ## write_auxdata() will reject a data.frame with negative values in
+    ## the "cdr3_end" column, so we replace them with NAs.
+    bad_idx <- which(auxdata[ , "cdr3_end"] < 0L)
+    auxdata[bad_idx, "cdr3_end"] <- NA_integer_
     write_auxdata(auxdata, auxdata_path)
 
     if (verbose)
@@ -363,7 +367,7 @@ create_J_region_db <- function(fasta_files, destdir,
         auxdata <- mcols(allele_set, use.names=FALSE)
         stopifnot(identical(allele_names, auxdata[ , "allele_name"]))
         auxdata <- as.data.frame(auxdata)
-        add_computed_auxdata_to_db(auxdata, destdir, verbose=verbose)
+        .add_computed_auxdata_to_db(auxdata, destdir, verbose=verbose)
     }
 
     allele_set

R/download_OGRDB_germline_json.R

@@ -139,6 +139,7 @@ download_OGRDB_germline_json <- function(organism, germline_sets,
     setNames(as.character(allele_description[fields]), fields)
 }
 
+### Returns a character matrix with colnames on it.
 .make_matrix_from_allele_data_list <- function(data, colnames, what)
 {
     stopifnot(is.list(data), is.character(colnames))
@@ -163,6 +164,7 @@ download_OGRDB_germline_json <- function(organism, germline_sets,
 ### extract_intdata_from_ogrdb_json()
 ###
 
+### Returns a named character vector.
 .extract_IMGT_v_gene_delineation <-
     function(allele_description, extra_fields=NULL)
 {
@@ -229,7 +231,7 @@ extract_intdata_from_ogrdb_json <- function(json_path, extra_fields=NULL)
     df <- matrix2df(m[ , names(col2class)], col2class)
     df <- cbind(df, coding_frame_start=integer(nrow(df)))
     if (!is.null(extra_fields))
-        df <- cbind(df, m[ , extra_fields])
+        df <- cbind(df, m[ , extra_fields, drop=FALSE])
     df
 }
 
@@ -362,6 +364,8 @@ validate_OGRDB_intdata <- function(organism, germline_set, source_set=FALSE,
 ### first 6 amino acids (AEYFQH) are the last 6 amino acids of the CDR3,
 ### confirming that the 18th nucleotide in the coding sequence is where the
 ### CDR3 ends.
+###
+### Returns a named character vector.
 .extract_j_annotation <-
     function(allele_description, extra_fields=NULL)
 {
@@ -408,8 +412,10 @@ extract_auxdata_from_ogrdb_json <- function(json_path, extra_fields=NULL)
     what <- c("J allele descriptions found in JSON file: ", json_path)
     m <- .make_matrix_from_allele_data_list(data, expected_names, what)
     df <- matrix2df(m[ , names(AUXDATA_COL2CLASS)], AUXDATA_COL2CLASS)
+    what <- "CDR3 end position reported in the JSON file"
+    warn_if_negative_cdr3_end(df, what)
     if (!is.null(extra_fields))
-        df <- cbind(df, m[ , extra_fields])
+        df <- cbind(df, m[ , extra_fields, drop=FALSE])
     df
 }
 

R/igblastn.R

@@ -180,9 +180,18 @@ print.igblastn_raw_output <- function(x, ...) cat(x, sep="\n")
     stderr_lines <- readLines(stderr_file)
     if (length(stderr_lines) == 0L && status == 0L)
         return()
-    ## When the 'igblastn' executable crashes, we typically get a 0 status
-    ## and an empty out file.
-    if (status == 0L && isTRUE(out_is_empty))
+    ## When the 'igblastn' executable crashes, we typically get an empty
+    ## out file. But is an empty out file always indicating a crash? We're
+    ## not sure. So we add a few more conditions before concluding that we
+    ## got a crash:
+    ## - In our experience an empty out file and a 0 status seem to always
+    ##   be associated with a crash.
+    ## - An empty out file, a non-0 status, and a single stderr line
+    ##   (e.g. "Segmentation fault (core dumped)") also seem to always
+    ##   be associated with a crash.
+    crash <- isTRUE(out_is_empty) &&
+             (status == 0L || length(stderr_lines) == 1L)
+    if (crash)
         .stop_on_igblastn_exe_errors(stderr_lines, crash=TRUE)
     warn_msgs <- .parse_errors_or_warnings(stderr_lines, "warning:")
     for (msg in warn_msgs)