Add 'source_set' argument to download_OGRDB_germline_sequences()

Author: hpages

R/OGRDB-utils.R

@@ -215,6 +215,8 @@ download_OGRDB_germline_set <-
     vapply(query, .encode_OGRDB_URL_component, character(1))
 }
 
+
+### No longer used. Superseded by download_OGRDB_germline_set() above.
 ### Retrieves:
 ###     /germline/set/
 ###       {species}/{set_name}/{release_version}/{format}
@@ -314,255 +316,3 @@ download_OGRDB_germline_set <-
     ans  # named list
 }
 
-
-### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-### download_germline_sequences_from_OGRDB()
-###
-
-.infer_region_type_from_seqid <- function(seqid, locus)
-{
-    stopifnot(is.character(seqid), isSingleNonWhiteString(locus),
-              all(has_prefix(seqid, locus)))
-    pos <- nchar(locus) + 1L
-    region_type <- substr(seqid, pos, pos)
-    stopifnot(all(region_type %in% VDJ_REGION_TYPES))
-    region_type
-}
-
-.check_ungapped_seqs <- function(ungapped_seqs,
-                                 species, species_subgroup=NULL, set_name,
-                                 release_version="published",
-                                 extended=FALSE)
-{
-    format <- if (extended) "airr_ex" else "airr"
-    parsed_json <- .fetch_germline_set_from_OGRDB(species,
-                        species_subgroup=species_subgroup,
-                        set_name=set_name,
-                        release_version=release_version,
-                        format=format)
-    stop("paranoid.mode not ready yet")
-}
-
-### Downloads the germline sequences (ungapped) for the specified
-### species/species_subgroup/set_name/release_version.
-### Produces a subset of the following 7 FASTA files: IGHV.fasta,
-### IGHD.fasta, IGHJ.fasta, IGKV.fasta, IGKJ.fasta, IGLV.fasta, and
-### IGLJ.fasta. The exact subset produced depends on the species/set_name.
-### Returns the number of files produced.
-### NOTE: Superseded by download_OGRDB_germline_sequences()!
-download_germline_sequences_from_OGRDB <-
-    function(species, species_subgroup=NULL, set_name,
-             locus=c("IGH", "IGK", "IGL"),
-             release_version="published",
-             extended=FALSE, destdir=".",
-             paranoid.mode=FALSE)
-{
-    stopifnot(isSingleNonWhiteString(release_version),
-              isTRUEorFALSE(extended),
-              isSingleNonWhiteString(destdir),
-              isTRUEorFALSE(paranoid.mode))
-    locus <- match.arg(locus)
-
-    format <- if (extended) "ungapped_ex" else "ungapped"
-    fasta_lines <- .fetch_germline_set_from_OGRDB(species,
-                        species_subgroup=species_subgroup,
-                        set_name=set_name,
-                        release_version=release_version,
-                        format=format)
-    filepath <- tempfile(fileext=".fasta")
-    writeLines(fasta_lines, filepath)
-    ungapped_seqs <- readDNAStringSet(filepath)
-    seq_region_types <- .infer_region_type_from_seqid(names(ungapped_seqs),
-                                                      locus)
-
-    if (paranoid.mode)
-        .check_ungapped_seqs(ungapped_seqs,
-                             species, species_subgroup=species_subgroup,
-                             set_name=set_name,
-                             release_version=release_version,
-                             extended=extended)
-
-    from <- paste0("germline set ", set_name, " (", species, ")")
-    file_count <- 0L
-    for (region_type in VDJ_REGION_TYPES) {
-        selected_seqs <- ungapped_seqs[seq_region_types == region_type]
-        if (length(selected_seqs) == 0L)
-            next
-        filename <- paste0(locus, region_type, ".fasta")
-        destfile <- file.path(destdir, filename)
-        message("Write ", length(selected_seqs), " ", region_type, " regions ",
-                "from ", from, " to ", filename, " ... ", appendLF=FALSE)
-        if (file.exists(destfile))
-            stop(wmsg(filename, " file already exists in ", destdir, "/"))
-        writeXStringSet(selected_seqs, destfile)
-        message("ok")
-        file_count <- file_count + 1L
-    }
-    file_count
-}
-
-
-### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-### download_human_germline_sequences_from_OGRDB()
-###
-
-### Downloads the germline sequences (ungapped) for all the Human germline
-### sets. See list of Human germline sets here:
-###   https://ogrdb.airr-community.org/germline_sets/Homo%20sapiens
-### Returns the number of files produced (should be 7).
-###
-### NOTE: No longer needed! Use download_OGRDB_germline_sequences() instead.
-download_human_germline_sequences_from_OGRDB <-
-    function(release_version="published",
-             extended=FALSE, destdir=".",
-             paranoid.mode=FALSE)
-{
-    stopifnot(isSingleNonWhiteString(release_version),
-              isTRUEorFALSE(extended),
-              isSingleNonWhiteString(destdir),
-              isTRUEorFALSE(paranoid.mode))
-
-    ## Produces files: IGHV.fasta, IGHD.fasta, IGHJ.fasta
-    file_count <-
-        download_germline_sequences_from_OGRDB("Human",
-                    set_name="IGH_VDJ", locus="IGH",
-                    release_version=release_version, extended=extended,
-                    destdir=destdir, paranoid.mode=paranoid.mode)
-    ## Produces files: IGKV.fasta, IGKJ.fasta
-    file_count <- file_count +
-        download_germline_sequences_from_OGRDB("Human",
-                    set_name="IGKappa_VJ", locus="IGK",
-                    release_version=release_version, extended=extended,
-                    destdir=destdir, paranoid.mode=paranoid.mode)
-    ## Produces files: IGLV.fasta, IGLJ.fasta
-    file_count <- file_count +
-        download_germline_sequences_from_OGRDB("Human",
-                    set_name="IGLambda_VJ", locus="IGL",
-                    release_version=release_version, extended=extended,
-                    destdir=destdir, paranoid.mode=paranoid.mode)
-    file_count
-}
-
-
-### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
-### download_mouse_germline_sequences_from_OGRDB()
-###
-
-### 'set_names' is expected to be a vector of Set Names from
-### the Germline Sets table displayed at
-###   https://ogrdb.airr-community.org/germline_sets/Mus%20musculus
-### except special Set Names "IGKJ (all strains)" and "IGLJ (all strains)".
-### Assumes that all the Set Names in 'set_names' are of the
-### form "<Species subgroup> <group>", where <group> is
-### one of IGH, IGKV, an IGLV.
-### Returns a 4-column matrix with colum names: set_name, species_subgroup,
-### group, and locus.
-.set_names_as_matrix <- function(set_names)
-{
-    if (!is.character(set_names) || anyNA(set_names))
-        stop(wmsg("'set_names' must be a character vector with no NAs"))
-    set_names <- trimws2(set_names)
-    if (anyDuplicated(set_names))
-        stop(wmsg("'set_names' cannot contain duplicates"))
-    FORBIDDEN_SET_NAMES <- c("IGKJ (all strains)", "IGLJ (all strains)")
-    if (any(tolower(set_names) %in% tolower(FORBIDDEN_SET_NAMES))) {
-        in1string <- paste0('"', FORBIDDEN_SET_NAMES, '"', collapse=" or ")
-        stop(wmsg("'set_names' cannot contain ", in1string))
-    }
-    split_set_names <- strsplit(set_names, split=" ", fixed=TRUE)
-    if (!all(lengths(split_set_names) == 2L))
-        stop(wmsg("each Set Name in 'set_names' must contain exactly 1 space"))
-    data <- unlist(split_set_names)
-    if (is.null(data))
-        data <- character(0)
-    m <- matrix(data, ncol=2L, byrow=TRUE,
-                dimnames=list(NULL, c("species_subgroup", "group")))
-    loci <- substr(m[ , "group"], 1L, 3L)
-    bad_idx <- which(!(loci %in% c("IGH", "IGK", "IGL")))
-    if (length(bad_idx) != 0L) {
-        in1string <- paste(set_names[bad_idx], collapse=", ")
-        stop(wmsg("Bad Set Name(s): ", in1string, "."),
-             "\n  ",
-             wmsg("The first 3 letters of the 2nd part of the name ",
-                  "is not a valid locus name."))
-    }
-    set_name <- matrix(set_names, dimnames=list(NULL, "set_name"))
-    locus    <- matrix(loci, dimnames=list(NULL, "locus"))
-    ans <- cbind(set_name, m, locus)
-    species_subgroup <- unique(ans[ , "species_subgroup"])
-    if (length(species_subgroup) > 1L)
-        warning(wmsg("the supplied Set Names are from ",
-                     "more than one Species subgroup"))
-    ans
-}
-
-### Downloads the germline sequences (ungapped) for the specified Mouse
-### germline sets. See list of Mouse germline sets here:
-###   https://ogrdb.airr-community.org/germline_sets/Mus%20musculus
-### Returns the number of files produced (will be typically between 4 and 7).
-###
-### To download all the germline sequences for Mouse strain A/J:
-###
-###   set_names <- c("A/J IGKV", "A/J IGLV")
-###   download_mouse_germline_sequences_from_OGRDB(set_names)
-###
-### --> produces 4 FASTA files (partial germline db).
-###
-### To download all the germline sequences for Mouse strain C57BL/6:
-###
-###   download_mouse_germline_sequences_from_OGRDB("C57BL/6 IGH")
-###
-### --> produces 5 FASTA files (partial germline db).
-###
-### To download all the germline sequences for Mouse strain C57BL/6J:
-###
-###   set_names <- c("C57BL/6J IGKV", "C57BL/6J IGLV")
-###   download_mouse_germline_sequences_from_OGRDB(set_names)
-###
-### --> produces 4 FASTA files (partial germline db).
-###
-### To download all the germline sequences for Mouse strain CAST/EiJ:
-###
-###   set_names <- c("CAST/EiJ IGH", "CAST/EiJ IGKV", "CAST/EiJ IGLV")
-###   download_mouse_germline_sequences_from_OGRDB(set_names)
-###
-### --> produces 7 FASTA files (full germline db).
-###
-### NOTE: No longer needed! Use download_OGRDB_germline_sequences() instead.
-download_mouse_germline_sequences_from_OGRDB <-
-    function(set_names,
-             release_version="published",
-             extended=FALSE, destdir=".",
-             paranoid.mode=FALSE)
-{
-    stopifnot(isSingleNonWhiteString(release_version),
-              isTRUEorFALSE(extended),
-              isSingleNonWhiteString(destdir),
-              isTRUEorFALSE(paranoid.mode))
-
-    m <- .set_names_as_matrix(set_names)
-    set_names         <- m[ , "set_name"]
-    species_subgroups <- m[ , "species_subgroup"]
-    loci              <- m[ , "locus"]
-    file_count <- 0L
-    for (i in seq_len(nrow(m))) {
-        file_count <- file_count +
-            download_germline_sequences_from_OGRDB("Mouse",
-                        species_subgroup=species_subgroups[[i]],
-                        set_name=set_names[[i]], locus=loci[[i]],
-                        release_version=release_version, extended=extended,
-                        destdir=destdir, paranoid.mode=paranoid.mode)
-    }
-    file_count <- file_count +
-        download_germline_sequences_from_OGRDB("Mouse",
-                    set_name="IGKJ (all strains)", locus="IGK",
-                    release_version=release_version, extended=extended,
-                    destdir=destdir, paranoid.mode=paranoid.mode)
-    file_count <- file_count +
-        download_germline_sequences_from_OGRDB("Mouse",
-                    set_name="IGLJ (all strains)", locus="IGL",
-                    release_version=release_version, extended=extended,
-                    destdir=destdir, paranoid.mode=paranoid.mode)
-    file_count
-}
-

R/download_OGRDB_germline_sequences.R

@@ -159,10 +159,9 @@
 ### Returns the list of FASTA files that were produced in an invisible
 ### character vector that carries the names of the corresponding germline
 ### sets.
-### TODO: Add 'source_set' argument like in download_V_ndm_data_from_OGRDB()
-### below.
 download_OGRDB_germline_sequences <- function(organism="Homo sapiens",
                                               germline_sets, gapped=TRUE,
+                                              source_set=FALSE,
                                               destdir=".", overwrite=FALSE,
                                               recache=FALSE, ...)
 {
@@ -172,6 +171,8 @@ download_OGRDB_germline_sequences <- function(organism="Homo sapiens",
     set_loci <- .infer_loci_from_OGRDB_set_names(set_names)
     if (!isTRUEorFALSE(gapped))
         stop(wmsg("'gapped' must be TRUE or FALSE"))
+    if (!isTRUEorFALSE(source_set))
+        stop(wmsg("'source_set' must be TRUE or FALSE"))
     if (!isSingleNonWhiteString(destdir))
         stop(wmsg("'destdir' must be a single (non-empty) string"))
     if (!dir.exists(destdir)) {
@@ -198,7 +199,8 @@ download_OGRDB_germline_sequences <- function(organism="Homo sapiens",
             ## already there.
             local_file <- .download_OGRDB_germline_set_to_OGRDB_store(organism,
                                           set_name, set_version,
-                                          format=format, recache=recache, ...)
+                                          format=format, source_set=source_set,
+                                          recache=recache, ...)
             .split_OGRDB_fasta_file(local_file, set_name, set_locus,
                                     destdir=tmp_destdir)
         }

TODO

@@ -197,30 +197,15 @@ o Great resource for how to use AIRR Community Reference germline sets with
   In particular, the author seems to be using an OGRDB REST API version 2:
     https://ogrdb.airr-community.org/api_v2
   but where is this API documented?
-  All the download utilities implemented in igblastr/R/OGRDB-utils.R use
-  the OGRDB REST API at
+  Note that internal utility igblastr:::.fetch_germline_set_from_OGRDB()
+  (implemented in R/OGRDB-utils.R) uses the OGRDB REST API at
     https://ogrdb.airr-community.org/api
-  which is poorly documented and is somewhat confusing (see below).
-
-o Investigate the following mysteries about the germline sets provided
-  by AIRR/OGRDB:
-
-  1. The OGRDB API at https://ogrdb.airr-community.org/api/ allows downloading
-     the germline sequences in 2 formats: ungapped or ungapped_ex.
-     Which format is appropriate to use with IgBLAST?
-     Note that downloading germline sets directly by clicking on
-     the "FASTA Ungapped" links here
-       https://ogrdb.airr-community.org/germline_sets/Homo%20sapiens
-     or
-       https://ogrdb.airr-community.org/germline_sets/Mus%20musculus
-     seems to retrieve the "ungapped_ex" sequences for Human and the "ungapped"
-     sequences for Mouse. Confusing!
-
-  2. For some Mouse strains, OGRDB seems to provide germline sequences
-     only for a limited number of loci/groups. For example for strain A/J,
-     only sequences from the light chain (i.e. groups IGKV, IGKJ, IGLV,
-     and IGLJ) seem to be available.
-     See https://ogrdb.airr-community.org/germline_sets/Mus%20musculus
+  which is poorly documented and is somewhat confusing.
+
+o Why are some loci/groups missing for some mouse strains at OGRDB?
+  For example for strain A/J, only sequences from the light chain (i.e.
+  groups IGKV, IGKJ, IGLV, and IGLJ) seem to be available.
+  See https://ogrdb.airr-community.org/germline_sets/Mus%20musculus
 
 o Implement install_OGRDB_germline_db(). Will download the germline sequences
   from https://ogrdb.airr-community.org/ (link provided by Kellie).

inst/extdata/germline_sets/AIRR/human/202309/ref/README.md

@@ -25,16 +25,8 @@ However, the FASTA files in this folder were obtained programmatically
 by running the following code in the folder on Oct 20, 2025:
 ```r
 library(igblastr)
-download_germline_sequences_from_OGRDB <- igblastr:::download_germline_sequences_from_OGRDB
-
-Hfile_count <- download_germline_sequences_from_OGRDB("Human", set_name="IGH_VDJ", locus="IGH", release_version="7", extended=TRUE)
-stopifnot(identical(Hfile_count, 3L))
-
-Kfile_count <- download_germline_sequences_from_OGRDB("Human", set_name="IGKappa_VJ", locus="IGK", release_version="2", extended=TRUE)
-stopifnot(identical(Kfile_count, 2L))
-
-Jfile_count <- download_germline_sequences_from_OGRDB("Human", set_name="IGLambda_VJ", locus="IGL", release_version="1", extended=TRUE)
-stopifnot(identical(Jfile_count, 2L))
+germline_sets <- c(IGH_VDJ=7, IGKappa_VJ=2, IGLambda_VJ=1)
+download_OGRDB_germline_sequences("Homo sapiens", germline_sets, gapped=FALSE)
 ```
 
 Files `*.ndm.imgt` were obtained programmatically by running the

inst/extdata/germline_sets/AIRR/human/202309/src/README.md

@@ -22,15 +22,7 @@ However, the FASTA files in this folder were obtained programmatically
 by running the following code in the folder on Oct 20, 2025:
 ```r
 library(igblastr)
-download_germline_sequences_from_OGRDB <- igblastr:::download_germline_sequences_from_OGRDB
-
-Hfile_count <- download_germline_sequences_from_OGRDB("Human", set_name="IGH_VDJ", locus="IGH", release_version="7")
-stopifnot(identical(Hfile_count, 3L))
-
-Kfile_count <- download_germline_sequences_from_OGRDB("Human", set_name="IGKappa_VJ", locus="IGK", release_version="2")
-stopifnot(identical(Kfile_count, 2L))
-
-Jfile_count <- download_germline_sequences_from_OGRDB("Human", set_name="IGLambda_VJ", locus="IGL", release_version="1")
-stopifnot(identical(Jfile_count, 2L))
+germline_sets <- c(IGH_VDJ=7, IGKappa_VJ=2, IGLambda_VJ=1)
+download_OGRDB_germline_sequences("Homo sapiens", germline_sets, gapped=FALSE, source_set=TRUE)
 ```