将 NMR ascii 文件转换为峰列表

Posted 2023-02-25

【中文标题】将 NMR ascii 文件转换为峰列表

【英文标题】：Converting NMR ascii file to peak list

【发布时间】：2012-01-30 07:18:29

【问题描述】：

我有一些 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 'BEGINdydx = 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'

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

另外，我问了this 相关的问题。

【问题讨论】：

如果您要在 R 中执行其余的分析，我不明白为什么需要 awk 脚​​本来完成这个简单的任务。为什么不学习如何在R?

可以在 R 中完成是的。我最终可能会用 R 编码。如果用 R 编码，我的代码可能也会更快。这只是在与朋友讨论时提出的问题，我想知道是否可以这样做。跨度>

@baptiste 题外话 - 为什么我的代码中的 Work 这个词是蓝色的？

不知道，我的猜测是第一个字母是大写的，语法荧光笔认为它是某种类或其他东西。或者也许这是一个微妙的潜意识暗示，我们应该回去工作。

wrt 3)：无需担心空格（它们实际上是空格还是单个制表符？）。在任何一种情况下，read.table 都可以轻松导入数据，如果您想要处理光谱的其他可能性，您可以查看包 hyperSpec（它不进行峰值查找，但各种光谱绘图等.).

【参考方案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")

【参考方案2】：

光谱中的包裹PRocess has a function to find peaks。还有更多，如果您搜索“peak Finding R”之类的内容

【讨论】：

这似乎是一个很好的资源。但我仍然想知道 ascii 文件中的信息是如何准确映射到频谱上的。

最近一期 R 杂志专门讨论磁共振成像；我认为您可以通过阅读一些文章来了解一些背景信息。

当我搜索 Peak Finding R 时，我发现的只是邮件存档，其中指定了算法和 Matlab 包。有什么建议吗？

RSiteSearch 和 sos::findFn 都为我提供了带有“峰值查找”或“峰值检测”等搜索词的有希望的结果