Преобразование файла NMR ascii в список пиков

Question

У меня есть несколько спектров ЯМР Брукера, которые я использую для создания программы в рамках проекта.Моя программа должна работать на реальном спектре.Поэтому я преобразовал 1r файлы спектров ЯМР Брукера в ASCII.Для карнитина это то, как выглядит файл ascii (это не полный список. Полный список состоит из тысяч строк. Это всего лишь снимок):

-0.807434   -23644  
-0.807067   -22980  
-0.806701   -22967  
-0.806334   -24513  
-0.805967   -27609  
-0.805601   -31145  
-0.805234   -33951  
-0.804867   -35553  
-0.804501   -35880  
-0.804134   -35240  
-0.803767   -34626  
-0.8034  -34613 
-0.803034   -34312  
-0.802667   -32411  
-0.8023  -28925 
-0.801934   -25177  
-0.801567   -22132  
-0.8012  -19395 

и это то, что представляет собой спектр: альтернативный текст http://www.bmrb.wisc.edu/metabolomics/standards/L_carnitine/nmr/bmse000211/spectra_png/1H.png

Моя программа должна идентифицировать пики по этим данным.Поэтому мне нужно знать, как интерпретировать эти цифры.И как именно они преобразуются в их соответствующие значения в спектре. До сих пор это то, что я узнал:

1.) Первый столбец представляет положение спектральной точки (ppm)

2.) Во втором столбце представлена ​​интенсивность каждого пика.

3.) Обратите внимание, что во втором столбце есть некоторые числа, которые не идеально выровнены, но находятся ближе к первому столбцу.Например: -34613, -28925, -19395.Я думаю, что это важно.

Ради полного раскрытия - я занимаюсь программированием в R.

ПРИМЕЧАНИЕ: я также спрашивал об этом в Biostar, но я думаю, что у меня больше шансовполучить ответ здесь, чем там, потому что не многие люди, кажется, отвечают там на вопросы.

РЕДАКТИРОВАТЬ: Это одно решение, которое я нашел, правдоподобно:

Друг дал мне идеюиспользуйте скрипт awk, чтобы проверить, где именно интенсивность файла меняется с положительной на отрицательную, чтобы найти локальные максимумы.Вот рабочий скрипт:

awk 'BEGIN{dydx = 0;}
{ 
  if(NR > 1)
     { dydx = ($2 - y0)/($1 - x0); } 
  if(NR > 2 && last * dydx < 0)
     { printf( "%.4f  %.4f\n", (x0 + $1)/2, log((dydx<0)?-dydx:dydx)); } ;
  last=dydx; x0=$1; y0=$2
}' /home/chaitanya/Work/nmr_spectra/caffeine/pdata/1/spectrumtext.txt  | awk '$2 > 17'

Скажите, если вы этого не понимаете.Я улучшу объяснение.

Кроме того, есть этот связанный вопрос, который я задал.

Answer 1

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

find_peaks <- function (x, y, n.fine = length(x), interval = range(x), ...) {
  maxdif <- max(diff(x)) # longest distance between successive points

  ## selected interval for the search
  range.ind <- seq(which.min(abs(x - interval[1])),
                   which.min(abs(x - interval[2])))
  x <- x[range.ind]
  y <- y[range.ind]

  ## smooth the data
  spl <- smooth.spline(x, y, ...)
  ## finer x positions
  x.fine <- seq(range(x)[1], range(x)[2], length = n.fine)
  ## predicted y positions
  y.spl <- predict(spl, x.fine, der = 0)$y
  ## testing numerically the second derivative
  test <- diff(diff((y.spl), 1) > 0, 1)
  maxima <- which(test == -1) + 1

  ## according to this criterion, we found rough positions
  guess <- data.frame(x=x.fine[maxima], y=y.spl[maxima])

  ## cost function to maximize 
  obj <- function(x) predict(spl, x)$y

  ## optimize the peak position around each guess
  fit <- data.frame(do.call(rbind,
          lapply(guess$x, function(g) {
            fit <- optimize(obj, interval = g + c(-1,1) * maxdif, maximum=TRUE)
            data.frame(x=fit$maximum,y=fit$objective)
          })))

  ## return both guesses and fits
  invisible(list(guess=guess, fit=fit))
}

set.seed(123)
x <- seq(1, 15, length=100)
y <- jitter(cos(x), a=0.2)

plot(x,y)
res <- find_peaks(x,y)
points(res$guess,col="blue")
points(res$fit,col="red")

Answer 2

Пакет PRocess имеет функцию для нахождения пиков в спектрах. Есть много других, если вы ищете что-то вроде "пика находки R"