LFQ demo - Githubissues

code for tweet https://twitter.com/hb9feb/status/1014602529034915840
#R

library(protViz)
library(rawDiag)

f <- function(rawfile, pepSeq, dt = 0.1){
  mass2Hplus <- (parentIonMass(pepSeq) + 1.008) / 2
  X <- readXICs(rawfile = rawfile, masses = mass2Hplus)
  S <- read.raw(rawfile)

  idx <- lapply(mass2Hplus, function(m){
    which(abs(S$PrecursorMass - m) < 0.1)
  })

  scanNumbers <- lapply(idx, function(x){S$scanNumber[x]})

  bestMatchingMS2Scan <- sapply(1:length(pepSeq), function(i){
    peakList <- readScans(rawfile, scans = scanNumbers[[i]])

    peptideSpecMatch <- lapply(peakList,
                               function(x){
                                 psm(pepSeq[i], x, FUN = function (b, y){cbind(b, y)}, plot = FALSE)})
    score <- sapply(1:length(peptideSpecMatch), 
                    function(j){
                      sum(peakList[[j]]$intensity[abs(peptideSpecMatch[[j]]$mZ.Da.error) < 0.1])})
    bestFirstMatch <- which(max(score, na.rm = TRUE) == score)[1]
    scanNumbers[[i]][bestFirstMatch]
  })

  peakList <- readScans(rawfile, scans = bestMatchingMS2Scan)

  pp <- lapply(1:length(pepSeq), function(j){
    jpeg(filename = paste("~/Desktop/rawDiag_", pepSeq[j],".jpeg", sep=''), quality = 100, height = 640)
    op<-par(mfrow = c(2,1), mar = c(5,4,4,3))
    peakplot(pepSeq[j], peakList[[j]], FUN = function (b, y){cbind(b, y)})

    t <- S$StartTime[bestMatchingMS2Scan[j]];

    peak.idx  <- which((t - dt) < X[[j]]$times & X[[j]]$times < (t + dt))

    plot(X[[j]], xlim = c(t - 0.2, t + 0.2), main = paste("RT =", round(t * 60), 'seconds', "[m+2H]2+ =", mass2Hplus[j] ),
         xlab = 'RT [min]', ylab = 'intensity');
    abline(v = t, col = rgb(0.8, 0.1, 0.1, alpha = 0.5), lwd = 3)

    # peak fitting
    xx <- X[[j]]$times[peak.idx]
    yy <- X[[j]]$intensities[peak.idx]
    points(xx, yy, pch = 16, col = rgb(0.0, 0.1, 0.8, alpha = 0.5))
    # text(xx, yy, peak.idx, pos = 1)
    peak <- data.frame(logy = log(yy), x = xx)
    x.mean <- mean(peak$x)
    peak$xc <- peak$x - x.mean
    (fit <- lm(logy ~ xc + I(xc^2), data = peak))
    xx <- with(peak, seq(min(xc) - 0.2, max(xc) + 0.2, length = 100))
    lines(xx + x.mean, exp(predict(fit, data.frame(xc = xx))), col=rgb(0.25, 0.25, 0.25, alpha = 0.3), lwd = 5)
    dev.off()
  })

}
# https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3918884/

f(rawfile = "/Users/cp/Downloads/20180220_14_autoQC01.raw", 
  c('GAGSSEPVTGLDAK', 'VEATFGVDESNAK', 
    'TPVISGGPYEYR', 'TPVITGAPYEYR', 'DGLDAASYYAPVR',
    'ADVTPADFSEWSK', 'GTFIIDPGGVIR')
    )
fgcz / rawDiag

LFQ demo #25
