MWE to demonstrate a bug in the calculation of non unique matches in synapter 1.5.3

griddetails_bug.R

library("synapter")
library("synapterdata")

###########################################################################
## some boring stuff to prepare the Synapter object
###########################################################################

fasfile <- synapterdata::getFasta()
mstrfile <- synapterdata::getMaster()
pepfile <- synapterdata::getMSeFinalPeptide()[2]
pep3dfile <- synapterdata::getMSePep3D()[2]

inlist <- list(quantpeptide = pepfile,
               quantpep3d = pep3dfile,
               identpeptide = mstrfile,
               fasta = fasfile)    

s <- Synapter(inlist, master=TRUE)

## apply default filtering
filterUniqueDbPeptides(s, verbose=TRUE) 

fdr <- fpr <- 0.05
ppm <- 20

## we could ignore the Quant/MSE part for now
filterIdentPepScore(s, method="BH", fdr=fdr)
filterQuantPepScore(s, method="BH", fdr=fdr)
filterIdentProtFpr(s, fpr=fpr)
filterQuantProtFpr(s, fpr=fpr)
filterIdentPpmError(s, ppm=ppm)
filterQuantPpmError(s, ppm=ppm)

mergePeptides(s)

modelRt(s, span=0.05)

###########################################################################
## END OF preprocessing
###########################################################################

## original getAssignmentDetails (some debugging output added)
getAssignmentDetails <- function(matchedEmrts) {
  x <- matchedEmrts[!is.na(matchedEmrts$precursor.leID.quant),
                    c("precursor.leID.quant",
                      "spectrumID",
                      "matched.quant.spectrumIDs")]
  table(apply(x, 1,
              function(x) {
                k <- unlist(strsplit(x["matched.quant.spectrumIDs"], ","))
                lk <- length(k)
                if (lk == 0) {
                  return(0)
                } else if (lk == 1) {
                  ifelse(x["spectrumID"] == x["precursor.leID.quant"], 1, -1)
                } else {
                  .x <- as.numeric(unlist(strsplit(x["matched.quant.spectrumIDs"], ",")))
                  res <- ifelse(x["precursor.leID.quant"] %in% .x, 2, -2)

                  ## add this lines for DEBUGGING purpose
                  id <- x["precursor.leID.quant"]
                  nid <- as.numeric(id)
                  nres <- ifelse(nid %in% .x, 2, -2)

                  if (nres != res) {
                    cat("precursor.leID.quant (", class(id), "):\n", sep="")
                    print(id)
                    cat("matched.quant.spectrumIDs:\n")
                    print(.x)
                    cat("reported result:", res, "(should be", nres, ")\n\n")
                  }

                  res
                }
              }))
}

###########################################################################

matchedEMRTs <- synapter:::findMSeEMRTs(s$IdentPeptideData,
                                        s$QuantPep3DData,
                                        s$MergedFeatures,
                                        nsd=3, ppm=9,
                                        s$RtModel, "transfer")

getAssignmentDetails(matchedEMRTs)
# precursor.leID.quant (character):
# precursor.leID.quant 
#              " 5622" 
# matched.quant.spectrumIDs:
# [1] 5469 5622
# reported result: -2 (should be 2 )
# 
# precursor.leID.quant (character):
# precursor.leID.quant 
#              "  252" 
# matched.quant.spectrumIDs:
# [1] 225 252
# reported result: -2 (should be 2 )
# 
# precursor.leID.quant (character):
# precursor.leID.quant 
#              " 9466" 
# matched.quant.spectrumIDs:
# [1] 9466 9650
# reported result: -2 (should be 2 )
# 
# 
#   -2   -1    0    1    2 
#    4  107   92 3470    9

## use the new implementation
source("utils.R")
gridSearch2(s$RtModel, s$IdentPeptideData, s$QuantPep3DData, 
            s$MergedFeatures[TRUE, ], ppms=9, nsds=3, verbose=FALSE)$details
# $`3.9`
#   -2   -1    0    1    2 
#    1  107   92 3470   12