Build XCalibur Inclusion List

build_inclusion_list.R

# BUILD THERMO XCALIBUR INCLUSION LIST

# returns inclusion list dataframe
build_inclusion_list <- function(
  input,                        # dataframe, with columns: 
                                # `Raw file`, `Modified sequence`, `Retention time`, `Charge`
  file='inclusion_list.csv',    # output file
  rt_width=3,                   # width, in minutes, of the RT window on either side
  mz_decimals=5,                # number of decimals to round m/z to
  rt_decimals=1,                # number of decimals to round RT to
  transform_rt=NA,              # number (RT = RT + x) or function ( RT = f(RT) ) to transform RTs
  rt_alignment=FALSE,           # perform RT alignment? (simple linear model)
  rt_alignment_reference=NA,    # name of the raw file to align to (reference experiment). 
                                # if NA, then use last experiment as reference.
  TMT=TRUE,                     # is this a TMT experiment?
  polarity='Positive'           # positive or negative ionization mode?
) {
  
  library(tidyverse)
  library(stringr)
  
  # element monoisotopic masses
  element_mass <- c(
    H=1.0078250321,
    D=2.014101779,
    C=12.0,
    Cx=13.0033548378,
    N=14.0030740052,
    Nx=15.0001088984,
    O=15.9949146221,
    Ox=17.9991603,
    F=18.99840320,
    P=30.97376151,
    S=31.97207069,
    Cl=34.96885271,
    Na=22.98976967,
    K=38.9637069,
    Ca=39.9625912,
    Fe=53.9396148
  )
  
  # amino acid monoisotopic masses
  aa_mass = c(
    A=71.0371137878,
    C=103.0091844778,
    D=115.026943032,
    E=129.0425930962,
    F=147.0684139162,
    G=57.0214637236,
    H=137.0589118624,
    I=113.0840639804,
    K=128.0949630177,
    L=113.0840639804,
    M=131.0404846062,
    N=114.0429274472,
    P=97.052763852,
    Q=128.0585775114,
    R=156.1011110281,
    S=87.0320284099,
    T=101.0476784741,
    U=150.95363, # Selenocysteine
    V=99.0684139162,
    W=186.0793129535,
    Y=163.0633285383,
    X=113.0840639804 # Xle - Isoleucine or Leucine
  )
  
  mods <- c(
    `p`=79.9663304084, # phosphorylation (HPO3)
    `_(ac)`=42.0105646863, # n-terminal acetylation (C2H2O)
    `(ox)`=15.9949146221, # methionine oxidation (O)
    `(ar)`=10.0082685996, # superheavy arginine (Cx6Nx4) - (C6N4)
    `(ly)`=8.01419881320 # superheavy lysine (Cx6Nx2) - (C6N2)
  )
  
  proton_mass = 1.0072764666
  electron_mass = 0.0005485799
  
  water_mass = 18.0105646863
  tmt_mass = 229.162932141
  

# Perform RT alignment? ---------------------------------------------------
  
  input$RT_aligned <- input$`Retention time`
  
  if(rt_alignment) {
    
    # if no reference experiment given, use the last one as the reference
    ref_exp_name <- NA
    if(is.na(rt_alignment_reference)) {
      ref_exp_name <- last(sort(unique(input$`Raw file`)))
    } else {
      ref_exp_name <- rt_alignment_reference
    }
    # pull reference experiment data from input dataframe
    ref_exp <- input %>%
      filter(grepl(ref_exp_name, `Raw file`)) %>%
      distinct(`Raw file`, `Modified sequence`, .keep_all=T) %>%
      arrange(`Modified sequence`)
    
    input_aligned <- data.frame()
    for(exp in exps) {
      
      e <- input %>%
        filter(`Raw file` == exp) %>%
        distinct(`Raw file`, `Modified sequence`, .keep_all=T) %>%
        arrange(`Modified sequence`)
      
      if(exp != ref_exp_name) {
        # find common peptides
        common_peps <- intersect(ref_exp$`Modified sequence`, e$`Modified sequence`)
        
        # linear alignment
        y <- e %>% filter(`Modified sequence` %in% common_peps) %>% pull(`Retention time`)
        x <- ref_exp %>% filter(`Modified sequence` %in% common_peps) %>% pull(`Retention time`)
        fit <- lm(y ~ x)
        
        e$RT_aligned <- (e$`Retention time` - fit$coefficients[1]) / fit$coefficients[2]
      }
      
      input_aligned <- rbind(input_aligned, e)
    }
    
    input <- input_aligned
  }
  

# Begin building inclusion list dataframe ---------------------------------
  
  inc_list <- input %>%
    mutate(Sequence=gsub('([^A-Z])', '', `Modified sequence`)) %>% # derive unmodified AA sequence
    group_by(`Modified sequence`) %>%
    summarise(RT=mean(RT_aligned), 
              #RT_spread=max(RT_aligned)-min(RT_aligned),
              Charge=first(Charge),
              Sequence=first(Sequence))

# Recalculate peptide mass ------------------------------------------------

  inc_list$Mass <- sapply(inc_list$`Modified sequence`, function(s) {
    mass <- 0.0
    
    # Add amino acid masses
    for(i in 1:nchar(s)) {
      char <- substr(s, i, i)
      if(char %in% names(aa_mass)) { mass <- mass + aa_mass[char] }
    }
    
    # Add modification masses
    for(mod in names(mods)) {
      mod_count <- str_count(s, mod)
      mass <- mass + (mod_count * mods[mod])
    }
    
    if(TMT) {
      # add a TMT for each lysine
      mass <- mass + (str_count(s, 'K') * tmt_mass)
      
      # add a TMT if the n-terminus is not acetylated
      if(!grepl('_(ac)', s)) {
        mass <- mass + tmt_mass
      }
    }
    
    # add a water
    mass <- mass + water_mass
    
  })
  
  # Calculate M+H and m/z
  if(polarity == 'Positive') {
    inc_list$`M+H` <- inc_list$Mass + (inc_list$Charge * proton_mass)
  }
  inc_list$`m/z` <- inc_list$`M+H` / inc_list$Charge
  

# Apply RT Transform? -----------------------------------------------------

  if(!is.na(transform_rt)) {
    if(class(transform_rt) == 'numeric' | class(transform_rt) == 'integer') {
      inc_list$RT <- inc_list$RT + transform_RT
    } else if (class(transform_rt) == 'function') {
      inc_list$RT <- transform_rt(inc_list$RT)
    } else {
      stop('Argument "transform_rt" has invalid type. "transform_rt" needs to be an integer, numeric, or a function f(RT) that takes in an RT and returns an RT.')
    }
  }
  

# Final Data Wrangling ----------------------------------------------------

  # collapse duplicate ions, wrangle columns
  inc_list <- inc_list %>%
    group_by(Sequence, `m/z`) %>%
    summarise(RT=mean(RT), Charge=first(Charge)) %>%
    ungroup() %>%
    mutate(
      `Mass [m/z]`=round(`m/z`, mz_decimals),
      `Formula [M]`=NA, `Formula type`=NA, `Species`=NA,
      `CS [z]`=Charge,
      `Polarity`=polarity,
      `Start [min]`=round(RT - rt_width, rt_decimals),
      `End [min]`=round(RT + rt_width, rt_decimals),
      `(N)CE`=NA, `(N)CE type`=NA, `MSX ID`=NA, `Comment`=NA
    ) %>%
    select(-c('Sequence', 'Charge', 'm/z', 'RT')) %>%
    arrange(`Start [min]`)

# Save file ---------------------------------------------------------------

  write_csv(inc_list, path=file, na='')  
  
  return(inc_list)
}

example_inclusion_list.csv