Calculate_Isotopic_Distribution

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)

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)

png(filename="distribution_image.png", height=402, width=650, bg="white")
plot(0:(MAX_MASS-1), distribution, type = "h", xlab="integer mass", ylab="frequency")   # plot of isotopic distribution
plot_from <- which.max(riptA)-10
plot_to <- which.max(riptA)+20
plot(plot_from:plot_to, riptA[(plot_from+1):(plot_to+1)]/sum(riptA[(plot_from+1):(plot_to+1)]), type = "h", xlab="integer mass", ylab="frequency")
dev.off()

R Server

Script

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