inst/extcode/PeptideEnrich.R

In WubingZhang/PhageR: Functions for phage display data analysis

#' Select enriched peptides
#'
#' @docType methods
#' @name PeptideEnrich
#' @rdname PeptideEnrich
#'
#' @param rawcount Peptide count matrix.
#' @param ctl Index of controls.
#' @param trt Index of treatments.
#' @param pvalue The pvalue cutoff.
#' @param lfc The LogFC cutoff.
#'
#' @return Enriched peptide list.
#'
#' @author Wubing Zhang
#' @import ggView
#' @export

PeptideEnrich <- function(rawcount, ctl = 1, trt = 2, lfc = 3){

  require(ggView)
  if(length(c(ctl, trt))>3){
    SampleAnn = data.frame(Condition = rep(c("ctl", "trt"), c(length(ctl), length(trt))),
                           Sibs = rep(paste0("r", 1:length(ctl)), 2),
                           row.names = colnames(rawcount)[c(ctl,trt)], stringsAsFactors = FALSE)
    EDS = DEAnalyze(rawcount, SampleAnn, "RNASeq", "DESeq2", paired = TRUE)
    kmerNormCount = EDS@normlized
    kmerNormCount[kmerNormCount<0] = 0
    sigres = EDS@DEGRes
    # saveRDS(EDS, paste0(outdir, "/", proj, "_", k, "mer_summary.rds"))
    # write.table(sigres, paste0(outdir, "/", proj, "_", k, "mer_summary.txt"),
    #             row.names = TRUE, sep = "\t", quote = FALSE)
    sigres = sigres[order(-sigres$stat), ]
    enriched_peptides = rownames(sigres)[sigres$log2FC>lfc]
  }else{
    kmerNormCount = log2(1e6 * t(t(rawcount) / colSums(rawcount)) + 1)
    kmerNormCount = kmerNormCount[order(kmerNormCount[,ctl]-kmerNormCount[,trt]), ]
    idx = kmerNormCount[, trt]-kmerNormCount[, ctl] > lfc
    enriched_peptides = rownames(kmerNormCount)[idx]
  }
  return(list(normcount = kmerNormCount, enriched_peptides = enriched_peptides))
}