################################################################## ##### Point to the working directory containing the data ######### setwd("path") ################################################################## ################################################################## ################################################################## ############# load needed libraries and functions ################ library(xlsx) ## to open the spectra directly from the Excel file library(hyperSpec) ## other baseline modelling ## Function to range data between 0 and 1 range01 <- function(x){(x-min(x))/(max(x)-min(x))} ## Gaussian function gaussian <- function(x, peak.intensity, peak.pos, peak.fwhm){ peak.intensity * exp(- (4 *log(2)*(x-peak.pos)^2)/(peak.fwhm)^2) } ## increase the control of iterations in nls new.nls.control <- nls.control(maxiter = 100000, tol = 1e-05, minFactor = 1e-10, printEval = FALSE, warnOnly = FALSE) ################################################################## ################################################################## ################################################################## ################################################################## ### Fig. 1b, emission spectroscopy ################################################################## ################################################################## # read the spectra emission_spectra <- read.xlsx2("Toussaint-etal-spectra.xls", 1, colClasses = NA) # plot all, not scaled, with a vertical shift quartz() plot(x=emission_spectra[,1], y=emission_spectra[,2], t="l", xlab="Wavelength (nm)", ylab="Intensity", main="", ylim=c(0,1.6E6)) for(i in 3:10) { lines(x=emission_spectra[,1], y=((i-2)*1E4)+emission_spectra[,i]) } #### Fig 1b #### # plot all, ranged between 0 and 1, with a vertical shift quartz() plot(x=emission_spectra[,1], y=range01(emission_spectra[,2]), t="l", xlab="Wavelength (nm)", ylab="Norm. intensity", main="", ylim=c(0,9)) for(i in 3:10) { lines(x=emission_spectra[,1], y=((i-2))+range01(emission_spectra[,i])) } ################################################################## ################################################################## ### Fig. 2a, excitation spectroscopy ################################################################## ################################################################## # read the spectra excitation_spectra <- read.xlsx2("Toussaint-etal-spectra.xls", 2, colClasses = NA) ## plot all, not scaled, with a vertical shift quartz() plot(x=excitation_spectra[,1], y=excitation_spectra[,2], t="l", xlab="Wavelength (nm)", ylab="Intensity", main="", ylim=c(0,1.5E5)) for(i in 3:10) { lines(x=excitation_spectra[,1], y=((i-2)*1E4)+excitation_spectra[,i]) } #### Fig 2a #### ## plot all, ranged between 0 and 1, with a vertical shift quartz() plot(x=excitation_spectra[,1], y=range01(excitation_spectra[,2]), t="l", xlab="Wavelength (nm)", ylab="Norm. intensity", main="", ylim=c(0,4)) for(i in 3:10) { lines(x=excitation_spectra[,1], y=((i-2)*0.35)+range01(excitation_spectra[,i])) } ### add the position of Q bands for reference spectra ## (arbitrary intensities) ## protoporphyrin IX (after our protoporphyrin IX reference spectrum) for(i in 2:10) { segments(x0=c(505, 539, 575, 630), y0=c(0, 0, 0, 0)+((i-2)*0.35), x1 = c(505, 539, 575, 630), y1 = c(0.3, 0.2, 0.2, 0.2)+((i-2)*0.35), col = "red", lwd=1) } ## uroporphyrin I or III (DiNello & Chang 1978) for(i in 2:10) { segments(x0=c(502, 536, 572, 627), y0=c(0, 0, 0, 0)+((i-2)*0.35), x1 = c(502, 536, 572, 627), y1 = c(0.2, 0.2, 0.2, 0.2)+((i-2)*0.35), col = "green", lwd=1) } ## porphyrin c (DiNello & Chang 1978) for(i in 2:10) { segments(x0=c(504, 538, 570, 625.6), y0=c(0, 0, 0, 0)+((i-2)*0.35), x1 = c(504, 538, 570, 625.6), y1 = c(0.2, 0.2, 0.2, 0.2)+((i-2)*0.35), col = "darkgreen", lwd=1) } ## porphyrin S-411 (DiNello & Chang 1978) for(i in 2:10) { segments(x0=c(505, 553, 578, 640), y0=c(0, 0, 0, 0)+((i-2)*0.35), x1 = c(505, 553, 578, 640), y1 = c(0.2, 0.2, 0.2, 0.2)+((i-2)*0.35), col = "blue", lwd=1) } ## add the position of the Soret band for the first spectrum (hedgehog) abline(v=410, col="burlywood", lty=2) legend("topright",legend=c("protoporphyrin IX","uroporphyrin I or III","porphyrin c","porphyrin S-411"),lwd=c(1.5,1.5,1.5,1.5),lty=c(1,1,1,1), col=c("red","green","darkgreen","blue")) ################################################################## ################################################################## ### Fig. 2a, assessing the position of Q bands by decomposing ### the Q-bands domain of the spectrum into Gaussian ################################################################## ################################################################## ## To ensure the good success of the Gaussian decomposition, we apply a baseline correction to the Q-bands domain of the spectrum (As such, or Gaussian will fit the Q-bands and not a backgrounf contribution.) ## Define the Q-bands domain, which is 480-655 nm = excitation_spectra[121:296,1] Qbands_domain <- excitation_spectra[121:296,1] ## Apply a baseline subtraction in this domain by (i) creating a hyperSpec Object from a Data Matrix, which is fitted by (ii) a rubberband baseline method with a bending for(i in 2:10) { spc <- new ("hyperSpec", spc=t(as.matrix(range01(excitation_spectra[121:296,i]))), wavelength=as.vector(Qbands_domain)) bend <- 1 * wl.eval (spc, function(x) x^2, normalize.wl=normalize01) bl_b <- spc.rubberband (spc + bend) - bend assign(paste("blsub_V", i, sep = ""), range01(as.vector(spc$spc-bl_b$spc))) } # plot the baseline-subtracted spectra to check that everything went fine. quartz() plot(x=Qbands_domain, y=blsub_V2, t="l", xlab="Wavelength (nm)", ylab="Norm. intensity", main="", ylim=c(0,9)) for(i in 3:10) { name <- paste("blsub_V", i, sep = "") lines(x=Qbands_domain, y=((i-2)+get(name))) } ####### ### Gaussian decomposition of Erinaceus europaeus NRM_594242 ## V10 ####### # 4 gaussian ## give starting Gaussian parameters for the decomposition (nls model) gau1_int <- 1 gau1_pos <- 505 gau1_fwhm <- 20 gau2_int <- 1 gau2_pos <- 539 gau2_fwhm <- 20 gau3_int <- 1 gau3_pos <- 577 gau3_fwhm <- 20 gau4_int <- 1 gau4_pos <- 630 gau4_fwhm <- 20 ### nls V10_nls4Qbands <- nls(blsub_V10 ~ gaussian(Qbands_domain, gau1.int, gau1.pos, gau1.fwhm) + gaussian(Qbands_domain, gau2.int, gau2.pos, gau2.fwhm) + gaussian(Qbands_domain, gau3.int, gau3.pos, gau3.fwhm) + gaussian(Qbands_domain, gau4.int, gau4.pos, gau4.fwhm), start = list(gau1.int = gau1_int, gau1.pos = gau1_pos, gau1.fwhm = gau1_fwhm, gau2.int = gau2_int, gau2.pos = gau2_pos, gau2.fwhm = gau2_fwhm, gau3.int = gau3_int, gau3.pos = gau3_pos, gau3.fwhm = gau3_fwhm, gau4.int = gau4_int, gau4.pos = gau4_pos, gau4.fwhm = gau4_fwhm), control=new.nls.control) # V10_nls4Qbands # gau1.int gau1.pos gau1.fwhm gau2.int gau2.pos gau2.fwhm gau3.int gau3.pos gau3.fwhm # 1.0107 504.5375 20.2012 0.7230 539.5708 13.7989 0.6401 578.4509 20.4872 # gau4.int gau4.pos gau4.fwhm # 0.4680 629.4065 14.5041 # residual sum-of-squares: 0.1895 ### predict the result of the model predict_V10_nls4Qbands <- predict(V10_nls4Qbands,list(x=Qbands_domain)) ### plot the result of the model along with the spectrum and baseline quartz() plot(x=Qbands_domain, y=blsub_V10, t="l", xlab="wavelength (nm)", ylab="Norm. intensity (a.u.)") lines(x=Qbands_domain, y=predict_V10_nls4Qbands, col="red") ## not bad at all! ## add the Guassians lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V10_nls4Qbands)[1], coef(V10_nls4Qbands)[2], coef(V10_nls4Qbands)[3]), col="blue") lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V10_nls4Qbands)[4], coef(V10_nls4Qbands)[5], coef(V10_nls4Qbands)[6]), col="darkblue") lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V10_nls4Qbands)[7], coef(V10_nls4Qbands)[8], coef(V10_nls4Qbands)[9]), col="dodgerblue") lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V10_nls4Qbands)[10], coef(V10_nls4Qbands)[11], coef(V10_nls4Qbands)[12]), col="lightblue") ####### ### Gaussian decomposition of Monodelphis brevicaudata JAGUARS-M2838 ## V7 ####### # 5 gaussian ## As most starting Gaussian parameters can be reused from the previous starting values, it is needed here to at least have 5 different Gaussian positions. Done here by "splitting Guassian 4" into two distinct Gaussian positions 4a and 4b. This is anyway just starting values for the decomposition to model 5 instead of 4 Gaussian gau4a_pos <- 620 gau4b_pos <- 640 ### nls V7_nls5Qbands <- nls(blsub_V7 ~ gaussian(Qbands_domain, gau1.int, gau1.pos, gau1.fwhm) + gaussian(Qbands_domain, gau2.int, gau2.pos, gau2.fwhm) + gaussian(Qbands_domain, gau3.int, gau3.pos, gau3.fwhm) + gaussian(Qbands_domain, gau4a.int, gau4a.pos, gau4a.fwhm) + gaussian(Qbands_domain, gau4b.int, gau4b.pos, gau4b.fwhm), start = list(gau1.int = gau1_int, gau1.pos = gau1_pos, gau1.fwhm = gau1_fwhm, gau2.int = gau2_int, gau2.pos = gau2_pos, gau2.fwhm = gau2_fwhm, gau3.int = gau3_int, gau3.pos = gau3_pos, gau3.fwhm = gau3_fwhm, gau4a.int = gau4_int, gau4a.pos = gau4a_pos, gau4a.fwhm = gau4_fwhm, gau4b.int = gau4_int, gau4b.pos = gau4b_pos, gau4b.fwhm = gau4_fwhm), control=new.nls.control) # V7_nls5Qbands # gau1.int gau1.pos gau1.fwhm gau2.int gau2.pos gau2.fwhm gau3.int gau3.pos # 0.9860 501.4081 20.8976 0.5161 535.3957 12.9000 0.5090 570.5935 # gau3.fwhm gau4a.int gau4a.pos gau4a.fwhm gau4b.int gau4b.pos gau4b.fwhm # 22.1583 0.4529 620.0669 15.8316 0.5390 638.3423 17.2826 # residual sum-of-squares: 0.1543 ### predict the result of the model predict_V7_nls5Qbands <- predict(V7_nls5Qbands,list(x=Qbands_domain)) ### plot the result of the model along with the spectrum and baseline quartz() plot(x=Qbands_domain, y=blsub_V7, t="l", xlab="wavelength (nm)", ylab="Norm. intensity (a.u.)") lines(x=Qbands_domain, y=predict_V7_nls5Qbands, col="red") ## not bad at all! ## add the Gaussian lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V7_nls5Qbands)[1], coef(V7_nls5Qbands)[2], coef(V7_nls5Qbands)[3]), col="blue") lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V7_nls5Qbands)[4], coef(V7_nls5Qbands)[5], coef(V7_nls5Qbands)[6]), col="darkblue") lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V7_nls5Qbands)[7], coef(V7_nls5Qbands)[8], coef(V7_nls5Qbands)[9]), col="dodgerblue") lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V7_nls5Qbands)[10], coef(V7_nls5Qbands)[11], coef(V7_nls5Qbands)[12]), col="lightblue") lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V7_nls5Qbands)[13], coef(V7_nls5Qbands)[14], coef(V7_nls5Qbands)[15]), col="darkslategray3") ####### ### Gaussian decomposition of Marmosa murina JAGUARS-M535 ## V6 ####### # 5 gaussian ### nls V6_nls5Qbands <- nls(blsub_V6 ~ gaussian(Qbands_domain, gau1.int, gau1.pos, gau1.fwhm) + gaussian(Qbands_domain, gau2.int, gau2.pos, gau2.fwhm) + gaussian(Qbands_domain, gau3.int, gau3.pos, gau3.fwhm) + gaussian(Qbands_domain, gau4a.int, gau4a.pos, gau4a.fwhm) + gaussian(Qbands_domain, gau4b.int, gau4b.pos, gau4b.fwhm), start = list(gau1.int = gau1_int, gau1.pos = gau1_pos, gau1.fwhm = gau1_fwhm, gau2.int = gau2_int, gau2.pos = gau2_pos, gau2.fwhm = gau2_fwhm, gau3.int = gau3_int, gau3.pos = gau3_pos, gau3.fwhm = gau3_fwhm, gau4a.int = gau4_int, gau4a.pos = gau4a_pos, gau4a.fwhm = gau4_fwhm, gau4b.int = gau4_int, gau4b.pos = gau4b_pos, gau4b.fwhm = gau4_fwhm), control=new.nls.control) # V6_nls5Qbands # gau1.int gau1.pos gau1.fwhm gau2.int gau2.pos gau2.fwhm gau3.int gau3.pos # 0.6921 502.1731 19.5834 0.3688 536.0730 12.4375 0.6680 574.1728 # gau3.fwhm gau4a.int gau4a.pos gau4a.fwhm gau4b.int gau4b.pos gau4b.fwhm # 35.5587 0.6386 624.7477 30.1555 0.6439 639.5688 16.2770 # residual sum-of-squares: 0.2593 ### predict the result of the model predict_V6_nls5Qbands <- predict(V6_nls5Qbands,list(x=Qbands_domain)) ### plot the result of the model along with the spectrum and baseline quartz() plot(x=Qbands_domain, y=blsub_V6, t="l", xlab="wavelength (nm)", ylab="Norm. intensity (a.u.)") lines(x=Qbands_domain, y=predict_V6_nls5Qbands, col="red") ## not bad at all! ## add the Gaussian lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V6_nls5Qbands)[1], coef(V6_nls5Qbands)[2], coef(V6_nls5Qbands)[3]), col="blue") lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V6_nls5Qbands)[4], coef(V6_nls5Qbands)[5], coef(V6_nls5Qbands)[6]), col="darkblue") lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V6_nls5Qbands)[7], coef(V6_nls5Qbands)[8], coef(V6_nls5Qbands)[9]), col="dodgerblue") lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V6_nls5Qbands)[10], coef(V6_nls5Qbands)[11], coef(V6_nls5Qbands)[12]), col="lightblue") lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V6_nls5Qbands)[13], coef(V6_nls5Qbands)[14], coef(V6_nls5Qbands)[15]), col="darkslategray3") ####### ### Gaussian decomposition of Glaucomys volans ZMB_Mam_60634 (ventral) ## V5 ####### # 4 gaussian ### nls V5_nls4Qbands <- nls(blsub_V5 ~ gaussian(Qbands_domain, gau1.int, gau1.pos, gau1.fwhm) + gaussian(Qbands_domain, gau2.int, gau2.pos, gau2.fwhm) + gaussian(Qbands_domain, gau3.int, gau3.pos, gau3.fwhm) + gaussian(Qbands_domain, gau4.int, gau4.pos, gau4.fwhm), start = list(gau1.int = gau1_int, gau1.pos = gau1_pos, gau1.fwhm = gau1_fwhm, gau2.int = gau2_int, gau2.pos = gau2_pos, gau2.fwhm = gau2_fwhm, gau3.int = gau3_int, gau3.pos = gau3_pos, gau3.fwhm = gau3_fwhm, gau4.int = gau4_int, gau4.pos = gau4_pos, gau4.fwhm = gau4_fwhm), control=new.nls.control) # V5_nls4Qbands # gau1.int gau1.pos gau1.fwhm gau2.int gau2.pos gau2.fwhm gau3.int gau3.pos gau3.fwhm # 0.9658 506.9384 21.9889 0.5623 550.7035 20.7611 0.6417 580.3184 27.8472 # gau4.int gau4.pos gau4.fwhm # 0.4794 643.0718 16.4452 # residual sum-of-squares: 0.3518 ### predict the result of the model predict_V5_nls4Qbands <- predict(V5_nls4Qbands,list(x=Qbands_domain)) ### plot the result of the model along with the spectrum and baseline quartz() plot(x=Qbands_domain, y=blsub_V5, t="l", xlab="wavelength (nm)", ylab="Norm. intensity (a.u.)") lines(x=Qbands_domain, y=predict_V5_nls4Qbands, col="red") ## not bad at all! ## add the Gaussian lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V5_nls4Qbands)[1], coef(V5_nls4Qbands)[2], coef(V5_nls4Qbands)[3]), col="blue") lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V5_nls4Qbands)[4], coef(V5_nls4Qbands)[5], coef(V5_nls4Qbands)[6]), col="darkblue") lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V5_nls4Qbands)[7], coef(V5_nls4Qbands)[8], coef(V5_nls4Qbands)[9]), col="dodgerblue") lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V5_nls4Qbands)[10], coef(V5_nls4Qbands)[11], coef(V5_nls4Qbands)[12]), col="lightblue") ####### ### Gaussian decomposition of Glaucomys volans ZMB_Mam_60634 (ventral, non-glowing) ## V4 ####### # 4 gaussian ### nls V4_nls4Qbands <- nls(blsub_V4 ~ gaussian(Qbands_domain, gau1.int, gau1.pos, gau1.fwhm) + gaussian(Qbands_domain, gau2.int, gau2.pos, gau2.fwhm) + gaussian(Qbands_domain, gau3.int, gau3.pos, gau3.fwhm) + gaussian(Qbands_domain, gau4.int, gau4.pos, gau4.fwhm), start = list(gau1.int = gau1_int, gau1.pos = gau1_pos, gau1.fwhm = gau1_fwhm, gau2.int = gau2_int, gau2.pos = gau2_pos, gau2.fwhm = gau2_fwhm, gau3.int = gau3_int, gau3.pos = gau3_pos, gau3.fwhm = gau3_fwhm, gau4.int = gau4_int, gau4.pos = gau4_pos, gau4.fwhm = gau4_fwhm), control=new.nls.control) # V4_nls4Qbands # gau1.int gau1.pos gau1.fwhm gau2.int gau2.pos gau2.fwhm gau3.int gau3.pos gau3.fwhm # 0.9286 507.2414 23.9021 0.7319 548.1625 27.7205 0.9355 578.6463 28.2186 # gau4.int gau4.pos gau4.fwhm # 0.5656 639.0158 25.4653 # residual sum-of-squares: 0.3106 ### predict the result of the model predict_V4_nls4Qbands <- predict(V4_nls4Qbands,list(x=Qbands_domain)) ### plot the result of the model along with the spectrum and baseline quartz() plot(x=Qbands_domain, y=blsub_V4, t="l", xlab="wavelength (nm)", ylab="Norm. intensity (a.u.)") lines(x=Qbands_domain, y=predict_V4_nls4Qbands, col="red") ## not bad at all! ## add the Gaussian lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V4_nls4Qbands)[1], coef(V4_nls4Qbands)[2], coef(V4_nls4Qbands)[3]), col="blue") lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V4_nls4Qbands)[4], coef(V4_nls4Qbands)[5], coef(V4_nls4Qbands)[6]), col="darkblue") lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V4_nls4Qbands)[7], coef(V4_nls4Qbands)[8], coef(V4_nls4Qbands)[9]), col="dodgerblue") lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V4_nls4Qbands)[10], coef(V4_nls4Qbands)[11], coef(V4_nls4Qbands)[12]), col="lightblue") ####### ### Gaussian decomposition of Glaucomys sabrinus NRM_875246 (ventral, glowing) ## V3 ####### # 4 gaussian ### nls V3_nls4Qbands <- nls(blsub_V3 ~ gaussian(Qbands_domain, gau1.int, gau1.pos, gau1.fwhm) + gaussian(Qbands_domain, gau2.int, gau2.pos, gau2.fwhm) + gaussian(Qbands_domain, gau3.int, gau3.pos, gau3.fwhm) + gaussian(Qbands_domain, gau4.int, gau4.pos, gau4.fwhm), start = list(gau1.int = gau1_int, gau1.pos = gau1_pos, gau1.fwhm = gau1_fwhm, gau2.int = gau2_int, gau2.pos = gau2_pos, gau2.fwhm = gau2_fwhm, gau3.int = gau3_int, gau3.pos = gau3_pos, gau3.fwhm = gau3_fwhm, gau4.int = gau4_int, gau4.pos = gau4_pos, gau4.fwhm = gau4_fwhm), control=new.nls.control) # V3_nls4Qbands # gau1.int gau1.pos gau1.fwhm gau2.int gau2.pos gau2.fwhm gau3.int gau3.pos gau3.fwhm # 0.9529 505.7477 23.7981 0.4877 551.2204 22.5437 0.8117 578.5712 34.2168 # gau4.int gau4.pos gau4.fwhm # 0.5196 644.5041 19.8417 # residual sum-of-squares: 1.571 ### predict the result of the model predict_V3_nls4Qbands <- predict(V3_nls4Qbands,list(x=Qbands_domain)) ### plot the result of the model along with the spectrum and baseline quartz() plot(x=Qbands_domain, y=blsub_V3, t="l", xlab="wavelength (nm)", ylab="Norm. intensity (a.u.)") lines(x=Qbands_domain, y=predict_V3_nls4Qbands, col="red") ## not bad at all! ## add the Gaussian lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V3_nls4Qbands)[1], coef(V3_nls4Qbands)[2], coef(V3_nls4Qbands)[3]), col="blue") lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V3_nls4Qbands)[4], coef(V3_nls4Qbands)[5], coef(V3_nls4Qbands)[6]), col="darkblue") lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V3_nls4Qbands)[7], coef(V3_nls4Qbands)[8], coef(V3_nls4Qbands)[9]), col="dodgerblue") lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V3_nls4Qbands)[10], coef(V3_nls4Qbands)[11], coef(V3_nls4Qbands)[12]), col="lightblue") ####### ### Gaussian decomposition of Glaucomys volans ZMB_Mam_60634 (dorsal, non-glowing) ## V2 ####### # 4 gaussian ### nls V2_nls4Qbands <- nls(blsub_V2 ~ gaussian(Qbands_domain, gau1.int, gau1.pos, gau1.fwhm) + gaussian(Qbands_domain, gau2.int, gau2.pos, gau2.fwhm) + gaussian(Qbands_domain, gau3.int, gau3.pos, gau3.fwhm) + gaussian(Qbands_domain, gau4.int, gau4.pos, gau4.fwhm), start = list(gau1.int = gau1_int, gau1.pos = gau1_pos, gau1.fwhm = gau1_fwhm, gau2.int = gau2_int, gau2.pos = gau2_pos, gau2.fwhm = gau2_fwhm, gau3.int = gau3_int, gau3.pos = gau3_pos, gau3.fwhm = gau3_fwhm, gau4.int = gau4_int, gau4.pos = gau4_pos, gau4.fwhm = gau4_fwhm), control=new.nls.control) # V2_nls4Qbands # gau1.int gau1.pos gau1.fwhm gau2.int gau2.pos gau2.fwhm gau3.int gau3.pos gau3.fwhm # 0.5662 504.8255 22.5606 0.2074 537.0763 13.3881 0.4913 578.4355 50.2898 # gau4.int gau4.pos gau4.fwhm # 0.7542 638.8328 24.1951 # residual sum-of-squares: 2.241 ### predict the result of the model predict_V2_nls4Qbands <- predict(V2_nls4Qbands,list(x=Qbands_domain)) ### plot the result of the model along with the spectrum and baseline quartz() plot(x=Qbands_domain, y=blsub_V2, t="l", xlab="wavelength (nm)", ylab="Norm. intensity (a.u.)") lines(x=Qbands_domain, y=predict_V2_nls4Qbands, col="red") ## not bad at all! ## add the Gaussian lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V2_nls4Qbands)[1], coef(V2_nls4Qbands)[2], coef(V2_nls4Qbands)[3]), col="blue") lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V2_nls4Qbands)[4], coef(V2_nls4Qbands)[5], coef(V2_nls4Qbands)[6]), col="darkblue") lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V2_nls4Qbands)[7], coef(V2_nls4Qbands)[8], coef(V2_nls4Qbands)[9]), col="dodgerblue") lines(x=Qbands_domain, y=gaussian(Qbands_domain, coef(V2_nls4Qbands)[10], coef(V2_nls4Qbands)[11], coef(V2_nls4Qbands)[12]), col="lightblue")