R/noiseq.R

In SoniaTC/NOISeq: Exploratory analysis and differential expression for RNA-seq data

noiseq <-
function (input, k = 0.5, norm = c("rpkm","uqua","tmm","n"),
          replicates = c("technical","biological","no"),
          factor=NULL, conditions = NULL, pnr = 0.2, nss = 5, v = 0.02, lc = 0)

# input: Object containing gene counts and as many columns as samples.

# k:      When counts = 0, 0 will be changed to k. By default, k = 0.5.
  
# norm:   Normalization method. It can be one of "rpkm" (default), "uqua"
#         (upper quartile), "tmm" (trimmed mean of M) or "n" (no normalization).

# factor: String with the factor to choose which conditions you are going
#         to compare.

# conditions: String with the conditions to compare in case one factor contains more
#         than 2 conditions.

# pnr:    Percentage of total reads (seq.depth) for each simulated sample.
#         Only needed when no replicates available. By default, pnr = 0.2.

# nss:    Number of simulated samples (>= 2). By default, nss = 5.
#         If nss = 0, real samples are used to compute noise.

# v:      Variability in sample total reads used to simulate samples.
#         By default, v = 0.02. Sample total reads is computed as a
#         random value from a uniform distribution in the interval
#         [(pnr-v)*sum(counts), (pnr+v)*sum(counts)]

# lc:     Length correction in done by dividing expression by length^lc.
#         By default, lc = 0. 

{
  
  if (inherits(input,"eSet") == FALSE)
    stop("Error. You must give an object generated by the readData function\n")

  if (is.null(factor))
    stop("Error. You must specify the factor to know which conditions you wish to compare.
         Please, give the argument 'factor'.\n")

  replicates <- match.arg(replicates)
  if (replicates == "biological") print("WARNING: Your experiment has biological replicates. You should consider using NOISeqBIO instead of NOISeq.")
  norm <- match.arg(norm)
  
  # (M,D) for signal and noise
  print("Computing (M,D) values...")
  miMD <- allMD(input, factor, k = k, replicates = replicates, 
                norm = norm, conditions = conditions,
                pnr = pnr, nss = nss, v = v, lc = lc)
 
    #------------------------------------------------------------------#
    
  
  
    ## Probability of differential expression

  print("Computing probability of differential expression...")
  prob.concounts <- probdeg(miMD$Ms, miMD$Ds, miMD$Mn, miMD$Dn)$prob

  if (!is.null(assayData(input)$exprs))
    todos <- rownames(as.matrix(assayData(input)$exprs))
  else
    todos <- rownames(as.matrix(assayData(input)$counts))
  
  prob <- prob.concounts[todos]
  names(prob) <- todos

    ## Results
  resultat <- data.frame("level_1" = miMD$Level1, "level_2" = miMD$Level2, 
                         "M" = miMD$Ms, "D" = miMD$Ds, "prob" = prob)
  rownames(resultat) <- todos

  # We change the name of the conditions to "name_mean"
  colnames(resultat)[1] <- paste(unlist(strsplit(miMD$comp," "))[1],"mean",sep="_")
  colnames(resultat)[2] <- paste(unlist(strsplit(miMD$comp," "))[3],"mean",sep="_")
  
  resultat <- data.frame(resultat, "ranking" = ranking(resultat)$statistic)
  if (!is.null(featureData(input)@data$Length))
    resultat <- data.frame(resultat, "Length" = as.numeric(as.character(featureData(input)@data[todos,"Length"])),
                           stringsAsFactors = FALSE)  
  if (!is.null(featureData(input)@data$GC))
    resultat <- data.frame(resultat, "GC" = as.numeric(as.character(featureData(input)@data[todos,"GC"])),
                           stringsAsFactors = FALSE)  
  if (!is.null(featureData(input)@data$Chromosome))
    resultat <- data.frame(resultat, "Chrom" = as.character(featureData(input)@data$Chromosome),
                           "GeneStart" = as.numeric(as.character(featureData(input)@data$GeneStart)),
                           "GeneEnd" = as.numeric(as.character(featureData(input)@data$GeneEnd)),
                           stringsAsFactors = FALSE)
  if (!is.null(featureData(input)@data$Biotype))
    resultat <- data.frame(resultat, "Biotype" = as.character(featureData(input)@data[todos,"Biotype"]),
                           stringsAsFactors = FALSE)
  #resultat[order(resultat[,5], decreasing = TRUE),]
  
  Output(data = list(resultat), method=norm,k=miMD$k,lc=lc,factor=factor,v=v,nss=nss,pnr=pnr,
         comparison=miMD$comp,replicates=replicates)
}