500 Raman spectra were acquired for 50 positions of the gratings.
## list the data
load("fluo532.rda")
names(fluo532)[1:3]
## [1] "fluo532/shift_0_l532nm_0-15secx500spectra_fluorescein_boratebuffer_pH9-1_para_pol_triple_P68mW_coupled_sub_low_gain_1800lmm_binningrange_90_160.txt"
## [2] "fluo532/shift_1_l532nm_0-15secx500spectra_fluorescein_boratebuffer_pH9-1_para_pol_triple_P68mW_coupled_sub_low_gain_1800lmm_binningrange_90_160.txt"
## [3] "fluo532/shift_10_l532nm_0-15secx500spectra_fluorescein_boratebuffer_pH9-1_para_pol_triple_P68mW_coupled_sub_low_gain_1800lmm_binningrange_90_160.txt"
## rename and reorder according to shift number
raw <- fluo532
orig_names <- names(raw)
fnames <- gsub("(.*)shift_([0-9]+)_(.*)","\\2", orig_names)
names(raw) <- fnames
raw <- raw[order(as.numeric(names(raw)))]
fnames <- names(raw)
ldply(raw, each(ncol, nrow, min, max, mean, sd = function(x) sd(colMeans(x))))
## .id ncol nrow min max mean sd
## 1 0 503 1024 1 15647 11438 738.2
## 2 1 502 1024 1 18080 11413 737.2
## 3 2 502 1024 1 15533 11377 734.8
## 4 3 503 1024 1 17008 11345 731.7
## 5 4 503 1024 1 42081 11305 728.8
## 6 5 503 1024 1 15488 11263 725.8
## 7 6 503 1024 1 15373 11215 722.2
## 8 7 503 1024 1 15647 11163 718.7
## 9 8 502 1024 1 15385 11109 715.8
## 10 9 502 1024 1 17453 11047 711.6
## 11 10 503 1024 1 15076 10998 707.3
## 12 11 503 1024 1 15016 10954 704.3
## 13 12 503 1024 1 14954 10914 701.6
## 14 13 503 1024 1 14929 10862 697.9
## 15 14 503 1024 1 14889 10817 694.6
## 16 15 503 1024 1 14781 10783 692.1
## 17 16 502 1024 1 14730 10745 690.3
## 18 17 502 1024 1 14716 10711 687.8
## 19 18 503 1024 1 14665 10673 684.4
## 20 19 503 1024 1 20759 10637 682.0
## 21 20 503 1024 1 14576 10611 679.8
## 22 21 503 1024 1 14504 10569 677.1
## 23 22 503 1024 1 14498 10543 675.3
## 24 23 503 1024 1 14494 10515 672.9
## 25 24 502 1024 1 14369 10481 671.6
## 26 25 503 1024 1 14767 10449 668.4
## 27 26 503 1024 1 14313 10424 666.8
## 28 27 503 1024 1 14278 10401 665.0
## 29 28 503 1024 1 15739 10370 663.0
## 30 29 503 1024 1 15480 10338 661.0
## check intensities
check <- llply(raw, average_intensity, exclude=1:2) # remove first spectra
m <- melt(check, id="spectrum")
m$grating <- factor(m$L1, levels=fnames)
p0 <-
ggplot( data=m, aes(spectrum, value, colour=grating)) +
facet_wrap(~grating,ncol=5) +
geom_path() + guides(colour="none")+
labs(x = expression("spectrum #"),
y = expression("average counts"))
p0
## take the average of each spectra matrix
avg <- llply(raw, average_spectra, exclude=1:2)
positions <- sapply(avg, function(d) mean(range(d$wavenumber)))
## remove spikes
avgc <- llply(avg, spike_cleanup, 1e-1)
## plot cleaned averages
m <- melt(avgc, id=c("wavenumber", "pixel"), meas="value")
m$grating <- as.numeric(factor(m$L1, levels=fnames))
p1 <-
ggplot( data=m, aes(wavenumber, value, colour=grating, group=grating)) +
geom_path() +
scale_colour_gradientn(breaks=c(0, seq(10, 50, by=10)), limit=c(0, 50), expand=c(0, 0),
colours = hue_pal(h = c(60, 240), c = 100, l = 65, h.start = 0,
direction = 1)(50)) +
guides(colour=guide_colourbar(direction = "horizontal", barheight = unit(2, "mm"),
title.position = "top"))+
labs(x = expression("Raman shift ["*cm^-1*"]"),
y = expression("Intensity [counts]"),
colour = expression(Shift/cm^-1)) +
theme(legend.position="top")
p1
p2 <- p1 + aes(x=pixel) +
scale_x_continuous("pixel #", breaks=c(1,seq(0,1024,by=256)),
minor_breaks=seq(0,1024,by=128),
expand=c(0,0), lim=c(1, 1024))
p2
m2 <- subset(m, pixel < 900 & pixel > 400)
structure <- ddply(m, "grating",
transform, value=background_suppression(value, wavenumber, 6))
p3 <- p2 %+% structure +
guides(colour="none") +
scale_y_continuous("")
p3
all <- melt(avgc, meas="value")
## all <- subset(all, pixel > 200 & pixel < 700)
naive <- kiefer_average(all,"L1")
naive <- transform(subset(naive, wavenumber < 1600 & wavenumber > 1200),
subtracted = background_suppression(value, wavenumber, 4))
p4 <-
ggplot(naive, aes(wavenumber, subtracted)) +
geom_path() +
labs(x = expression("Raman shift ["*cm^-1*"]"),
y = expression("Intensity [counts]"),
colour = expression(Shift/cm^-1))
p4
## madeup peak to test the analysis
## all <- ddply(all, "L1", mutate, value=value+lorentz(wavenumber, 650, peak=10, width=6))
all$shift <- as.numeric(all$L1)
all_structure <- ddply(all, .(L1), mutate,
subtracted = background_suppression(value, pixel, 4))
all_average <- ddply(all, .(pixel), summarize, value=mean(value))
all_average <- transform(all_average,
subtracted = background_suppression(value, pixel, 4))
sub <- subset(all_structure, L1 %in% c("0", "49"))
p5 <-
ggplot(all_structure)+
geom_path(aes(pixel, subtracted, group=L1), colour="grey") +
geom_path(aes(pixel, subtracted), data=all_average, colour="#E41A1C") +
scale_x_continuous("pixel #", breaks=c(1,seq(0,1024,by=256)),
minor_breaks=seq(0,1024,by=128),
expand=c(0,0), lim=c(1, 1024))
p5
differences_all <- ddply(all, .(L1), weighted_difference, reference=all_average,
n=4, exclude=c(580, 610))
differences_all <- ddply(differences_all, .(L1), mutate,
difference=background_suppression(difference, wavenumber, 6))
differences_all$L1 <- sort_factor(differences_all$L1)
p6 <-
ggplot(differences_all)+
## facet_grid(L1~., scales="free", space="free") +
geom_path(aes(pixel, difference + 8*L1, colour=L1, group=L1)) +
scale_y_continuous(expand=c(0, 0), breaks=NULL)+
scale_x_continuous("pixel #", breaks=c(1,seq(0,1024,by=256)),
minor_breaks=seq(0,1024,by=128),
expand=c(0,0), lim=c(1, 1024)) +
scale_colour_gradientn(breaks=c(0, seq(10, 50, by=10)), limit=c(0, 50), expand=c(0, 0),
colours = hue_pal(h = c(60, 240) , c = 100, l = 60, h.start = 0,
direction = 1)(50)) +
guides(colour="none")
p6
combined <- sum_spectra(differences_all)
final <- subset(combined, type=="average", select=c("wavenumber", "value"))
p7 <-
ggplot(final)+
geom_path(aes(wavenumber, value)) +
labs(x = expression("Raman shift ["*cm^-1*"]"),
y = expression("Intensity [counts]"),
colour = expression(Shift/cm^-1))
p7