R/querries.R

In eprifti/momr: Mining Metaomics Data (MetaOMineR)

#' \code{projectOntoMGS}
#' @title projectOntoMGS
#' @export
#' @description This function takes a list of genes and projects it to a mgs catalogue
#' @author Edi Prifti & Emmanuelle Le Chatelier
#' @param genebag : vector of gene_ids
#' @param list.mgs : this is the structured MGS information formatted as a list of geneids each corresponding to an MGS
#' @param res.filt.mode : the filtering method, either 'perc' or 'size', default will be size.
#'        perc : genes projected onto mgs are only kept if they exceed a threshold percentage of the mgs size
#'        size : genes projected onto mgs are only kept if they exceed a threshold number of genes
#' @param res.filt.threshold : the threshold used to retain an mgs according to the filtering mode, either 'perc' or 'size'
#'        if 'perc' : the minimal percentage of genes projected onto mgs (number of genes projected / mgs size) needed to select a mgs
#'        if 'size' : the minimal number of genes projected onto mgs
#' @param not_projected : whether a last genebag containing not projected genes is to be added to the list of selected mgs
#' @return a list of selected mgs with the geneids according to the projected genes
projectOntoMGS <- function (genebag,
                            list.mgs,
                            res.filt.mode = "size",
                            res.filt.threshold = 50,
                            not_projected = TRUE) {
  res <- list()
  if (res.filt.mode != "perc" & res.filt.mode != "size") {
    stop("Choose 'perc' or 'size' as a filtering mode of the result")
  }
  else {
    for (m in 1:length(list.mgs)) {
      tmp <- genebag[(genebag %in% list.mgs[[m]])]
      if (res.filt.mode == "size") {
        if (length(tmp) >= res.filt.threshold) {
          res[[names(list.mgs)[m]]] <- tmp
        }
      }
      else {
        if (length(tmp) >= (length(list.mgs[[m]]) * res.filt.threshold / 100)) {
          res[[names(list.mgs)[m]]] <- tmp
        }
      }
    }
    if (not_projected == TRUE) {
      if (length(unlist(res)) == 0) {
        res[["not_projected"]] <- genebag
      } else {
        res[["not_projected"]] <- genebag[-which(genebag %in% unlist(res))]
      }
    }
    return(res)
  }
}

#' \code{extractProfiles}
#' @title extractProfiles
#' @export
#' @description This function extracts the profiles from a gene profile matrix of a group of genes or a list
#' of groups of genes. It can also restrict the size of the result
#' @author Edi Prifti & Emmanuelle Le Chatelier
#' @param genebag : vector or list of gene_ids
#' @param data : raw count or frequency matrix with genes_ids as rawnames
#' @param size.max : default 15000, maximum number of profiles to be extracted
#' @param size.min : default 1, the minimal size threshold above which a group of genes is selected.
#' This is used to extract multiple profiles and filtering the list with a minimal number of genes
#' @param silent : default TRUE, detailling and following computation progress
#' @return a matrix or a list of profile matrixes
extractProfiles <- function (genebag,
                             data,
                             size.max = 15000,
                             size.min = 1,
                             silent = TRUE)
{
  if (!is.list(genebag))
  {
    print("Profile extraction from a single gene vector")
    if (length(genebag) < size.max)
    {
      ind <- match(genebag, rownames(data))
      if (sum(is.na(ind)) > 0)
      {
        warning(paste(
          "There ",
          sum(is.na(ind)),
          "are genes from the genebag not found in the data !"
        ))
        ind <- ind[!is.na(ind)]
      }
      res <- data[ind, ]
    } else
    {
      ind <- match(genebag[1:size.max], rownames(data))
      if (sum(is.na(ind)) > 0)
      {
        warning(paste(
          "There ",
          sum(is.na(ind)),
          "are genes from the genebag not found in the data !"
        ))
        ind <- ind[!is.na(ind)]
      }
      res <- data[ind, ]
    }
  }
  else
  {
    print("Multiple profile extraction")
    res <- list()
    for (i in 1:length(genebag))
    {
      if (i %% 10 == 0 & silent == FALSE)
        print(i)
      if (length(genebag[[i]]) < size.max)
      {
        ind <- match(genebag[[i]], rownames(data))
        if (sum(is.na(ind)) > 0)
        {
          warning(paste(
            "There ",
            sum(is.na(ind)),
            "are genes from the genebag not found in the data !"
          ))
          ind <- ind[!is.na(ind)]
        }
        res[[i]] <- data[ind, ]
      } else
      {
        ind <- match(genebag[[i]][1:size.max], rownames(data))
        if (sum(is.na(ind)) > 0)
        {
          warning(paste(
            "There ",
            sum(is.na(ind)),
            "are genes from the genebag not found in the data !"
          ))
          ind <- ind[!is.na(ind)]
        }
        res[[i]] <- data[ind, ]
      }
    }
    names(res) <- names(genebag)
    res <- res[as.numeric(summary(genebag)[, "Length"]) >= size.min]
  }
  return(res)
}


#' \code{aggregateProfiles}
#' @title aggregateProfiles
#' @export
#' @description This function takes a list of profile matrixes and returns an aggregated big matrix.
#' The individual matrixes can be filtered in size so that the first X rows are used for each of them.
#' This function is used to prepare the data and plot different MGS as barcodes.
#' @author Emmanuelle Le Chatelier
#' @param list.profiles : list of matrix profiles
#' @param max.size : the maximum number of rows to be selected in the final agregated matrix.
#' default   max.size=25.
#' @param min.size : this is the minimum number of rows rows to be selected in the final aggregated matrix.
#' If a group has less it will be discarded. By default min.size=max.size.
#' @return an aggregated profile matrix.
aggregateProfiles <- function(list.profiles,
                              max.size = 25,
                              min.size = max.size)
{
  if (!is.list(list.profiles))
  {
    stop("Error! The object is not a list")
  }
  # list of matrices to be returned
  res <- c()
  for (i in 1:length(list.profiles))
  {
    if (nrow(list.profiles[[i]]) >= min.size)
    {
      res <- rbind(res, list.profiles[[i]][1:min(max.size, nrow(list.profiles[[i]])), ])
    }
  }
  return(res)
}

#' \code{computeFilteredVectors}
#' @export
#' @title computeFilteredVectors
#' @description filters and computes verctors based on gene profiles from a single matrix or a
#' list of matrix profiles
#' @author Edi Prifti & Emmanuelle Le Chatelier
#' @param profile : list of (or unique) matrix profiles
#' @param type : vectorisation method, vectors are calculated from the mean, median or sum of a given list of genes;
#' default type="mean" otherwise it will be "median" or  "sum"
#' @param filt : filtering threshold in % , rate of positive value in an individual under which matrix is cleaned
#' and sparse values put to 0 default filt= 0, no filtering
#' @param debug : default is FALSE, when TRUE information on advancement is printed
#' @return a filtered vector or a matrix of filtered vectors
computeFilteredVectors <- function (profile,
                                    type = "mean",
                                    filt = 0,
                                    debug = FALSE) {
  if (is.list(profile)) {
    res <- matrix(
      data = NA,
      ncol = ncol(profile[[1]]),
      nrow = length(profile)
    )
    for (i in 1:length(profile)) {
      if (debug)
        if (i %% 100 == 0) {
          print(i)
        }
      if (type == "mean") {
        # mean
        res[i, ] <- apply(filterMat(as.matrix(profile[[i]]), filt = filt), 2, mean)
      } else {
        if (type == "median") {
          # median
          res[i, ] <- apply(filterMat(as.matrix(profile[[i]]), filt = filt), 2, median)
        } else {
          #sum
          res[i, ] <- apply(filterMat(as.matrix(profile[[i]]), filt = filt), 2, sum)
        }
      }
    }
    rownames(res) <- names(profile)
    colnames(res) <- colnames(profile[[1]])
  }
  else {
    if (type == "mean") {
      # mean
      res <- apply(filterMat(as.matrix(profile), filt = filt), 2, mean)
    } else {
      if (type == "median") {
        # median
        res <- apply(filterMat(as.matrix(profile), filt = filt), 2, median)
      } else {
        # sum
        res <- apply(filterMat(as.matrix(profile), filt = filt), 2, sum)
      }
    }
  }
  return(res)
}

#' \code{buildMgsFinal}
#' @title buildMgsFinal
#' @export
#' @description This function will take a vector of genes (to be transformed into a list of genebags)
#'      or a list of genebags and will extract the profiles. Next genes well be ordered by connectivity
#'      which is to be computed for each group and the 50 most connected are selected to consitute the
#'      marker genes. These will be then used to compute the mean vectors. A final object containing all
#'      this information along with taxonomical annotation will be returned
#' @author Edi Prifti
#' @param genebag : a vector of genes to be projected onto mgs or a list of genebags, default = NULL.
#' @param mgs.cat : MGS catalogue to be used
#' @param mgs.taxo : taxonomy table for the MGS catalogue
#' @param profiles :  the data profile matrix to extract the profiles
#' @param conn : if TRUE the connectivity of a group is to be computed and ordered, default = TRUE.
#' @param silent : print detailled information on progress, default = FALSE.
#' @param filt : filtering based on percentage of prevalence to avoid noise for no signal samples by computeFilteredVectors.
#' @param save : Save files after each step
#' @return a list containing the final elements such as the 50 most connected genes, the mean vectors etc
buildMgsFinal <- function (genebag = NULL,
                           mgs.cat,
                           mgs.taxo,
                           profiles,
                           conn = TRUE,
                           silent = TRUE,
                           filt = 10,
                           save = FALSE)
{
  if (!is.list(genebag))
  {
    if (!silent)
    {
      print("           Genebag is not a list, projecting onto the MGS catalogue ...")
    }
    genebag.list <- projectOntoMGS(
      genebag = genebag,
      list.mgs = mgs.cat,
      res.filt.mode = "size",
      res.filt.threshold = 50
    )
  } else {
    genebag.list <- genebag
  }
  
  genebag.list.dat <- extractProfiles(genebag.list, profiles, silent = silent)
  
  if (!silent)
  {
    print("           Profiles have been extracted for each MGS ...")
  }
  
  if (save)
  {
    saveRDS(genebag.list.dat, file = "genebag.list.dat.rda", compress = FALSE)
  }
  
  if (!silent & save) {
    print("           Profiles file has been saved in the working directory ...")
  }
  
  genebag.list.ordered <- list()
  genebag.list.50 <- list()
  
  if (!silent & conn)
  {
    print("           Profiles will be reordered for each MGS using connectivity ...")
  }
  
  for (i in 1:length(genebag.list))
  {
    if (i %% 100 == 0 & !silent)
      print(i)
    
    dat <- genebag.list.dat[[i]]
    
    if (is.null(dim(dat))) {
      nr <- 1
    }
    else {
      nr <- nrow(dat)
    }
    
    if (nr == 0)
    {
      next()
    }
    
    if (conn == TRUE) {
      dat <- filterMat(as.matrix(dat), filt = filt)
      con <- connectivity(prof = dat,
                          soft = TRUE,
                          method = "spearman")
      dat <- dat[order(con, decreasing = T), ]
    }
    
    genebag.list.ordered[[i]] <- dat
    
    size.max <- min(50, nr)
    if (nr == 1)
    {
      genebag.list.50[[i]] <- genebag.list.ordered[[i]]
    } else
    {
      genebag.list.50[[i]] <- genebag.list.ordered[[i]][1:size.max, ]
    }
  }
  
  
  if (!silent & conn)
  {
    print("           Profiles re-ordering, filtering has been finished ...")
  }
  if (!silent & !conn)
  {
    print("           Profiles filtering has been finished ...")
  }
  names(genebag.list.ordered) <- names(genebag.list.dat)
  names(genebag.list.50) <- names(genebag.list.dat)
  if (save & conn)
  {
    saveRDS(genebag.list.ordered,
            file = "genebag.list.ordered.rda",
            compress = FALSE)
  }
  if (save)
  {
    saveRDS(genebag.list.50, file = "genebag.list.50.rda", compress = FALSE)
  }
  
  without.data <- unlist(lapply(lapply(genebag.list.50, nrow), is.null))
  
  res <- list()
  res$genebags <- genebag.list[!without.data]
  res$genebags.all <- mgs.cat[names(genebag.list)[!without.data]]
  res$mgs.profiles <- genebag.list.ordered[!without.data]
  res$mgs.50 <- genebag.list.50[!without.data]
  res$mean_vectors <- computeFilteredVectors(profile = res$mgs.50,
                                             type = "mean",
                                             filt = filt)
  if (!silent)
  {
    print("           The filtered tracer vectors have been created successfully ...")
  }
  res$taxonomy <- as.data.frame(mgs.taxo[names(res$genebags), ])
  name <- as.character(res$taxonomy$species)
  name <- paste(names(res$genebags), name)
  name <- gsub(" NA", "", name)
  res$taxonomy$name <- name
  if (!silent)
  {
    print("           All operations are finished ...")
  }
  return(res)
}


#' \code{connectivity}
#' @title connectivity
#' @description This function computes the intra-row correlation and applied a threshold to compute the connections of each row.
#' A connectivity vector is returned.
#' @author Emmanuelle le Chatelier & Edi Prifti
#' @param prof : a profile matrix. This data is used to compute correlations.
#' @param method : the correlation method, default = pearson.
#' @param th : default 0, if >0 than this threshold will be applied to compute a hard threholded connectivity.
#' @param soft : default = FALSE, if TRUE and when th > 0, the connectivity is computed as the soft threholded, sum of correlations
#' above the threshold
#' @return a vector of connectivity
#' @note this connectivity does not use a hard thresholding but is based on a total correlation score
connectivity <- function (prof,
                          method = "pearson",
                          th = 0,
                          soft = FALSE)
{
  if (th < 0 | th >= 1) {
    stop("The threhold should be positive and smaller than 1.")
  }
  else {
    data.cor <- Hmisc::rcorr(t(prof), type = method)$r
    if (soft) {
      data.cor[data.cor <= th] <- 0                 #
      con <- colSums(data.cor, na.rm = TRUE)        #
    }
    else {
      con <- colSums(data.cor > th, na.rm = TRUE)   #
    }
  }
  return(con)
}


#' \code{selectListSize}
#' @title selectListSize
#' @export
#' @description This function will sextract a part of a list based on the length of its components. Typically
#' a genebag list can be used. This function will work for uniclass lists of vectors and data frames.
#' @author Edi Prifti
#' @param l : a list of vectors or matrixes.
#' @param min.size : default = 0 and this has no action. This is the lower threshold for the element's size.
#' @param max.size : default = 0 and this has no action. This is the upper threshold for the element's size.
#' @param names : default = FALSE the function will return a subset of the list, if TRUE the function will
#' return only the names of the elements of the list as identifiers.
#' @return the trimmed subset of the list.
selectListSize <- function(l,
                           min.size = 0,
                           max.size = 0,
                           names = FALSE) {
  if (min.size == 0 & max.size == 0) {
    stop("nothing to be done")
  }
  if (min.size > max.size &
      max.size != 0) {
    stop("minimum value should be lower than maximum value")
  }
  if (min.size != 0 | max.size != 0) {
    type <- table(unlist(lapply(l, class)))
    if (length(type) == 1) {
      if (names(type) == "matrix" | names(type) == "data.frame") {
        l.size <- as.numeric(lapply(l, nrow))
      } else {
        l.size <- as.numeric(lapply(l, length))
      }
    } else {
      stop("multi class list case. not handled")
    }
    
    if (min.size != 0 & max.size != 0) {
      l.ind <- which(l.size >= min.size & l.size <= max.size)
    } else if (max.size != 0) {
      l.ind <- which(l.size <= max.size)
    } else if (min.size != 0) {
      l.ind <- which(l.size >= min.size)
    }
    if (names) {
      return(names(l[l.ind]))
    } else{
      return(l[l.ind])
    }
  }
}


#' \code{profiles2Genebags}
#' @title profiles2Genebags
#' @description This function will extract from a list of group profiles a list of identifiers of the matrix.
#' @author Edi Prifti
#' @param profiles : a list of dataframes
#' @return a list of genebags
profiles2Genebags <- function(profiles) {
  if (!is.list(profiles)) {
    stop("Error ! A list of profiles is needed.")
  } else{
    if (dim(profiles[[1]]) < 2) {
      stop("Error ! the elements of the list should be data frames or matrixes.")
    } else{
      return(lapply(profiles, rownames))
    }
  }
}

#' End of section and file