-1

プロジェクトの一部としてプログラムを作成するために使用している Bruker NMR スペクトルがいくつかあります。私のプログラムは、実際のスペクトルで動作する必要があります。そこで、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.) 2 番目の列は、各ピークの強度を表します。

3.) 2 番目の列には、完全には整列していないが、最初の列に近い数値がいくつかあることに注意してください。例:-34613、-28925、-19395。これは重要なことだと思います。

完全な開示のために-私はRでプログラミングを行っています.

注: 私は Biostar でもこれを尋ねましたが、ここで質問に答えている人はあまりいないので、ここよりもここで答えを得る可能性が高いと思います。

編集:これは、私が見つけた1つの解決策がもっともらしいです:

友人から、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

再現可能なコードを使用した実際の例を次に示します。戦略やコーディングに関しては良いとは言えませんが、始めるきっかけにはなるでしょう。

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 に答える
2

パッケージPRocess には、スペクトルのピークを見つける機能があります。「ピーク発見R」などで検索すると、他にもたくさんあります。

于 2012-01-30T08:25:55.887 に答える