REMOVE OLD INMF IMPLEMENTATION

CMakeLists.txt

@@ -13,7 +13,7 @@ set(CMAKE_CXX_FLAGS_RELEASE "")
 set(CMAKE_VERBOSE_MAKEFILE TRUE)
 
 set(srclist ModularityOptimizer.cpp RModularityOptimizer.cpp RcppExports.cpp
-data_processing.cpp fast_wilcox.cpp feature_mat.cpp quantile_norm.cpp run_nmf.cpp snn.cpp)
+data_processing.cpp fast_wilcox.cpp feature_mat.cpp quantile_norm.cpp snn.cpp)
 list(TRANSFORM srclist PREPEND "${PROJECT_SOURCE_DIR}/src/")
 add_library(rliger2 SHARED ${srclist})
 set_target_properties(rliger2 PROPERTIES PREFIX "")

DESCRIPTION

@@ -27,7 +27,7 @@ Author: Joshua Welch [aut],
         Robert Lee [ctb],
         Yichen Wang [aut, cre],
         Andrew Robbins [ctb]
-Maintainer: Chao Gao <gchao@umich.edu>
+Maintainer: Yichen Wang <wayichen@umich.edu>
 BugReports: https://github.com/welch-lab/liger/issues
 URL: https://github.com/welch-lab/liger
 License: GPL-3 | file LICENSE

NAMESPACE

@@ -11,9 +11,12 @@ S3method(as.ligerDataset,matrixLike)
 S3method(c,liger)
 S3method(cbind,ligerDataset)
 S3method(fortify,liger)
-S3method(runBPPINMF,Seurat)
-S3method(runBPPINMF,liger)
-S3method(runBPPINMF,list)
+S3method(runINMF,Seurat)
+S3method(runINMF,liger)
+S3method(runINMF,list)
+S3method(runOnlineINMF,Seurat)
+S3method(runOnlineINMF,liger)
+S3method(runOnlineINMF,list)
 S3method(runUINMF,Seurat)
 S3method(runUINMF,liger)
 S3method(runUINMF,list)
@@ -39,6 +42,7 @@ export(calcDatasetSpecificity)
 export(calcNormLoadings)
 export(calcPurity)
 export(cellMeta)
+export(closeAllH5)
 export(commandDiff)
 export(commands)
 export(convertOldLiger)
@@ -73,7 +77,6 @@ export(normPeak)
 export(normalize)
 export(normalizePeak)
 export(online_iNMF)
-export(optimizeALS)
 export(optimizeNewData)
 export(optimizeNewK)
 export(optimizeNewLambda)
@@ -117,14 +120,15 @@ export(removeMissing)
 export(restoreH5Liger)
 export(retrieveCellFeature)
 export(reverseMethData)
-export(runBPPINMF)
 export(runCluster)
 export(runDoubletFinder)
 export(runGOEnrich)
 export(runGSEA)
 export(runGeneralQC)
 export(runINMF)
+export(runINMF_R)
 export(runMarkerDEG)
+export(runOnlineINMF)
 export(runPairwiseDEG)
 export(runPseudoBulkDEG)
 export(runTSNE)
@@ -180,12 +184,11 @@ exportMethods(length)
 exportMethods(names)
 exportMethods(normData)
 exportMethods(normPeak)
-exportMethods(optimizeALS)
 exportMethods(quantileNorm)
 exportMethods(quantile_norm)
 exportMethods(rawData)
 exportMethods(rawPeak)
-exportMethods(runINMF)
+exportMethods(runINMF_R)
 exportMethods(scaleData)
 exportMethods(scaleUnsharedData)
 exportMethods(show)

R/RcppExports.R

@@ -115,10 +115,6 @@ max_factor <- function(H, dims_use, center_cols = FALSE) {
     .Call(`_rliger2_max_factor`, H, dims_use, center_cols)
 }
 
-solveNNLS <- function(C, B) {
-    .Call(`_rliger2_solveNNLS`, C, B)
-}
-
 ComputeSNN <- function(nn_ranked, prune) {
     .Call(`_rliger2_ComputeSNN`, nn_ranked, prune)
 }

R/h5Utility.R

@@ -195,11 +195,14 @@ restoreH5Liger <- function(object, filePath = NULL) {
         }
         h5.meta <- h5fileInfo(object)
         if (is.null(filePath)) filePath <- h5.meta$filename
+        if (is.null(filePath)) {
+            stop("No filename identified")
+        }
         if (!file.exists(filePath)) {
             stop("HDF5 file path does not exist:\n",
                  filePath)
         }
-        .log("Restoring ligerDataset from: ", filePath)
+        .log("filename identified: ", filePath)
         h5file <- hdf5r::H5File$new(filePath, mode = "r+")
         h5.meta$filename <- h5file$filename
         pathChecks <- unlist(lapply(h5.meta[4:10], function(x) {
@@ -228,16 +231,15 @@ restoreH5Liger <- function(object, filePath = NULL) {
         rawData(object, check = FALSE) <- h5file[[h5.meta$rawData]]
         if (!is.null(h5.meta$normData))
             normData(object, check = FALSE) <- h5file[[h5.meta$normData]]
-        if (!is.null(h5.meta$scaleData))
+        if (!is.null(h5.meta$scaleData)) {
             scaleData(object, check = FALSE) <- h5file[[h5.meta$scaleData]]
+        }
         methods::validObject(object)
     } else {
         # Working for liger object
         if (!is.null(filePath)) {
             if (!is.list(filePath) || is.null(names(filePath)))
                 stop("`filePath` has to be a named list for liger object.")
-            if (!any(names(filePath) %in% names(object)))
-                stop("names of `filePath` must be found in `names(object)`.")
         }
         for (d in names(object)) {
             if (isH5Liger(object, d)) {
@@ -249,7 +251,8 @@ restoreH5Liger <- function(object, filePath = NULL) {
                                 filePath[[d]])
                     else path <- filePath[[d]]
                 }
-                dataset(object, d, qc = FALSE) <-
+                .log("Restoring dataset \"", d, "\"")
+                datasets(object, check = FALSE)[[d]] <-
                     restoreH5Liger(dataset(object, d), filePath[[d]])
             }
         }
@@ -270,3 +273,39 @@ restoreH5Liger <- function(object, filePath = NULL) {
         list(folder = NULL, data = path)
     }
 }
+
+#' Close all links (to HDF5 files) of a liger object
+#' @description When need to interact with the data embedded in HDF5 files out
+#' of the currect R session, the HDF5 files has to be closed in order to be
+#' available to other processes.
+#' @param object liger object.
+#' @return \code{object} with links closed.
+#' @export
+closeAllH5 <- function(object) {
+    if (!isH5Liger(object)) return(object)
+    for (dn in names(object)) {
+        if (!isH5Liger(object, dn)) next
+        ld <- dataset(object, dn)
+        # if (!is.null(rawData(ld))) rawData(ld)$close()
+        # if (!is.null(normData(ld))) normData(ld)$close()
+        # if (!is.null(scaleData(ld))) scaleData(ld)$close()
+        # if (!is.null(scaleUnsharedData(ld))) scaleUnsharedData(ld)$close()
+        h5f <- getH5File(object, dn)
+        h5f$close_all()
+    }
+    return(object)
+}
+
+.H5GroupToH5SpMat <- function(obj, dims) {
+    groupPath <- obj$get_obj_name()
+    RcppPlanc::H5SpMat(filename = obj$get_filename(),
+                       valuePath = paste0(groupPath, "/data"),
+                       rowindPath = paste0(groupPath, "/indices"),
+                       colptrPath = paste0(groupPath, "/indptr"),
+                       nrow = dims[1], ncol = dims[2])
+}
+
+.H5DToH5Mat <- function(obj) {
+    RcppPlanc::H5Mat(filename = obj$get_filename(),
+                     dataPath = obj$get_obj_name())
+}

R/import.R

@@ -324,21 +324,24 @@ readLiger <- function(
         clusterName = "clusters",
         h5FilePath = NULL,
         update = TRUE) {
-    oldObj <- readRDS(filename)
-    if (!inherits(oldObj, "liger"))
+    obj <- readRDS(filename)
+    if (!inherits(obj, "liger"))
         stop("Object is not of class \"liger\".")
-    oldVer <- oldObj@version
-    if (oldVer >= package_version("1.99.0")) return(oldObj)
-    .log("Older version (", oldVer, ") of liger object detected.")
+    ver <- obj@version
+    if (ver >= package_version("1.99.0")) {
+        if (isH5Liger(obj)) obj <- restoreH5Liger(obj)
+        return(obj)
+    }
+    .log("Older version (", ver, ") of liger object detected.")
     if (isTRUE(update)) {
         .log("Updating the object structure to make it compatible ",
              "with current version (", utils::packageVersion("rliger2"), ")")
-        return(convertOldLiger(oldObj, dimredName = dimredName,
+        return(convertOldLiger(obj, dimredName = dimredName,
                                clusterName = clusterName,
                                h5FilePath = h5FilePath))
     } else {
         .log("`update = FALSE` specified. Returning the original object.")
-        return(oldObj)
+        return(obj)
     }
 }
 

R/liger-class.R

@@ -93,7 +93,6 @@ liger <- setClass(
         for (d in names(x)) {
             ld <- dataset(x, d)
             if (!is.null(ld@V)) {
-                print(all.equal(rownames(ld@V), varFeatures(x)))
                 if (!identical(rownames(ld@V), varFeatures(x)))
                     return(paste("Variable features do not match dimension",
                                  "of V matrix in dataset", d))
@@ -105,9 +104,18 @@ liger <- setClass(
                         return(paste("Variable features do not match dimension",
                                      "of scaleData in dataset", d))
                 } else {
-                    if (scaleData(ld)$dims[1] != length(varFeatures(x)))
-                        return(paste("Variable features do not match dimension",
-                                     "of scaleData in dataset (H5)", d))
+                    if (inherits(scaleData(ld), "H5D")) {
+                        if (scaleData(ld)$dims[1] != length(varFeatures(x)))
+                            return(paste("Variable features do not match ",
+                                         "dimension of scaleData in dataset ",
+                                         "(H5)", d))
+                    } else if (inherits(scaleData(ld), "H5Group")) {
+                        if (scaleData(ld)[["featureIdx"]]$dims != length(varFeatures(x))) {
+                            return(paste("Variable features do not match ",
+                                         "dimension of scaleData in dataset ",
+                                         "(H5)", d))
+                        }
+                    }
                 }
             }
         }
@@ -116,7 +124,7 @@ liger <- setClass(
 }
 
 .valid.liger <- function(object) {
-    message("Checking liger object validity")
+    # message("Checking liger object validity")
     res <- .checkLigerBarcodes(object)
     if (!is.null(res)) return(res)
     res <- .checkLigerVarFeature(object)
@@ -372,7 +380,7 @@ setGeneric("dataset", function(x, dataset = NULL) standardGeneric("dataset"))
 
 #' @export
 #' @rdname liger-class
-setGeneric("dataset<-", function(x, dataset, type = NULL, qc = TRUE, value) {
+setGeneric("dataset<-", function(x, dataset, type = NULL, qc = FALSE, value) {
     standardGeneric("dataset<-")
 })
 
@@ -404,6 +412,56 @@ setMethod("dataset", signature(x = "liger", dataset = "numeric"),
               datasets(x)[[dataset]]
           })
 
+.mergeCellMeta <- function(cm1, cm2) {
+    newDF <- S4Vectors::DataFrame(row.names = c(rownames(cm1), rownames(cm2)))
+    for (var in names(cm1)) {
+        value <- cm1[[var]]
+        if (var %in% names(cm2)) {
+            # TODO: check type
+            tryCatch(
+                expr = {
+                    if (is.null(dim(value))) {
+                        value <- c(value, cm2[[var]])
+                    } else {
+                        value <- rbind(value, cm2[[var]])
+                    }
+                },
+                finally = {
+                    cm2Idx <- seq(nrow(cm1) + 1, nrow(cm1) + nrow(cm2))
+                    if (is.null(dim(value))) value[cm2Idx] <- NA
+                    else {
+                        empty <- matrix(NA, nrow = nrow(cm2), ncol = ncol(value))
+                        value <- rbind(value, empty)
+                    }
+                }
+            )
+        } else {
+            cm2Idx <- seq(nrow(cm1) + 1, nrow(cm1) + nrow(cm2))
+            if (is.null(dim(value))) value[cm2Idx] <- NA
+            else {
+                empty <- matrix(NA, nrow = nrow(cm2), ncol = ncol(value))
+                value <- rbind(value, empty)
+            }
+        }
+        newDF[[var]] <- value
+    }
+    for (var in names(cm2)[!names(cm2) %in% names(cm1)]) {
+        value <- cm2[[var]]
+        if (is.null(dim(value))) {
+            if (is.factor(value)) {
+                value <- factor(c(rep(NA, nrow(cm1)), value),
+                                levels = levels(value))
+            } else {
+                value <- c(rep(NA, nrow(cm1)), value)
+            }
+        } else {
+            empty <- matrix(NA, nrow = nrow(cm1), ncol = ncol(value))
+            value <- rbind(empty, value)
+        }
+        newDF[[var]] <- value
+    }
+    return(newDF)
+}
 .expandDataFrame <- function(df, idx) {
     dfList <- as.list(df)
     dfList <- lapply(dfList, function(x, idx) {
@@ -428,19 +486,20 @@ setMethod("dataset", signature(x = "liger", dataset = "numeric"),
 setReplaceMethod("dataset", signature(x = "liger", dataset = "character",
                                       type = "missing", qc = "ANY",
                                       value = "ligerDataset"),
-                 function(x, dataset, type = NULL, qc = TRUE, value) {
+                 function(x, dataset, type = NULL, qc = FALSE, value) {
                      if (dataset %in% names(x)) {
                          dataset(x, dataset) <- NULL
                      }
                      new.idx <- seq(ncol(x) + 1, ncol(x) + ncol(value))
                      x@datasets[[dataset]] <- value
-                     # TODO also add rows to cellMeta and H.norm
                      cm <- x@cellMeta
                      remainingRowname <- rownames(cm)
                      cm <- .expandDataFrame(cm, new.idx)
                      rownames(cm) <- c(remainingRowname, colnames(value))
+
                      levels(cm$dataset) <- c(levels(cm$dataset), dataset)
                      cm$dataset[new.idx] <- dataset
+                     cm$barcode[new.idx] <- colnames(value)
                      x@cellMeta <- cm
                      # x@W is genes x k, no need to worry
                      if (!is.null(x@H.norm)) {
@@ -462,11 +521,12 @@ setReplaceMethod("dataset", signature(x = "liger", dataset = "character",
                                       value = "matrixLike"),
                  function(x, dataset,
                           type = c("rawData", "normData", "scaleData"),
-                          qc = TRUE,
+                          qc = FALSE,
                           value) {
                      type <- match.arg(type)
                      if (type == "rawData") {
                          ld <- createLigerDataset(rawData = value)
+                         colnames(ld) <- paste0(dataset, "_", colnames(ld))
                      } else if (type == "normData") {
                          ld <- createLigerDataset(normData = value)
                      } else if (type == "scaleData") {