R/IsotopicDistributionHDX.R

In OrgMassSpecR: Organic Mass Spectrometry

IsotopicDistributionHDX <- function(sequence, incorp, charge = 1, 
                                    custom = list(code = NULL, elements = NULL)) {

    if(length(custom$elements != 0)) {
        custom_elements <- c(C = 0, H = 0, N = 0, O = 0, S = 0, P = 0)
        custom_elements[names(custom$elements)] <- custom$elements
    }

    if(charge < 0 | charge > 3) stop("charge must be 1, 2, or 3")

    seq_vector <- strsplit(sequence, split = "")[[1]]
    x <- c(C = 0, H = 0, N = 0, O = 0, S = 0, P = 0)

    for(i in 1:(length(seq_vector))) {
        if(seq_vector[i] == "A") x <- x + c(C = 3, H = 5, N = 1, O = 1, S = 0, P = 0)
        if(seq_vector[i] == "R") x <- x + c(C = 6, H =12, N = 4, O = 1, S = 0, P = 0)
        if(seq_vector[i] == "N") x <- x + c(C = 4, H = 6, N = 2, O = 2, S = 0, P = 0)
        if(seq_vector[i] == "D") x <- x + c(C = 4, H = 5, N = 1, O = 3, S = 0, P = 0)
        if(seq_vector[i] == "C") x <- x + c(C = 3, H = 5, N = 1, O = 1, S = 1, P = 0) 
        if(seq_vector[i] == "E") x <- x + c(C = 5, H = 7, N = 1, O = 3, S = 0, P = 0)
        if(seq_vector[i] == "Q") x <- x + c(C = 5, H = 8, N = 2, O = 2, S = 0, P = 0)
        if(seq_vector[i] == "G") x <- x + c(C = 2, H = 3, N = 1, O = 1, S = 0, P = 0)
        if(seq_vector[i] == "H") x <- x + c(C = 6, H = 7, N = 3, O = 1, S = 0, P = 0)
        if(seq_vector[i] == "I") x <- x + c(C = 6, H =11, N = 1, O = 1, S = 0, P = 0)
        if(seq_vector[i] == "L") x <- x + c(C = 6, H =11, N = 1, O = 1, S = 0, P = 0)
        if(seq_vector[i] == "K") x <- x + c(C = 6, H =12, N = 2, O = 1, S = 0, P = 0)
        if(seq_vector[i] == "M") x <- x + c(C = 5, H = 9, N = 1, O = 1, S = 1, P = 0)
        if(seq_vector[i] == "F") x <- x + c(C = 9, H = 9, N = 1, O = 1, S = 0, P = 0)
        if(seq_vector[i] == "P") x <- x + c(C = 5, H = 7, N = 1, O = 1, S = 0, P = 0)
        if(seq_vector[i] == "S") x <- x + c(C = 3, H = 5, N = 1, O = 2, S = 0, P = 0)
        if(seq_vector[i] == "T") x <- x + c(C = 4, H = 7, N = 1, O = 2, S = 0, P = 0)
        if(seq_vector[i] == "W") x <- x + c(C =11, H =10, N = 2, O = 1, S = 0, P = 0)
        if(seq_vector[i] == "Y") x <- x + c(C = 9, H = 9, N = 1, O = 2, S = 0, P = 0)
        if(seq_vector[i] == "V") x <- x + c(C = 5, H = 9, N = 1, O = 1, S = 0, P = 0)
            
        if(length(custom$elements != 0))
            if(seq_vector[i] == custom$code) x <- x + custom_elements    
    }

    ## add N-terminal H and C-terminal OH
    elements <- x + c(C = 0, H = 2, N = 0, O = 1, S = 0, P = 0) 

    num_exch_sites <- ExchangeableAmides(sequence)

    simulation <- function(elements) {
						
        mz <- vector(mode="numeric")
		
        ## mass of carbons
        mc <- sum(sample(c(12.0000000, 13.0033548378), 
                  size = elements["C"], 
                  replace = TRUE, 
                  prob = c(0.9893, 0.0107)))   
    
        ## mass of unexchanged (non-backbone) hydrogens
        mh <- sum(sample(c(1.0078250321, 2.0141017780), 
                  size = elements["H"] - num_exch_sites,
                  replace = TRUE, 
                  prob = c(0.999885, 0.000115)))
                  
        ## mass of exchanged (backbone) hydrogens
        md <- sum(sample(c(1.0078250321, 2.0141017780),
                  size = num_exch_sites,
                  replace = TRUE,
                  prob = c(1 - incorp, incorp)))

        ## mass of nitrogens
        mn <- sum(sample(c(14.0030740052, 15.0001088984), 
                  size = elements["N"], 
                  replace = TRUE, 
                  prob = c(0.99632, 0.00368)))

        ## mass of oxygens
        mo <- sum(sample(c(15.9949146221, 16.99913150, 17.9991604), 
                  size = elements["O"], 
                  replace = TRUE, 
                  prob=c(0.99757, 0.00038, 0.00205)))

        ## mass of sulfers
        ms <- sum(sample(c(31.97207069, 32.97145850, 33.96786683, 35.96708088), 
                  size=elements["S"], 
                  replace = TRUE, 
                  prob = c(0.9493, 0.0076, 0.0429, 0.0002)))

        ## mass of charge
        mch <- sum(sample(c(1.0072764522, 2.0135531981), 
                   size = charge, 
                   replace = TRUE, 
                   prob=c(0.999885, 0.000115)))

        ## m/z of molecule
        if(charge != 0) mz <- sum(mc, mh, md, mn, mo, ms, mch) / charge 
            else mz <- mz <- sum(mc, mh, md, mn, mo, ms, mch)		
        return(mz)	
    }

    ## run simulation
    sim <- replicate(10000, expr = simulation(elements))  
	
    ## bin ions
    b <- seq(from = min(sim) - (1 / (2 * charge)),
             to = max(sim) + 1, 
             by = 1 / charge)
    bins <- cut(sim, breaks = b)
    mz <- round(tapply(sim, bins, mean), digits = 2)
    intensity <- as.vector(table(bins))
    spec <- data.frame(mz, intensity)
    spec <- subset(spec, intensity != 0)
    spec <- transform(spec, percent = round(intensity / max(intensity) * 100, digits = 2))
    row.names(spec) <- 1:(nrow(spec))

    return(spec)

}