R/qPCR_rxns.R

In dimitriskokoretsis/qpcrR:

#' @import data.table
#'
#' @export

qPCR_rxns <- function(genes,
                      sample,biol.rep,tech.rep.unkn,
                      std.curve,tech.rep.std.dilution,std.curve.sample,
                      NTC) {

  if(typeof(genes)!="character") {
    genes <- as.character(genes)
    warning("Argument genes was coerced to 'character'.")
  }

  if(length(NTC)==1) {
    NTC <- rep(NTC,times=length(genes))
  }

  n.rxns <- sum(tech.rep.unkn) * sum(biol.rep) + # unknown reactions
    length(genes) * sum(tech.rep.std.dilution) + # std. curve reactions
    sum(NTC) # ntc reactions

  # First, make a data frame with one row for each reaction
  # Fields: sample, rnx.type(ntc,std,unkn) biol.rep OR dilution, tech.rep, plate, row, column

  rxns <- data.table(sample=character(n.rxns),
                     rxn.type=character(n.rxns),
                     target=character(n.rxns),
                     biol.repl=integer(n.rxns),
                     dilution=numeric(n.rxns),
                     tech.repl=integer(n.rxns))

  i <- 1

  for(gene in seq_len(length(genes))) {

    # Unknown reactions

    rxns[i:(i+tech.rep.unkn[gene]*sum(biol.rep)-1),
         c("rxn.type","target","sample"):=
           .("Unkn",
             rep(..genes[..gene],
                 times=..tech.rep.unkn[..gene]*sum(..biol.rep)),
             rep(rep(..sample,
                     times=..biol.rep),
                 times=..tech.rep.unkn[..gene]))]

    rxns[target==genes[gene] & rxn.type=="Unkn",
         c("biol.repl","dilution"):=
           .(rep(1:..biol.rep[..sample==sample],
                 times=..tech.rep.unkn[..gene]),
             NA),
         by=sample]

    rxns[target==genes[gene] & rxn.type=="Unkn",
         tech.repl:=
           1:..tech.rep.unkn[..gene],
         by=.(sample,biol.repl)]


    # Standard curve reactions

    rxns[(i+tech.rep.unkn[gene]*sum(biol.rep)):
           (i+tech.rep.unkn[gene]*sum(biol.rep) + sum(tech.rep.std.dilution)-1),
         c("rxn.type","target","sample"):=
           .("Std",
             rep(..genes[..gene],
                 times=sum(..tech.rep.std.dilution)),
             rep(..std.curve.sample,
                 times=sum(..tech.rep.std.dilution)))]

    rxns[target==genes[gene] & rxn.type=="Std",
         c("dilution","biol.repl"):=
           .(rep(..std.curve,
                 times=..tech.rep.std.dilution),
             NA)]

    for(dil in 1:length(std.curve)) {

      # Fill technical replicates for each dilution

      rxns[target==genes[gene] & rxn.type=="Std" & dilution==std.curve[dil],
           tech.repl:=1:..tech.rep.std.dilution[..dil]][]
    }


    # Non-template control reactions
    if(NTC[gene]==TRUE) {

      rxns[i+tech.rep.unkn[gene]*sum(biol.rep) + sum(tech.rep.std.dilution),
           c("rxn.type","target","sample","biol.repl","dilution","tech.repl"):=
             .("NTC",..genes[..gene],NA,NA,NA,NA)][]

    }

    i <- i +
      tech.rep.unkn[gene] * sum(biol.rep) +
      sum(tech.rep.std.dilution) +
      NTC[gene]

  }

  setorder(rxns,target,rxn.type,sample,biol.repl,-dilution,tech.repl)

  return(rxns)
}