Simulate_Mass_Spectrum

Script

#
# calculate isotopic distributions of molecules using the FFT
#
# (c) Magnus Palmblad, 1999
#
# translated from MATLAB to R (by Magnus Palmblad) 2010
#

MAX_ELEMENTS <- 5; MAX_MASS <- 2^13

M <- elemental_composition       # read elemental composition (empirical formula)
M[1]=M[1]+1;                     # add H+

mi_mass <- sum(c(1.0078250321,12,14.0030740052,15.9949146221,31.97207069)*M)
int_mass <- sum(c(1,12,14,16,32)*M)

A<-rep(0,times=MAX_ELEMENTS*MAX_MASS); dim(A)<-c(MAX_ELEMENTS,MAX_MASS)
A[1,2:3]   <- c(0.9998443,0.0001557)
A[2,13:14] <- c(0.98889,0.01111)
A[3,15:16] <- c(0.99634,0.00366)
A[4,17:19] <- c(0.997628,0.000372,0.002000)
A[5,33:37] <- c(0.95018,0.00750,0.04215,0,0.00017)

tA <- mvfft(t(A))                                  # FFT along each element's isotopic distribution

ptA <- t(rep(1,MAX_MASS))
for(i in 1:MAX_ELEMENTS) ptA <- ptA*(tA[,i]^M[i])  # multiply transforms (elementwise)

riptA <- Re(mvfft(t(ptA), inverse = TRUE))         # inverse FFT to get convolutions

distribution <- riptA/sum(riptA)

before<-round(int_mass/1000)+1                     # buffer distribution according to molecular mass
after<-round(int_mass/700)+5

main_distribution <- distribution[(which.max(riptA)-before):(which.max(riptA)+after)]
mz_axis <- (mi_mass/int_mass)*c((which.max(riptA)-before):(which.max(riptA)+after)); 

FWHM<-mi_mass/RP
stddev<-FWHM/2*sqrt(2*log(2))

spectrum_x <- (mi_mass/int_mass)*seq((which.max(riptA)-before),(which.max(riptA)+after),by=0.01)
spectrum_y <- rep(0,times=length(spectrum_x))
for(i in 1:length(mz_axis)) {spectrum_y[i*100-99]<-main_distribution[i]; if((i>1)&&(i

R Server

Script

int nH=2, nC, nN, nO=1, nS; // start with water
char [] sequence = peptide_sequence.toCharArray();
 
for (int i=0; i