functions.r

do_msgf_runs = function(parameters,fastapath, mgfpath = fastapath,out = './msgfoutput',modfile = './modification.txt',msgfpluspath = './MSGFPlus/MSGFPlus.jar'){
  ## takes all files with .fasta extension as fastafiles
  ## takes all files with .mgf extension as mgf files
  dbsnames = list.files(fastapath, pattern = 'fasta')
  dbspaths = list.files(fastapath,full.names = TRUE, pattern = 'fasta')
  mgfpaths = list.files(mgfpath,full.names = TRUE,pattern = '.mgf')

  dir.create(out)

  a = Sys.time()
  f = lapply(1:length(dbsnames),
             function(db){
               output = paste0(out,'/',dbsnames[db],'.mzid')

               systemCall = paste("java -Xmx27000M -jar",msgfpluspath ,"-s", mgfpaths,
                                  '-o', output, '-d', dbspaths[db], '-mod', modfile, parameters)
               system(systemCall)
      file.remove(list.files(fastapath,
                                 full.names = TRUE,pattern = 'canno|cnlcp|cseq|csarr|rev'))
             })
  b = Sys.time()
  cat("start:\n"); print(a)
  cat("end:\n"); print(b)
  print(b-a)
}

library(tidyverse)
library(mzR)
library(stringr)
parse_msgf_mzid = function(mzid_path){
  d = openIDfile(mzid_path)
  if (!any(grepl('MS-GF+',mzidInfo(d)$software)))
    stop('Only mzID files generated by MS-GF+ can be parsed')
  t =  cbind(psms(d),score(d)[,-1]) %>% as_data_frame %>%
    filter(rank == 1) %>%
    transmute(spec_id = acquisitionNum,
              sequence = as.character(sequence),
              protein_id = as.character(DatabaseAccess),
              protein_id = str_replace(protein_id, '(XXX_)',''), ##remove XXX_ in decoy protein
              score= MS.GF.SpecEValue,
              database = gsub(".*/+([^/]+).fasta","\\1",levels(database(d)[['location']])),
              decoy = isDecoy,
              database_size = database(d)$numDatabaseSequences) %>%
    distinct
}

library(stringr)
is_subset = function(protein_id,fastapath){
  fas = read_lines(fastapath)
  headers = fas[grepl('>',fas,fixed = TRUE)]
  ## reducing search time by only searching the unique
  protein_id_unique = unique(protein_id)
  is_subset = sapply(protein_id_unique,function(p){any(str_detect(headers,fixed(p)))})
  protein_id %in% protein_id_unique[is_subset]
}

preprocess = function(dat,remove_target_decoy_PSM = TRUE, remove_multiple_proteins_PSM = FALSE){
  if (remove_target_decoy_PSM == TRUE){
    dat = group_by(dat,spec_id) %>% filter(!(any(decoy) & any(!decoy)))
  }
  if (remove_multiple_proteins_PSM == TRUE){
    dat = group_by(dat,spec_id) %>% filter(n() > 1)
  }
  dat = group_by(dat,spec_id,sequence,decoy,database,database_size,score) %>%
    summarise(subset = any(subset),non_subset = any(!subset),
              protein_id = paste0(protein_id,collapse = ';')) %>%
    ungroup
  dat
}

makeplots = function(df, title = FALSE){
  if (title){
  par(mfrow = c(2,2),oma = c(0, 0, 3, 0),mar = c(4,5,3.5,1))
    mtxt = "Histograms of MS-GF+ target scores"
    mtext(paste0(" \n\n complete library: ",subset_complete_n," proteins",
                 "\nsubset library: ",subset_n,
                 ' proteins related to ',subset_name), outer = TRUE, cex = 1.2,line = 0)
  } else {
  par(mfrow = c(2,2),oma = c(0, 0, 0, 0),mar = c(4,5,3.5,1))
  }

  minline = .5
  fdr = .01

  d = filter(df, method == 'all_all')
  ylm = makeplot(d,fdr = .01)
  title('search-all-asess-all (all-all)',line = minline + 1)
  title(bquote(~pi[0] ~'on'~ .(sum(!d$decoy))~'target PSMs: '
               ~ .(round(sum(d$decoy)/sum(!d$decoy)*100,1))~'%'),
        line = minline)
  title('a', line = minline,cex.main = 2,adj = 0)

  d = filter(df,method == 'sub_sub')
  makeplot(d,ylim = ylm,fdr = fdr)
  title('search-subset-asess-subset (sub-sub)',line = minline+1)
  title(bquote(~pi[0] ~'on'~ .(sum(!d$decoy))~'target PSMs: '
               ~ .(round(sum(d$decoy)/sum(!d$decoy)*100,1))~'%'),
        line = minline)
  title('b', line = minline,cex.main = 2,adj = 0)

  d = filter(df, method == 'all_sub')
  pi0 = sum(d$decoy)/sum(!d$decoy)
  ylm = makeplot(d,fdr = fdr,y_offset = TRUE)
  title('search-all-asess-subset (all-sub)',line = minline+1)
  title(bquote(~pi[0] ~'on'~ .(sum(!d$decoy))~'target subset PSMs: '
               ~ .(round(pi0*100,1))~'%'),
        line = minline)
  title('c', line = minline,cex.main = 2,adj = 0)

  d = group_by(df,spec_id) %>%
    mutate(subsetpsm = any(method == 'all_sub') & any(method == 'sub_sub')) %>% ungroup %>%
    filter(method == 'sub_sub')
  pi0 = sum(d$decoy)/sum(!d$decoy)
  d = filter(d,!decoy)
  ylm = makeplot(d,ylim = ylm,sbst = TRUE,fdr = fdr,y_offset = TRUE)
  title('search-subset-asess-subset (sub-sub)',line = minline+1)
  title(bquote(~pi[0] ~'on'~ .(sum(!d$decoy))~'target PSMs: '
               ~ .(round(pi0*100,1))~'%'),
        line = minline)
  title('d', line = minline,cex.main = 2,adj = 0)

    d = filter(df, database != 'complete')
    subset_name = first(d$database)
  subset_n = first(d$database_size)
  subset_complete_n = first(filter(df, database == 'complete')$database_size)

  add.box(df,fdr,ylm)
}

add.box = function(df,fdr,ylm){
  swapped = filter(df,method %in% c('all_all', 'sub_sub')) %>% group_by(spec_id) %>%
    filter(any((method == 'sub_sub') & !decoy & (length(unique(score)) > 1) & (FDR <= fdr))) %>%
    ungroup
  if (nrow(swapped) != 0){
    swapped_comp = filter(swapped, method == 'all_all')$score
    swapped_sub= filter(swapped, method == 'sub_sub')$score
    boxplot(list(swapped_comp,swapped_sub),border = c('black','orange'),#boxfill = c('black','orange'),
            add = TRUE,horizontal = TRUE,at = c(-ylm[2]*.05,-ylm[2]*.14),
            boxlwd = ylm[2] *.0003,boxwex = ylm[2] *.06,axes = FALSE)

    ## segments(min(swapped_sub),-ylm[2]*.155,min(swapped_sub),0,lwd = 4,lty = 3,lend = 1)

    text(55,-ylm[2]*.055,
         labels = 'score complete search' ,pos = 4,xpd = TRUE)
    text(55,-ylm[2]*.115,
         labels = 'score subset search' ,pos = 4,xpd = TRUE,col = 'dark orange')
  }}

makeplot= function(d,binWidth = 2,fdr = .01,
                   sbst = FALSE,xlim = c(0,80),ylim = NA,y_offset = FALSE,cex = 1
                   ){
  tar = filter(d,!decoy)$score

  cutoff = min(filter(d,!decoy,FDR <= fdr)$score)
  sub = filter(d,!decoy,subset,FDR <= fdr)

  if (sbst){
  dec = filter(d,!subsetpsm)$score
  col1 ="dark orange"
  col2 = "dark green"
  col3 = 'dark gray'
  fdrlabel = paste0(fdr*100,'% FDR\n',sum(!sub$subsetpsm),
                   ' (',round(sum(!sub$subsetpsm)/nrow(sub)*100,1),'%) / ',
                   sum(sub$subsetpsm),' of ',nrow(sub),
                   ' subset PSMs\nhave different / same sequence')
  lgnd = c('Different sequence as complete search',
           'Same sequence as complete search',
           'All PSMs')
  } else {
  dec = filter(d,decoy)$score
  col1 = "blue"
  col2 = "red"
  col3 = 'dark gray'
  fdrlabel =paste0(fdr*100,'% FDR\n',nrow(sub),' subset PSMs')
  lgnd = c('False PSMs','Correct PSMs','All PSMs')
  }

  breaks = seq(0,ceiling(xlim[2])+1,binWidth)

  histTar=hist(tar,breaks=breaks,plot = FALSE)
  histDec=hist(dec,breaks=breaks,plot = FALSE)
  if (is.na(dplyr::first(ylim))) {ylim = c(0,max(histTar$counts))} else { ylim}
  if (y_offset){ylim[1] = -ylim[2]*.13 }
  histTar2 = histTar
  histTar2$counts = ifelse(histTar2$counts > ylim[2]*1.15,ylim[2]*1.15,histTar2$counts)
  plot(histTar2,xlab="MS-GF+ Score",ylab = '',
         border = 'gray87',col = 'gray87',
       ylim = c(ylim[1],ylim[2]*1.2),main = '',axes = FALSE,font.lab = 2)
  axis(1,c(seq(0,xlim[2],length.out = 3),cutoff),sprintf('%.1f',c(seq(0,xlim[2],length.out = 3),cutoff)))
  axis(2,floor(seq(0,ylim[2],length.out = 3)),las = 1)
  text( - 4, ylim[2]*1.07, adj = 0.8, labels = "PSMs count", xpd = TRUE,font=2)

  counts = ifelse(histDec$counts > ylim[2]*1.15,ylim[2]*1.15,histDec$counts)
  points(histTar$mids-(.2*binWidth),counts,type="h",
         col=col1,lty=1,lwd=3,lend = 1)

  counts=histTar$counts-histDec$counts
  counts[counts<0]=0
  counts = ifelse(counts > ylim[2],ylim[2],counts)
  points(histTar$mids+(.2*binWidth),counts,type = "h",col=col2,lty=1,lwd=3,lend = 1)
  legend("bottomright",lgnd, text.col = c(col1,col2,col3),box.col = 0,inset = c(0,.55))

  segments(cutoff,0,cutoff,ylim[2]*1.15,lwd = 4,lty = 3,lend = 1)
  text(cutoff,ylim[2]*1.05,
       labels = fdrlabel,pos = 4,xpd = FALSE)
  tohigh = range(histTar$mids[histTar$counts > ylim[2]*1.15])

  segments(tohigh[1]-binWidth*1.2,ylim[2]*1.15,tohigh[2]+binWidth*1.2,ylim[2]*1.15,lty = 2,lwd = 1.5)
  ylim
}

make_overview_table = function(dat){
  out = map_df(dat,function(d){
    summarise(d,
              db = str_replace_all(gsub('(^.*)_\\[.*\\]_.*$','\\1', last(database)),'_',' '),
              db_size = last(database_size),
              complete_search_total = sum((database == 'complete') & subset),
              complete_search_targets = sum((database == 'complete') & subset & !decoy),
              complete_search_decoys= sum((database == 'complete') & subset & decoy),
              subset_search_total = sum(!(database == 'complete')),
              subset_search_targets = sum(!(database == 'complete') & !decoy),
              subset_search_decoys= sum(!(database == 'complete') & decoy))}) %>%
    arrange(desc(db_size))

  d = first(dat) %>% filter(database == 'complete')
  summarise(d,
            db = first(database),
            db_size = last(database_size),
            complete_search_total = n(),
            complete_search_targets = sum(!decoy),
            complete_search_decoys= sum(decoy),
            subset_search_total = NA,
            subset_search_targets = NA,
            subset_search_decoys = NA) %>%
    bind_rows(out)
}