R/base_mv_tolerance.R

In huangwb8/GSClassifier: Gene Signature Classifier

#' @title mv_tolerance
#' @description Missing value tolerance analysis for \code{GSClassifier} models
#' @param X a data frame with the number of colume ≥ 2
#' @param gene.loss a vector of the number of genes masked
#' @param model a character of GSClassifier model or a \code{luckyModel} model
#' @param seed random seed
#' @param verbose whether to report the process
#' @importFrom pROC multiclass.roc
#' @return A list containing results of subtype identification and multi-roc analysis
#' @seealso \code{\link[pROC]{multiclass.roc}}.
#' @author Weibin Huang<\email{hwb2012@@qq.com}>
#' @examples
#' Empty
#' @export
mv_tolerance <- function(X,
                         gene.loss = c(2, 4, 6, 8, 10, 12),
                         levels = c(1, 2, 3, 4),
                         model = 'PAD.train_20220916',
                         seed = 487,
                         verbose = T) {
  # Test
  if (F) {
    # Internal validation cohort
    testData <-
      readRDS(system.file("extdata", "testData.rds", package = "GSClassifier"))
    expr_pad <- testData$PanSTAD_expr_part
    modelInfo <- modelData(
      design = testData$PanSTAD_phenotype_part,
      id.col = "ID",
      variable = c("platform", "PAD_subtype"),
      Prop = 0.7,
      seed = 19871
    )
    validInform <- modelInfo$Data$Valid
    X <- expr_pad[, validInform$ID]

    # Other parameters
    gene.loss = c(2, 4, 6, 8, 10, 12)
    model = 'PAD.train_20220916'
    seed = 487
    verbose = T
    levels = c(1, 2, 3, 4)

  }

  # MVI with quantile algorithm
  X <- GSClassifier:::na_fill(X,
                              method = 'quantile',
                              seed = 411,
                              verbose = verbose)

  # Get model
  if (is.character(model)) {
    # Model in GSClassifier
    m <- readRDS(system.file("extdata",
                             paste0(model, '.rds', collapse = ''),
                             package = "GSClassifier"))
  } else if (is.list(model)) {
    # Model in luckyModel
    m <- model
  } else {
    # Error report
    stop('Please use a right model format.')
  }

  # Subtype identification for raw matrix
  res0 <- parCallEnsemble(
    X = X,
    ens = m$ens$Model,
    geneAnnotation = m$geneAnnotation,
    geneSet = m$geneSet,
    scaller = m$scaller$Model,
    geneid = "ensembl",
    subtype = NULL,
    verbose = verbose,
    numCores = 4
  )

  # multi-ROC: time-consuming if a large matrix used
  set.seed(seed)
  seeds <- sample(1:100000,
                  length(gene.loss),
                  replace = F)
  mAUC <- list()
  model_res <- list()
  for (i in 1:length(gene.loss)) {
    # i=1

    # Masked matrix with zero value
    X2 <- masked_matrix(X,
                        gene.loss = gene.loss[i],
                        seed = 487)

    # Subtype identification
    res <- parCallEnsemble(
      X = X2,
      ens = m$ens$Model,
      geneAnnotation = m$geneAnnotation,
      geneSet = m$geneSet,
      scaller = m$scaller$Model,
      geneid = "ensembl",
      subtype = NULL,
      verbose = verbose,
      numCores = 4
    )

    mAUC[[paste0('GeneLoss=', gene.loss[i])]] <- multiclass.roc(
      response = res0$BestCall,
      predictor = res$BestCall,
      levels = levels,
      quiet = T
    )
    model_res[[paste0('GeneLoss=', gene.loss[i])]] <- res
  }
  # mymusic()

  # Output
  l <- list(multiAUC = mAUC,
            modelRes = model_res)
  return(l)

  # saveRDS(l, './data/results-mv_tolerance-GC.rds')

}

#' @inheritParams mv_tolerance
masked_matrix <- function(X,
                          gene.loss = 2,
                          seed = 487) {
  # Zero Matrix
  # Xzero <- matrix(
  #   data = rep(0, nrow(X)*ncol(X)),
  #   nrow = nrow(X), byrow = F,
  #   dimnames = list(rownames(X),
  #                   colnames(X)))

  # Random seeds
  set.seed(seed)
  seeds <- sample(1:100000, ncol(X), replace = F)
  X2 <- as.matrix(X)
  for (i in 1:ncol(X)) {
    # i=1
    set.seed(seeds[i])
    X2[sample(1:nrow(X2), gene.loss, replace = F), i] <- 0
  }
  return(X2)
}