Read BioAnalyzer data into R for custom plotting

readBioanalyzer.R

library(drc) #used for fitting Weibull retardation model to get Time vs. Size

clean.up<-function(data){
  data<-data[-(nrow(data):(nrow(data)-1)),]
  data$Value<-as.numeric(as.vector(data$Value))
  data$Time<-as.numeric(as.vector(data$Time))
  return(data)
}

find_peaks <- function (data, nt = c(25, 200, 500, 1000, 2000, 4000, 6000),  m = 50, k = 2){
  x<-data$Value
  x[x<k]<-0
  shape <- diff(sign(diff(x, na.pad = FALSE)))
  pks <- sapply(which(shape < 0), FUN = function(i){
    z <- i - m + 1
    z <- ifelse(z > 0, z, 1)
    w <- i + m + 1
    w <- ifelse(w < length(x), w, length(x))
    if(all(x[c(z : i, (i + 2) : w)] <= x[i + 1])) return(i + 1) else return(numeric(0))
  })
  pks <- unlist(pks)
  pks
}

#' readBioanalyzer
#'
#' Reads in bioAnalyzer data for custom plotting
#'
#' @param biosample full file path to the Sample CSV file, character.
#' @param ladder full file path to the Ladder CSV file, character.
#' @param nt the nucleotide size of the ladder, integer vector.
#' @param m sliding window in finding peaks, integer, default is 50 time samples.
#' @param k noise floor of the bioAnalyzer (everything below is set to 0), integer, default is 2.
#' @param plot if diagnostic plot should be plotted, boolean, default is TRUE.
#' 
#' @return data.frame with nucleotides nt and fluorescence (numeric).
#'
#' @examples
#' data <- readBioanalyzer(biosample = 'SAMPLE_2020-02-06_15-28-26_Sample3.csv', ladder = 'SAMPLE_2020-02-06_15-28-26_Ladder.csv')
#'
#' @export
readBioanalyzer<-function(biosample, ladder, nt = c(25, 200, 500, 1000, 2000, 4000, 6000),  m = 50, k = 2, plot = TRUE){
  
  sampleName <- basename(biosample)
  
  ladder<-read.table(ladder, skip = 17, fill=TRUE, header = TRUE, sep=',')
  
  biosample<-read.table(biosample, skip = 17, fill=TRUE, header = TRUE, sep=',')
  

  
  ladder<-clean.up(ladder)
  biosample<-clean.up(biosample)
  

  
  index<-find_peaks(ladder, m = m, k = k, nt = nt)
  nt<-c(0, nt)
  Time<-c(0, ladder$Time[index])
  
  
  
  model <- drm(nt ~ Time, fct = W2.4())
  par(yaxs="r", xaxs="r", mfrow=c(2,2))
  plot(model, ylab='Nucleotides', xlab='Time', pch=21, bg='red', main='Electrophoresis retardation model fit', cex.main=0.9 )
  x<-as.numeric(model$curve[[1]](ladder$Time) )-min(ladder$Time)
  
  plot(sqrt(x), ladder$Value, type='l', xlab='Nucleotides', axes = FALSE, ylab='Fluorescence [F.U.]', ylim=c(-2, max(ladder$Value)*1.1), main= 'Ladder', cex.main=0.9)
  axis(1, at=c(sqrt(c(nt, 10000))), labels=c(nt, 10000) )
  axis(2, las=1)
  points(sqrt(x)[index], ladder$Value[index], pch=21, bg='red')
  box()
  
  plot(sqrt(x), biosample$Value[seq_along(x)], type='l', xlab='Nucleotides', axes = FALSE, ylab='Fluorescence [F.U.]', ylim=c(-0.2, max(biosample$Value)*1.1), main= sampleName, cex.main=0.6)
  axis(1, at=c(sqrt(c(nt, 10000))), labels=c(nt, 10000) )
  axis(2, las=1)
  
  
  return( data.frame(nt = x, fluorescence = biosample$Value[seq_along(x)]) )

}