R/union.R

In gson: Base Class and Methods for 'gson' Format

##' Combine gene set collections
##'
##' `gson_union()` merges `GSON` objects into one collection. Two objects are the
##' same source when they declare the same `gsname`, and -- once every object
##' declares one -- the same `gsidtype`; the gene sets of one source are joined by
##' ID, so a gene set that appears in two files of one collection becomes one set
##' holding every gene either file gives it. A gene set ID that two objects from
##' different sources both use is a conflict, because the merged membership table
##' would say that one ID names two sets of genes, and `.conflict` decides what
##' happens: `"error"` refuses, `"prefix"` qualifies every gene set ID with the
##' source it comes from (`KEGG:hsa00010`), and `"rename"` qualifies only the IDs
##' that actually collide, leaving the rest of the collection alone. Where two
##' objects name the same gene set or the same gene differently, the first object
##' that says something wins.
##'
##' The metadata is merged by what the merged object can still honestly say:
##' `species`, `keytype` and `schema_version` have to agree or the call is an
##' error, because a membership table holds genes of one organism identified by
##' one identifier type, while `gsname`, `gsidtype`, `version`, `accessed_date`
##' and `info` are joined with `" + "` when the sources name themselves
##' differently. `urlpattern` is kept only when every object declares the same
##' one and no gene set ID was renamed -- a pattern that resolves to a wrong page
##' for part of the collection is worse than no pattern, and `gson_url()` says so.
##'
##' @param x,y A `GSON` or `GSONList` object. Either may be omitted.
##' @param ... More `GSON` or `GSONList` objects.
##' @param .conflict What to do about a gene set ID used by two different
##'   sources: `"error"`, `"prefix"` or `"rename"`.
##'
##' @return A `GSON` object holding every gene set of every input.
##' @examples
##' wpfile <- system.file('extdata', "wikipathways-20220310-gmt-Homo_sapiens.gmt", package='gson')
##' wp <- read.gmt.wp(wpfile, output = "GSON")
##' gson_union(wp[1:2], wp[3:4])
##' @importFrom methods is
##' @export
gson_union <- function(x = NULL, y = NULL, ...,
                       .conflict = c("error", "prefix", "rename")) {
    .conflict <- match.arg(.conflict)
    objects <- gson_union_inputs(c(list(x, y), list(...)))
    labels <- vapply(seq_along(objects), function(i)
        gson_source_label(objects[[i]], i), character(1))
    by_domain <- all(vapply(objects, function(o)
        length(union_values(list(o@gsidtype))) > 0, logical(1)))
    identities <- vapply(objects, gson_source_identity, character(1),
                         by_domain = by_domain)
    memberships <- lapply(objects, function(o) usable_memberships(o@gsid2gene))
    ids <- lapply(memberships, function(g) unique(g$gsid))

    colliding <- colliding_ids(ids, identities)
    if (.conflict == "error" && length(colliding) > 0) {
        stop("gene set ID(s) used by more than one source: ",
             paste(colliding[seq_len(min(5, length(colliding)))], collapse = ", "),
             " -- merge with .conflict = \"prefix\" to qualify every ID with its ",
             "source, or .conflict = \"rename\" to qualify only these",
             call. = FALSE)
    }
    distinct_label <- unique(data.frame(label = labels, identity = identities,
                                        stringsAsFactors = FALSE))
    if (.conflict != "error" && anyDuplicated(distinct_label$label) > 0) {
        ambiguous <- distinct_label$label[duplicated(distinct_label$label)][1]
        stop("gson_union() cannot qualify gene set IDs with their source: two of ",
             "the collections are both named '", ambiguous, "' but hold different ",
             "gene set ID domains. Give them different gsname values before ",
             "merging.", call. = FALSE)
    }

    maps <- lapply(seq_along(objects), function(i)
        stats::setNames(rename_ids(ids[[i]], labels[i], colliding, .conflict),
                        ids[[i]]))
    origin <- rep(identities, lengths(ids))
    merged_ids <- unlist(maps, use.names = FALSE)
    ## the same gene set ID in two files of one source is one gene set, so only
    ## a pair of sources that end up sharing an ID cannot be merged
    per_source <- unique(data.frame(identity = origin, gsid = merged_ids,
                                    stringsAsFactors = FALSE))
    if (anyDuplicated(per_source$gsid) > 0) {
        stop("qualifying the gene set IDs with their source still leaves the ID '",
             per_source$gsid[duplicated(per_source$gsid)][1],
             "' meaning two different gene sets, so the objects cannot be merged ",
             "into one collection", call. = FALSE)
    }

    gsid2gene <- bind_tables(memberships, maps, c("gsid", "gene"), "gsid")
    if (is.null(gsid2gene)) {
        gsid2gene <- data.frame(gsid = character(0), gene = character(0))
    }
    gsid2gene <- gsid2gene[!duplicated(gsid2gene), , drop = FALSE]
    gsid2name <- bind_tables(lapply(objects, function(o) o@gsid2name), maps,
                             c("gsid", "name"), "gsid")
    if (!is.null(gsid2name)) {
        gsid2name <- gsid2name[!duplicated(gsid2name$gsid), , drop = FALSE]
    }
    gene2name <- bind_gene_names(objects)
    if (!is.null(gene2name)) {
        ## a name for a gene no set contains is a dangling reference
        gene2name <- gene2name[gene2name$gene %in% gsid2gene$gene, , drop = FALSE]
    }

    renamed <- any(vapply(seq_along(objects), function(i)
        any(maps[[i]] != ids[[i]]), logical(1)))
    gson(gsid2gene = gsid2gene,
         gsid2name = gsid2name,
         gene2name = gene2name,
         schema_version = union_scalar(
             lapply(objects, function(o) o@schema_version), "schema_version",
             paste0("the collections are written to different GSON file ",
                    "schemas")),
         species = union_scalar(
             lapply(objects, function(o) o@species), "species",
             "a gene set collection holds the genes of one species"),
         gsname = union_text(lapply(objects, function(o) o@gsname)),
         version = union_text(lapply(objects, function(o) o@version)),
         accessed_date = union_text(lapply(objects, function(o) o@accessed_date)),
         keytype = union_scalar(
             lapply(objects, function(o) o@keytype), "keytype",
             "a membership table holds the genes of one identifier type"),
         gsidtype = union_text(lapply(objects, function(o) o@gsidtype)),
         urlpattern = union_urlpattern(objects, renamed),
         info = union_text(lapply(objects, function(o) o@info)))
}

## accept both loose objects and GSONLists, so that a union of a list of
## collections is the same call as a union of the collections themselves
gson_union_inputs <- function(args) {
    out <- list()
    for (arg in args) {
        if (is.null(arg)) {
            next
        }
        if (is(arg, "GSONList")) {
            pieces <- unclass(arg)
        } else if (is(arg, "GSON")) {
            pieces <- list(arg)
        } else {
            stop("gson_union() takes GSON and GSONList objects, not ",
                 class(arg)[1], call. = FALSE)
        }
        for (piece in pieces) {
            if (!is(piece, "GSON")) {
                stop("gson_union() takes GSON and GSONList objects, and one of ",
                     "them holds a ", class(piece)[1], call. = FALSE)
            }
            out[[length(out) + 1]] <- piece
        }
    }
    if (length(out) == 0) {
        stop("gson_union() needs at least one GSON object", call. = FALSE)
    }
    out
}

## gsname first, because that is what a reader recognises; gsidtype is the next
## best thing to say where an ID came from, and a positional label is the last
gson_source_label <- function(x, i) {
    for (value in list(x@gsname, x@gsidtype)) {
        if (!is.null(value) && length(value) > 0 &&
            !is.na(value[1]) && nzchar(value[1])) {
            return(value[1])
        }
    }
    paste0("source", i)
}

## two files of one source hold the same gene sets; two ID domains behind one
## name do not. The domain only tells them apart when every collection declares
## one, because a file written before gsidtype existed is not a different source
## from one written after it.
gson_source_identity <- function(x, by_domain) {
    part <- function(value) {
        declared <- union_values(list(value))
        paste(declared, collapse = "\r")
    }
    if (by_domain) {
        return(paste0(part(x@gsname), "\r", part(x@gsidtype)))
    }
    part(x@gsname)
}

colliding_ids <- function(ids, identities) {
    key <- data.frame(gsid = unlist(ids, use.names = FALSE),
                      identity = rep(identities, lengths(ids)),
                      stringsAsFactors = FALSE)
    if (nrow(key) == 0) {
        return(character(0))
    }
    distinct <- vapply(split(key$identity, key$gsid, drop = TRUE),
                       function(v) length(unique(v)) > 1, logical(1))
    names(distinct)[distinct]
}

rename_ids <- function(id, label, colliding, .conflict) {
    hit <- rep(TRUE, length(id))
    if (.conflict == "rename") {
        hit <- id %in% colliding
    } else if (.conflict != "prefix") {
        return(id)
    }
    ## an ID that already carries its source as a prefix -- every GO term --
    ## does not need it twice
    hit <- hit & !startsWith(id, paste0(label, ":"))
    id[hit] <- paste0(label, ":", id[hit])
    id
}

## one table per pair of columns, every row carrying the ID its source ended up
## with; a gene set ID that no table row holds is left out of a name table that
## could not be looked up anyway
bind_tables <- function(tables, maps, cols, key) {
    holds_rows <- !vapply(tables, is.null, logical(1))
    tables <- tables[holds_rows]
    maps <- maps[holds_rows]
    if (length(tables) == 0) {
        return(NULL)
    }
    out <- do.call(rbind, lapply(seq_along(tables), function(i) {
        table <- tables[[i]]
        table[[key]] <- maps[[i]][table[[key]]]
        table
    }))
    if (is.null(out)) {
        return(NULL)
    }
    out <- out[!is.na(out[[key]]), , drop = FALSE]
    colnames(out) <- cols
    out
}

bind_gene_names <- function(objects) {
    gene2name <- lapply(objects, function(o) o@gene2name)
    if (all(vapply(gene2name, is.null, logical(1)))) {
        return(NULL)
    }
    gene2name <- gene2name[!vapply(gene2name, is.null, logical(1))]
    out <- do.call(rbind, gene2name)
    if (is.null(out)) {
        return(NULL)
    }
    colnames(out) <- c("gene", "name")
    out[!duplicated(out$gene), , drop = FALSE]
}

union_values <- function(values) {
    values <- values[!vapply(values, is.null, logical(1))]
    if (length(values) == 0) {
        return(character(0))
    }
    distinct <- unique(unlist(values, use.names = FALSE))
    distinct[!is.na(distinct) & nzchar(distinct)]
}

union_text <- function(values) {
    distinct <- union_values(values)
    if (length(distinct) == 0) {
        return(NULL)
    }
    paste(distinct, collapse = " + ")
}

## a field the merged object cannot hold two values of: refuse rather than pick
## one, because which one it keeps is a claim about the whole collection
union_scalar <- function(values, what, why) {
    distinct <- union_values(values)
    if (length(distinct) == 0) {
        return(NULL)
    }
    if (length(distinct) == 1) {
        return(distinct)
    }
    stop("the collections to merge declare different ", what, " values (",
         paste(distinct, collapse = ", "), "); ", why, call. = FALSE)
}

union_urlpattern <- function(objects, renamed) {
    patterns <- vapply(objects, function(o) {
        pattern <- gson_url_pattern(o)
        if (is.null(pattern)) NA_character_ else pattern
    }, character(1))
    if (anyNA(patterns) || length(unique(patterns)) != 1) {
        return(NULL)
    }
    if (renamed) {
        return(NULL)
    }
    patterns[[1]]
}