-1

我有一些 Bruker NMR 光谱，我用它们来创建一个程序作为项目的一部分。我的程序需要在实际频谱上工作。所以我将布鲁克核磁共振光谱的 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

这就是频谱：（来源：wisc.edu）替代文字

我的程序必须从这些数据中识别出峰值。所以我需要知道如何解释这些数字。以及它们如何准确地转换为光谱中的适当值。到目前为止，这是我学到的：

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'

如果你不明白，请告诉我。我会改进解释。

另外，我问了这个相关的问题。

4

2 回答 2

2

包PRocess 具有查找光谱峰值的功能。还有更多，如果您搜索“峰值查找 R”之类的内容

于 2012-01-30T08:25:55.887 回答
2

这是一个带有可重现代码的工作示例。我不认为它在策略或编码方面有任何好处，但它可以让你开始。

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")

测试

于 2012-01-31T09:14:28.167 回答