Использование R для анализа FTIR-спектров

Question

Это мой первый пост, извините, если с ним возникнут проблемы. Я знаю, что это очень долго, но я не знал, как «разложить» свой вопрос.

Я пытаюсь написать код для анализа FTIR спектров. По сути, то, что я пытаюсь сделать, - это идентифицировать пики в спектрах и измерять их высоту.

Я оглядывался и не смог найти пример, поэтому сам разработал его с нуля.

Шаги:

1. Сглаживание спектра.
2. Определение максимумов (пиков) и минимумов (базовая линия для измерения высоты пиков).
3. Создание базовых линий.
4. Измерение пиков с использованием базовой линии в качестве эталона.

Я пробовал сглаживание с помощью средств прокатки, лёсса, сплайна. Сглаживания, гаммы ... и лучшего результаты получены из скользящих средних.

Как вы можете видеть на рисунке базовых линий и на последнем рисунке, основная проблема заключается в некоторых базовых линиях. Любые предложения о том, как это улучшить, будут приветствоваться.

И в целом, я был бы благодарен за любое предложение, как улучшить этот код.

Git repo: https://github.com/paracon/peakfinder

Итак, код выглядит так:

#LOADING LIBRARIES

library(ggplot2)


#LOADING DATA

ftir <- read.table(file="https://raw.githubusercontent.com/paracon/peakfinder/master/DATA/spectra.csv", 
                   header=TRUE, 
                   dec=".",
                   sep=",", 
                   stringsAsFactors = FALSE)

ftir$replicate <- as.factor(ftir$replicate)
ftir$treatment <- as.factor(ftir$treatment)
ftir$wavenumber <- as.numeric(ftir$wavenumber)


ftir <- ftir[ftir$treatment=="trt1" & ftir$replicate=="rep2",]


#ROLLING MEANS SMOOTHING

rolling_means_sensitivity <- 15

ftir$absorbance_smoothed <- zoo::rollmean(ftir$absorbance, k = rolling_means_sensitivity, fill = NA)



#PLOTTING SMOOTHED VS ORIGINAL DATA

ggplot(ftir, aes(x=wavenumber)) +
  geom_line(aes(y=absorbance, colour="blue"), size=.7) +
  geom_line(aes(y=absorbance_smoothed, colour="red"), size=.7) +
  theme_bw()

К сожалению, я не могу загружать изображения, поэтому я добавляю на них ссылки:

Исходный и сглаженный спектр: https://raw.githubusercontent.com/paracon/peakfinder/master/images/original_vs_smoothed.png

#DETECTING CHANGES OF SLOPE

lag <- c(0, head(ftir$absorbance_smoothed, -1))
ftir <- cbind(ftir, lag)

ftir$absorbance_smoothed_diff <- ftir$absorbance_smoothed - ftir$lag
try(ftir[is.na(ftir$absorbance_smoothed_diff),]$absorbance_smoothed_diff <- 0)
ftir$change <- NA
ftir[ftir$absorbance_smoothed_diff<0,]$change  <- "smaller"
ftir[ftir$absorbance_smoothed_diff>0,]$change  <- "bigger"

ftir_changes <- ftir[!is.na(ftir$change),]
ftir_changes <- ftir_changes[which(diff(sign(ftir_changes$absorbance_smoothed_diff))!=0),]


#PLOTTING PEAKS AND VALLEYS

ggplot(ftir, aes(x=wavenumber)) +
  geom_line(aes(y=absorbance, colour="blue"), size=.7) +
  #geom_line(aes(y=absorbance_smoothed, colour="red"), size=.7) +
  geom_point(data=ftir_changes, aes(y=absorbance_smoothed)) +
  #geom_point(size=.3) +
  theme_bw()

Пики и спады: https://raw.githubusercontent.com/paracon/peakfinder/master/images/peaks_and_valleys.png

#IDENTIFYING VALLEYS

ftir_valleys <- ftir_changes[which(diff(sign(ftir_changes$absorbance_smoothed_diff))>0),]
ftir_valleys <- rbind(ftir_valleys, ftir[ftir$wavenumber==min(ftir$wavenumber, na.rm=TRUE),])
ftir_valleys <- rbind(ftir_valleys, ftir[ftir$wavenumber==max(ftir$wavenumber, na.rm=TRUE),])


#GENERATING BASELINE FROM VALLEYS

ftir_valleys <- dplyr::left_join(ftir[,c("treatment","wavenumber")], ftir_valleys[,c("treatment","wavenumber", "absorbance")], by=c("treatment","wavenumber"))

ftir_valleys$absorbance_interp <- zoo::na.approx(ftir_valleys$absorbance, na.rm="FALSE")

Исходные данные: https://raw.githubusercontent.com/paracon/peakfinder/master/images/baselines.png

#IDENTIFYING PEAKS

ftir_peaks <- ftir_changes[which(diff(sign(ftir_changes$absorbance_smoothed_diff))<0),]


ftir_peaks$absorbance_smoothed_orig <- ftir_peaks$absorbance_smoothed

for (i in (1:nrow(ftir_peaks))){
  #i=10
  rownumber_i <- as.integer(rownames(ftir[ftir$absorbance==ftir_peaks[i,]$absorbance & ftir$wavenumber==ftir_peaks[i,]$wavenumber,]))
  #ftir_peaks[i,]$loess2 <- max(ftir[seq(min(rownumber_i, na.rm=TRUE)-(rolling_means_sensitivity/2), max(rownumber_i,na.rm=TRUE)+(rolling_means_sensitivity/2),1),]$absorbance)
  try(ftir_peaks[i,]$absorbance_smoothed <- max(
    ftir[rownames(ftir)>(min(rownumber_i,na.rm=TRUE)-(rolling_means_sensitivity/2)) & rownames(ftir)<(max(rownumber_i,na.rm=TRUE)+(rolling_means_sensitivity/2)),]$absorbance, na.rm=TRUE
  ), silent=TRUE)
  if (ftir_peaks[i,]$absorbance_smoothed==-Inf){
    ftir_peaks[i,]$absorbance_smoothed <- ftir_peaks[i,]$absorbance_smoothed_orig
  }
  if (ftir_peaks[i,]$absorbance_smoothed==Inf){
    ftir_peaks[i,]$absorbance_smoothed <- ftir_peaks[i,]$absorbance_smoothed_orig
  }
  if (is.na(ftir_peaks[i,]$absorbance_smoothed)){
    ftir_peaks[i,]$absorbance_smoothed <- ftir_peaks[i,]$absorbance_smoothed_orig
  }
}


#CALCULATING THE PEAK HEIGHT

ftir_peaks <- dplyr::left_join(ftir_peaks[,c("treatment","wavenumber", "absorbance", "absorbance_smoothed")], ftir_valleys[,c("treatment","wavenumber", "absorbance_interp")], by=c("treatment","wavenumber"))
ftir_peaks$height <- NA
ftir_peaks$height <- ftir_peaks$absorbance_smoothed - ftir_peaks$absorbance_interp


#PLOTTING PEAK HEIGHTS

ggplot(ftir, aes(x=wavenumber)) +
  scale_x_reverse() + 
  geom_line(aes(y=absorbance, colour="blue"), size=.7) +
  #geom_line(aes(y=absorbance_smoothed, colour="red"), size=.7) +
  geom_line(data=ftir_valleys, aes(y=absorbance_interp)) +
  theme_bw() +
  annotate("text", x=ftir_peaks$wavenumber, y=ftir_peaks$absorbance, label=round(ftir_peaks$height,3),check_overlap = TRUE)

Высота пиков: https://raw.githubusercontent.com/paracon/peakfinder/master/images/peak_heights.png