R/LOD_POM.R

In OncoSimulR: Forward Genetic Simulation of Cancer Progression with Epistasis

## Copyright 2016, 2017 Ramon Diaz-Uriarte

## This program is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.

## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.

## You should have received a copy of the GNU General Public License
## along with this program.  If not, see <http://www.gnu.org/licenses/>.


## Note that, in contrast to POM, LOD is not well defined if the population
## becomes extinct.

## Functions to obtain LOD and POM similar to Szendro et al., 2014, PNAS.
genot_max <- function(x) {
    x$GenotypesLabels[which.max(x$pops.by.time[nrow(x$pops.by.time), -1])]
}



## Filter the PhylogDF so we obtain LOD, sensu stricto.

## No longer necessary with LOD from C++ removing duplicates
## ## Now this is coming from LOD_DF, which already only has
## ## implicitly pop_size_child == 0
## Only used in one function use for testing
filter_phylog_df_LOD <- function(y) {
    keep <- !rev(duplicated(rev(y$child)))
    return(y[keep, ])
}



## ## Filter the PhylogDF so we obtain LOD, sensu stricto.
## ## For the old version
## filter_phylog_df_LOD_with_n <- function(y) {
##     y <- y[y$pop_size_child == 0, , drop = FALSE]
##     keep <- !rev(duplicated(rev(y$child)))
##     return(y[keep, ])
## }


## from phylogClone, key parts for the LOD strict structure
phcl_from_lod <- function(df, x) {
    ## no longer necessary
    ## ## the !keepEvents. Which I move here to speed things up.
    ## df <- df[!duplicated(df[, c(1, 2)]), , drop = FALSE]
    
    tG <- unique(c(as.character(df[, 1]), as.character(df[, 2])))
    ## ## Do as in phylogClone. So that we have the same nodes
    ## ## in LOD all and not LOD all?
    ## z <- which_N_at_T(x, N = 1, t = "last")
    ## tG <- x$GenotypesLabels[z]
    
    ## ## FIXME: aren't these two warnings redundant or aliased?
    ## ## yes, partlt
    ##  I think this can never happen now
    ## if ((length(tG) == 1) && (tG == "")) {
    ##     warning("There never was a descendant of WT")
    ## }
    if (nrow(df) == 0) {
        warning("LOD structure has 0 rows: no descendants of initMutant ever appeared. ")
        return(NA)
    }
    g <- igraph::graph.data.frame(df[, c(1, 2)])
    nodesInP <- unique(unlist(igraph::neighborhood(g, order = 1e+09, 
                                                   nodes = tG, mode = "in")))
    allLabels <- unique(as.character(unlist(df[, c(1, 2)])))
    nodesRm <- setdiff(allLabels, V(g)$name[nodesInP])
    g <- igraph::delete.vertices(g, nodesRm)
    tmp <- list(graph = g, df = df)
    class(tmp) <- c(class(tmp), "phylogClone")
    return(tmp)
}

LOD.internal <- function(x) {
    if(is.null(x$pops.by.time)) {
        warning("Missing needed components. This might be a failed simulation.",
                " Returning NA.")
        return(list(all_paths = NA, lod_single = NA))
    }
    if (!inherits(x, "oncosimul2")) 
        stop("LOD information is only stored with v >= 2")
    ## y <- filter_phylog_df_LOD(x$other$LOD_DF)
    y <- x$other$LOD_DF
    pc <- phcl_from_lod(y)
    ## need eval for oncoSimulPop calls and for LOD_as_path
    initMutant <- x$InitMutant

    ## No descendants means that: never a descendant.
    ## Note the same that the init mutant be the final state.
    if((length(pc) == 1) && (is.na(pc))) {
        lod <- "No_descendants"
        ## bail out here. We do not need the rest.
        if(initMutant != "")
            attributes(lod)$initMutant <- initMutant
        return(lod)
    }
    
    pcg <- pc$graph
    end <- genot_max(x)
    
    if(end == initMutant) {
        if(initMutant == "") {
            stinitm <- "WT"
        } else {
            stinitm <- paste0("initMutant(", initMutant, ")")
        }
        lod <- paste0(stinitm, "_is_end")
    } else {
        all_paths <- igraph::all_simple_paths(pcg, from = initMutant, to = end,
                                              mode = "out")
        if(length(all_paths) > 1)
            stop("length(all_paths) > 1???")
        lod <- igraph::as_ids(all_paths[[1]])
    }
    if(initMutant != "")
        attributes(lod)$initMutant <- initMutant
    return(lod)
}


POM.internal <- function(x) {
    if(is.null(x$pops.by.time)) {
        warning("Missing needed components. This might be a failed simulation.",
                " Returning NA.")
        return(NA)
    } else {
        x$other$POM
    }
}



diversityLOD <- function(llod) {
    ## nn <- names(llod[[1]])
    nn <- llod[[1]]
    if( is_null_na(nn) ||
        !(is.list(llod)))
        stop("Object must be a list of LODs")
    pathstr <- unlist(lapply(llod, function(x) paste(x,
                                                     collapse = "_")))
    shannonI(table(pathstr))
}

## diversityLOD <- function(llod) {
##     nn <- names(llod[[1]])
##     if( is_null_na(nn) ||
##         !(is.list(llod)))
##         stop("Object must be a list of LODs")
##     pathstr <- unlist(lapply(llod, function(x) paste(names(x),
##                                                      collapse = "_")))
##     shannonI(table(pathstr))
## }

LOD_as_path <- function(llod) {
    path_l <- function(u) {
        if(length(u) == 1) {
            if(is.null(attributes(u)$initMutant))
                initMutant <- ""
            else
                initMutant <- attributes(u)$initMutant
            if(initMutant == "") initMutant <- "WT"
            if(grepl("_is_end", u))
                return(initMutant)
            if(u == "No_descendants")
                return(initMutant)
        } else {
            ## Deal with "" meaning WT
            ## the_names <- names(u)
            the_names <- u
            the_names_wt <- which(the_names == "")
            
            if(length(the_names_wt)) {
                if(length(the_names_wt) > 1) stop("more than 1 WT?!")
                if(the_names_wt > 1) stop("WT in position not 1?!")
                the_names[the_names_wt] <- "WT"
            }
            return(paste(the_names, collapse = " -> ")) 
        }
    }
    if(!is.list(llod))
        pathstr <- path_l(llod)
    else
        pathstr <- unlist(lapply(llod, path_l))
    return(pathstr)
}
## We would just need a LOD_as_DAG

diversityPOM <- function(lpom) {
    if(!inherits(lpom, "list"))
        stop("Object must be a list of POMs")
    pomstr <- unlist(lapply(lpom, function(x) paste(x, collapse = "_")))
    shannonI(table(pomstr))
}

## a legacy
diversity_POM <- diversityPOM
diversity_LOD <- diversityLOD


POM <- function(x) {
    UseMethod("POM", x)
}


LOD <- function(x) {
    UseMethod("LOD", x)
}

POM.oncosimul2 <- function(x) return(POM.internal(x))
LOD.oncosimul2 <- function(x) return(LOD.internal(x))

POM.oncosimulpop <- function(x) return(lapply(x, POM.internal))
LOD.oncosimulpop <- function(x) return(lapply(x, LOD.internal))




## POM.oncosimul2 <- function(x) {
##     out <- POM.internal(x)
##     class(out) <- c(class(out), "oncosimul_pom")
##     return(out)
## }

## LOD.oncosimul2 <- function(x) {
##     out <- LOD.internal(x)
##     class(out) <- c(class(out), "oncosimul_lod")
##     return(out)
## }


## POM.oncosimulpop <- function(x) {
##     out <- lapply(x, POM.internal)
##     class(out) <- c(class(out), "oncosimul_pom_list")
##     return(out)
## }

## LOD.oncosimulpop <- function(x) {
##     out <- lapply(x, LOD.internal)
##     class(out) <- c(class(out), "oncosimul_lod_list")
##     return(out)
## }


## summary.oncosimul_lod_list <- function(x) {
##     cat("List of ", length(x), " simulations\n.")
##     cat("Shannon's diversity (entropy) = ", diversityLOD(x), "\n")
## }

## summary.oncosimul_pom_list <- function(x) {
##     cat("List of ", length(x), " simulations\n.")
##     cat("Shannon's diversity (entropy) = ", diversityPOM(x), "\n")
## }



POM_pre_2.9.2 <- function(x) {
    if(is.null(x$pops.by.time)) {
        warning("Missing needed components. This might be a failed simulation.",
                " Returning NA.")
        return(NA)
    }
    x$GenotypesLabels[rle(apply(x$pops.by.time[, -1, drop = FALSE],
                                1, which.max))$values]
}



LOD.internal_pre_2.9.2 <- function(x, strict) {
    ## Not identical to LOD of Szendro because:
    
    ##  a) I think we want all paths, not just their single LOD, which I
    ##  think they might use out of convenience.

    ##  b) For me it is a mess, a complicated mess, to use their LOD as
    ##  they define it and there are many ambiguities in how to define it
    ##  in continuous time.

    ## This also means that single simulation might yield multiple LODs

    ## keepEvents is FALSE to make this object as small as possible.
    if(is.null(x$pops.by.time)) {
        warning("Missing needed components. This might be a failed simulation.",
                " Returning NA.")
        return(list(all_paths = NA, lod_single = NA))
    }
    if(strict) {
        if (!inherits(x, "oncosimul2")) 
            stop("LOD information is only stored with v >= 2")
        y <- filter_phylog_df_LOD(x$other$LOD_DF)
        pc <- phcl_from_lod(y)
    } else {
        pc <- phylogClone(x, keepEvents = FALSE)
    }
    ## need eval for oncoSimulPop calls and for LOD_as_path
    initMutant <- x$InitMutant
    
    if((length(pc) == 1) && (is.na(pc))) {
        lodlist <- list(all_paths = NA,
                        lod_single = "No_descendants")
        ## bail out here. We do not need the rest.
        if(initMutant != "")
            attributes(lodlist)$initMutant <- initMutant
        return(lodlist)
    }
    
    pcg <- pc$graph
    end <- genot_max(x)
    
    ## if(!is.null(eval(attributes(x)$call$initMutant))) {
    ##     initMutant <- eval(attributes(s7)$call$initMutant)
    ## } else {
    ##     initMutant <- ""
    ## }
    ## browser()
    if(end == initMutant) {
        if(initMutant == "") {
            stinitm <- "WT"
        } else {
            stinitm <- paste0("initMutant(", initMutant, ")")
        }
        lod_single <- paste0(stinitm, "_is_end")
        all_paths <- list(lod_single)
    } else {
        all_paths <- igraph::all_simple_paths(pcg, from = initMutant, to = end,
                                              mode = "out")
       
        if(!strict) {
            ## the next is partially redundant
            ## graph_to_end <- igraph::make_ego_graph(pcg, order = 1e9, nodes = end,
            ##                                        mode = "in")
            ## if(length(graph_to_end) != 1) stop("length(graph_to_end) > 1")
            ## I am not sure if I should keep the last one. Redundant
            
            ## This gives a single path and it is the first entry into each
            ## destination. But we do not check there is no extinction afterwards.
            ## The closest to the single Szendro LOD
            singlep <- pc$df
            singlep[, 1] <- as.character(singlep[, 1])
            singlep[, 2] <- as.character(singlep[, 2])
            singlep <- singlep[ do.call(order, singlep[, c(2, 3)]), ]
            singlep <- singlep[!duplicated(singlep[, 2]), ]
            gsingle <- igraph::graph_from_data_frame(singlep)
            lod_single <- igraph::all_simple_paths(gsingle, from = initMutant,
                                                   to = end, mode = "out")
            if(length(lod_single) != 1) stop("lod_single != 1")
        }
    }
    if(strict) {
        if(length(all_paths) > 1)
            stop("length(all_paths) > 1???")
        lodlist <- list(all_paths = NA,
                    lod_single = all_paths[[1]])
    } else {
        lodlist <- list(all_paths = all_paths,
                    lod_single = lod_single[[1]])
    }
    if(initMutant != "")
        attributes(lodlist)$initMutant <- initMutant
    return(lodlist)
}



## LOD_as_path_pre_2.9.2 <- function(llod) {
##     path_l <- function(u) {
##         if(length(u$lod_single) == 1) {
##             if(is.null(attributes(u)$initMutant))
##                 initMutant <- ""
##             else
##                 initMutant <- attributes(u)$initMutant
##             if(initMutant == "") initMutant <- "WT"
##             if(grepl("_is_end", u$lod_single))
##                 return(initMutant)
##             if(u$lod_single == "No_descendants")
##                 return(initMutant)
##         } else {
##             ## Deal with "" meaning WT
##             the_names <- names(u$lod_single)
##             the_names_wt <- which(the_names == "")
            
##             if(length(the_names_wt)) {
##                 if(length(the_names_wt) > 1) stop("more than 1 WT?!")
##                 if(the_names_wt > 1) stop("WT in position not 1?!")
##                 the_names[the_names_wt] <- "WT"
##             }
##             return(paste(the_names, collapse = " -> ")) 
##             ## return(paste0("WT", paste(names(u$lod_single),
##             ##                           collapse = " -> ")) )
##         }
##     }
##     if(identical(names(llod), c("all_paths", "lod_single")))
##         pathstr <- path_l(llod)
##     else {
##         ## should be a list
##         pathstr <- unlist(lapply(llod, path_l))
##     }
##     return(pathstr)
##     ## pathstr <- unlist(lapply(llod, function(x) paste(names(x$lod_single),
##     ##                                                  collapse = " -> ")))
##     ## return(paste0("WT", pathstr))
## }



LOD.oncosimul2_pre_2.9.2 <- function(x, strict = TRUE)
    return(LOD.internal_pre_2.9.2(x, strict))

## LOD.oncosimulpop_pre_2.9.2 <- function(x, strict = TRUE)
##     return(lapply(x, LOD.internal_pre_2.9.2, strict))






## Note for self: we could get all the LODs per simulation in the strict
## sense of those never becoming extinct if we subset the phylogClone
## object to children in which if we arrive at the children at any two
## times t and t+k, we retain only rows where any time > t is such that
## the popsize is > 0. But this is not worth it now.