In fgcz/rawR: Direct Access to Orbitrap Data and Beyond

Configuiration parameters

| | | |-------------------------------|---------------------------------------| | file report release date: | r date() | | autoQC01 lm cache file: | r params[['input']] | | nodename: | r Sys.info()['nodename'] | | R version: | r R.version.string | | rawrr version: | r packageVersion('rawrr') | | rmarkdown version: | r packageVersion('rmarkdown') | | r.squared.cutoff.red: | r params[['r.squared.cutoff.red']] |

stopifnot(
require(knitr),
require(kableExtra),
require(lattice),
require(protViz),
require(colorspace)
)

.getIRTs <- function(){
  iRTpeptide <- c("LGGNEQVTR", "YILAGVENSK", "GTFIIDPGGVIR", "GTFIIDPAAVIR",
    "GAGSSEPVTGLDAK", "TPVISGGPYEYR", "VEATFGVDESNAK", "TPVITGAPYEYR",
    "DGLDAASYYAPVR", "ADVTPADFSEWSK", "LFLQFGAQGSPFLK")

  df <- protViz::iRTpeptides[protViz::iRTpeptides$peptide %in% iRTpeptide,]
  df$ssrc <- protViz::ssrc(as.character(df$peptide))
  df$mass <- protViz::parentIonMass(as.character(df$peptide))
  df$mass2Hplus <- (df$mass + 1.008) / 2
  df
}

.flm_autoQC01 <-
  function(filename, peptides = .getIRTs(), tol=10, r.squared.cutoff=0.98,
           fileprefix='/srv/www/htdocs'){

    rawfile <- file.path(fileprefix, filename)
    if(!file.exists(rawfile)){
      warning(paste0("file ", rawfile, "does not exist."))
      return()
    }
    start  <- as.numeric(Sys.time()) * 1000

    result <- tryCatch({
      message(paste0("fetching XICs for rawfile", rawfile, " ..."))
      XIC <- rawrr::readChromatogram(rawfile, peptides$mass2Hplus, tol = tol)
      t <- sapply(XIC, function(x){if(length(x$times) > 0){
        x$times[x$intensities == max(x$intensities)][1]}else{NA}})
      n <- length(t)
      intensity.max <- max(sapply(XIC, function(x){max(x$intensities)}))

      xx <- data.frame(rt = t, irtscore = peptides$rt)
      fm <- lm(rt ~ irtscore, data = xx, na.action=na.exclude)
      list(fm = fm,
           xx = xx,
           intensity.max = intensity.max,
           n=n,
           filename = filename,
           runtime = (as.numeric(Sys.time()) * 1000 - start) / 1000)
    },
    error = function(err) { NULL })
  } 

rv <- lapply("32b749642472_4590.raw", .flm_autoQC01,
    fileprefix = "/Users/cp/Library/Caches/org.R-project.R/R/ExperimentHub/") |>
  lapply(FUN=function(x){
    data.frame(r.squared = summary(x$fm)$r.squared,
               slope = x$fm$coefficients[2],
               intercept = x$fm$coefficients[1],
               n = x$n,
               intensity.max = x$intensity.max,
               filename = as.character(x$filename),
               runtime = x$runtime,
               row.names = NULL)
  }) |>
  Reduce(f=rbind)

.assignInstrument <- function(x){
  stopifnot(is.data.frame(x))

  x$instrument <- NA
  for (p in c("FUSION_1", "FUSION_2", "G2HD_1", "LC1100",
              "LTQ_1", "LTQFT_1", "ORBI_1", "ORBI_2", "PROTEONXPR36",
              "QEXACTIVE_1", "QEXACTIVE_2", "QEXACTIVE_3", "QEXACTIVEHF_1",
              "QEXACTIVEHF_2", "QEXACTIVEHF_4", "QEXACTIVEHFX_1", "QTRAP_1", "T100_1",
              "TOFTOF_2", "TRIPLETOF_1", "TSQ_1", "TSQ_2", "VELOS_1",
              "VELOS_2", "LUMOS_1", "LUMOS_2", "EXPLORIS_1")){
    x$instrument[grep(p, x$filename)] <- p
  }
  x
}

cv <- 1-2:7/10
t <- trellis.par.get("strip.background")
t$col <- (rgb(cv,cv,cv))
trellis.par.set("strip.background",t)

tp <- trellis.par.get("par.sub.text")
tp$cex <- 0.5
trellis.par.set("par.sub.text", tp)

if(file.exists(params[['input']])){
  autoQC01 <- read.table(params[['input']], header=TRUE, sep=';')
  autoQC01 <- assignInstrument(autoQC01)
  autoQC01$POSIXct <- as.POSIXct(autoQC01$time, origin="1970-01-01")
}

if(nrow(autoQC01) > 1){
  qc01.mean <- aggregate(slope ~ instrument, data=autoQC01, FUN=mean, subset=autoQC01$r.squared > params$r.squared.cutoff.red)
  qc01.sd <- aggregate(slope ~ instrument, data=autoQC01, FUN=sd, subset=autoQC01$r.squared > params$r.squared.cutoff.red)
  names(qc01.mean) <- c('instrument','instrument.slope.mean')
  names(qc01.sd) <- c('instrument','instrument.slope.sd')
  qc01.mean.sd <- merge( qc01.mean, qc01.sd)
  autoQC01 <- merge(autoQC01, qc01.mean.sd, by='instrument')
  autoQC01$slope.zscore <- (autoQC01$slope - autoQC01$instrument.slope.mean) / autoQC01$instrument.slope.sd
}

if(nrow(autoQC01) > 1){
  qc01.mean <- aggregate(intercept ~ instrument, data=autoQC01, FUN=mean, subset=autoQC01$r.squared > 0.9)
  qc01.sd <- aggregate(intercept ~ instrument, data=autoQC01, FUN=sd, subset=autoQC01$r.squared > 0.9)
  names(qc01.mean) <- c('instrument','instrument.intercept.mean')
  names(qc01.sd) <- c('instrument','instrument.intercept.sd')
  qc01.mean.sd <- merge( qc01.mean, qc01.sd)

  autoQC01 <- merge(autoQC01, qc01.mean.sd, by='instrument')
  autoQC01$intercept.zscore <- (autoQC01$intercept - autoQC01$instrument.intercept.mean) / autoQC01$instrument.intercept.sd
}

Last 24 hours

autoQC01.24h <- autoQC01[rev(which(Sys.time() - autoQC01$time < 3600 * 24 * 2)), ]

idx <- rev(order(autoQC01.24h$time))

autoQC01.24h <- autoQC01.24h[idx, c('POSIXct', 'filename', 'slope',
    'slope.zscore', 'intercept', 'intercept.zscore', 'r.squared')] 

QCred <- which(autoQC01.24h$r.squared < params$r.squared.cutoff.red)
QCyellow <- which(autoQC01.24h$r.squared < params$r.squared.cutoff.yellow)

kable(autoQC01.24h, row.names=FALSE) |>
  kable_styling("striped", full_width = FALSE) |> 
  row_spec(QCyellow, bold = TRUE, color = "black", background = params$r.squared.cutoff.yellow.rgb) |>
  row_spec(QCred, bold = TRUE, color = "white", background = params$r.squared.cutoff.red.rgb)

Fitted models

rv <- parallel::mclapply(autoQC01.24h$filename, FUN=.flm_autoQC01, mc.cores=16) 

rv <- lapply(rv, function(x, r.squared.cutoff=0.98){
    message(paste("plotting", x$filename, "..."))

        if (summary(x$fm)$r.squared < params$r.squared.cutoff.red){
             op <- par(mfrow = c(1,5), mar=c(5,5,5,1),bg=params$r.squared.cutoff.red.rgb, col='white')
        }else if (summary(x$fm)$r.squared < params$r.squared.cutoff.yellow){
             op <- par(mfrow = c(1,5), mar=c(5,5,5,1), bg=params$r.squared.cutoff.yellow.rgb, col='black')
        }else{
             op <- par(mfrow = c(1,5), mar=c(5,5,5,1))
        }

        plot(x$xx$rt ~ x$xx$irtscore, asp = 1, main = "rt ~ irtscore", type='n')
        legend("topleft", strsplit(as.character(x$filename), split = '/', perl = TRUE)[[1]][c(4,3,1)],box.col = 'transparent',cex = 0.75)
        legend.text<- paste(c('slope','intercept','r.squared'), as.numeric(round(c(coef(x$fm)["irtscore"], coef(x$fm)[1], summary(x$fm)$r.squared),3)),sep=": ")

        legend("bottomright", legend.text, cex=0.75, box.col = 'transparent')
        points(x$xx$irtscore, x$xx$rt, pch=16)
        abline(x$fm)

        plot(x$fm)
        par(op)
    })

---

Study slope and intercept

tp <- trellis.par.get("superpose.symbol")
tp$col <- colorspace::rainbow_hcl(8, alpha = 0.152)
trellis.par.set("superpose.symbol", tp)

xyplot(intercept ~ slope | instrument,
  group=instrument,
  data=autoQC01,
  subset=r.squared > params$r.squared.cutoff.yellow, 
  pch=16
  )

iRT peptides

df <- .getIRTs()
names(df) <- c("peptide", "iRT", "ssrc", "mass", "mZ")
df <- df[,c(1,2,5)]
kable(df) |>
  kable_styling("striped")

Session Info

sessionInfo()

fgcz/rawR documentation built on May 5, 2024, 3:46 p.m.