# R program for comparing mass spectra

msbin <- function(x, y, bin.width=0.25, x.max=max(x), x.min=min(x)){
# ---------------------------------------------------------------
# Simple binning of an MS/MS spectrum
#
# x             input array of m/z values
# y             input array of intensity values
# bin.width     m/z bin width of output
#
# Value:
# x     =       evenly sampled m/z values for the binned spectrum
# y     =       binned intensity values
# ---------------------------------------------------------------

  x.n <- length(x)

  n.bins <- round( (x.max - x.min) / bin.width ) + 1

  x.bin <- seq(from=x.min, to=x.max, length=n.bins)
  y.bin <- rep(0, n.bins)

  for(i in 1:x.n){
    idx <- round( (x[i] - x.min) / bin.width ) + 1
    y.bin[idx] <- y.bin[idx] + y[i]
  }

  out <- list(x = x.bin, y = y.bin)
  out
}

fftconvolve <- function(y1, y2){
# ----------------------------------------------------------
# Function to perform FFT based convolution of two sequences
# ----------------------------------------------------------

  y1.n <- length(y1)
  y2.n <- length(y2)
 
  n.out <- y1.n + y2.n - 1

  nf <- log(n.out) / log(2.0)
  if (nf == floor(nf)){
        n <- floor(nf)
  }
  else{
        n <- floor(nf) + 1
  }

  n.round <- 2^n

  if (n.round > y1.n){   
    y1 <- c(y1, rep(0.0, n.round - y1.n))
  }
  if (n.round > y2.n){
        y2 <- c(y2, rep(0.0, n.round - y2.n))
  }

  y1.fft <- fft(y1)
  y2.fft <- fft(y2)
     
  y.real <- (Re(y1.fft)*Re(y2.fft) - Im(y1.fft)*Im(y2.fft)) / n.round
  y.imag <- (Im(y1.fft)*Re(y2.fft) + Re(y1.fft)*Im(y2.fft)) / n.round
  y.ifft <- complex(length=n.round, real=y.real, imaginary=y.imag)
 
  y.out <- fft(y.ifft, inverse=T)
 
  out <- Re(y.out[1:n.out])
}

gmg <- function(x, a0, a1, a2, a3){
# -----------------------------------------------------
# Produce a Gaussian Modified Gaussian peak given x and
# a0 (peak area)
# a1 (mean)
# a2 (sigma)
# a3 (skewness)
# -----------------------------------------------------
  y <- a0 * exp(-0.5 * (x-a1)^2/(a3^2 + a2^2)) *
                2*pnorm(a3*(x-a1)/(a2*sqrt(a3^2 + a2^2))) /
                (sqrt(2 * 3.14159) * sqrt(a3^2 + a2^2))
  y
}

peak <- gmg(seq(-5,5,0.1),3,1,.2,0)

dmmp <- as.matrix(read.table("dmmp.dat", sep=","))
dmmp.bin <- msbin( dmmp[,1], dmmp[,2], bin.width=0.1, x.min=20, x.max=200 )
dmmp.m<-fftconvolve(dmmp.bin$y, peak)
dmmp.mz<-data.frame(dmmp.bin$x,dmmp.m[61:(length(dmmp.m)-40)])
names(dmmp.mz)<-c("mz","I")
names(dmmp.bin)<-c("mz","I")

demp <- as.matrix(read.table("demp.dat", sep=","))
demp.bin <- msbin( demp[,1], demp[,2], bin.width=0.1, x.min=20, x.max=200 )
demp.m<-fftconvolve(demp.bin$y, peak)
demp.mz<-data.frame(demp.bin$x,demp.m[61:(length(demp.m)-40)])
names(demp.mz)<-c("mz","I")
names(demp.bin)<-c("mz","I")

dimp <- as.matrix(read.table("dimp.dat", sep=","))
dimp.bin <- msbin( dimp[,1], dimp[,2], bin.width=0.1, x.min=20, x.max=200 )
dimp.m<-fftconvolve(dimp.bin$y, peak)
dimp.mz<-data.frame(dimp.bin$x,dimp.m[61:(length(dimp.m)-40)])
names(dimp.mz)<-c("mz","I")
names(dimp.bin)<-c("mz","I")

must <- as.matrix(read.table("must.dat", sep=","))
must.bin <- msbin( must[,1], must[,2], bin.width=0.1, x.min=20, x.max=200 )
must.m<-fftconvolve(must.bin$y, peak)
must.mz<-data.frame(must.bin$x,must.m[61:(length(must.m)-40)])
names(must.mz)<-c("mz","I")
names(must.bin)<-c("mz","I")

par(mfcol=c(4,1))
xr <- c(20,200)

plot(dmmp.mz, xlim=xr, type="l")
lines(dmmp, type="h", col="red")
title("DMMP [1]")

plot(demp.mz, xlim=xr,type="l")
lines(demp, type="h", col="red")
title("DEMP [2]")

plot(dimp.mz, xlim=xr, type="l")
lines(dimp, type="h", col="red")
title("DIMP [3]")

plot(must.mz, xlim=xr, type="l")
lines(must, type="h", col="red")
title("Mustard [4]")

dist<-matrix(rep(0,16),ncol=4)

dist[2,1]=acos(dmmp.mz$I %*% demp.mz$I/ sqrt(sum(demp.mz$I^2)*sum(dmmp.mz$I^2)))*2*pi
dist[3,1]=acos(dmmp.mz$I %*% dimp.mz$I/ sqrt(sum(dimp.mz$I^2)*sum(dmmp.mz$I^2)))*2*pi
dist[4,1]=acos(must.mz$I %*% demp.mz$I/ sqrt(sum(must.mz$I^2)*sum(dmmp.mz$I^2)))*2*pi
dist[3,2]=acos(demp.mz$I %*% dimp.mz$I/ sqrt(sum(demp.mz$I^2)*sum(dimp.mz$I^2)))*2*pi
dist[4,2]=acos(demp.mz$I %*% must.mz$I/ sqrt(sum(demp.mz$I^2)*sum(must.mz$I^2)))*2*pi
dist[4,3]=acos(must.mz$I %*% dimp.mz$I/ sqrt(sum(must.mz$I^2)*sum(dimp.mz$I^2)))*2*pi

dist.bin<-matrix(rep(0,16),ncol=4)

dist.bin[2,1]=acos(dmmp.bin$I %*% demp.bin$I/ sqrt(sum(demp.bin$I^2)*sum(dmmp.bin$I^2)))*2*pi
dist.bin[3,1]=acos(dmmp.bin$I %*% dimp.bin$I/ sqrt(sum(dimp.bin$I^2)*sum(dmmp.bin$I^2)))*2*pi
dist.bin[4,1]=acos(must.bin$I %*% demp.bin$I/ sqrt(sum(must.bin$I^2)*sum(dmmp.bin$I^2)))*2*pi
dist.bin[3,2]=acos(demp.bin$I %*% dimp.bin$I/ sqrt(sum(demp.bin$I^2)*sum(dimp.bin$I^2)))*2*pi
dist.bin[4,2]=acos(demp.bin$I %*% must.bin$I/ sqrt(sum(demp.bin$I^2)*sum(must.bin$I^2)))*2*pi
dist.bin[4,3]=acos(must.bin$I %*% dimp.bin$I/ sqrt(sum(must.bin$I^2)*sum(dimp.bin$I^2)))*2*pi