2

一連の地形プロファイル スキャンを組み合わせて、単一の連続プロファイルを作成したいと考えています。唯一の問題は、各スキャンが異なる高さから取得された場合とされていない場合があることです。そのため、異なるファイルにはカバーされる領域に関してかなりの量の重複がありますが、異なるデータには共通の基準点がない場合があります。絶対高さ。

以下は、4 つの異なるスキャンです。各スキャンには約 30 の測定値が含まれており、最後のいくつかの測定値は新しいデータを表し、残りは前のスキャンと重複しています。最初のスキャンには既知の絶対値のみが含まれているため、最初のスキャンは「ゴールド スタンダード」です。2 番目のスキャンはたまたま同じ高さから取得されるため、オーバーラップは (ほぼ) 完全に一致し、前のスキャンに 4 つの新しいポイントのみが追加されます。3 番目と 4 番目のスキャンは異なる高さから撮影されているため、オーバーラップが (相対的に) 同じ領域をカバーしていますが、前の 2 つのスキャンに単純につなぎ合わせることができません。

   Scan1<-c(5,6,7,8,15,16,18,20,25,23,20,17,15,10,10,9,8,9,11,10,13,16,17,19,20,25,28,30,29,30)
   Scan2<-c(15,16,18,20,25,23,20,16,15,10,10,9,8,9,11,10,13,16,17,19,20,25,28,30,29,30,32,35,38,37)
   Scan3<-c(28,25,23,18,18,17,16,17,19,18,21,23,25,27,26,33,36,37,37,38,40,43,46,45,43,42,40,38,32,30)
   Scan4<-c(27,30,29,36,39,39,40,41,43,46,49,48,46,45,43,41,35,33,30,29,28,30,31,32,35)

R を使用して、これら 4 つのスキャンをつなぎ合わせて連続した地形プロファイルを作成する方法はありますか? 絶対的な高さは、最初のスキャンに基づいて、連続する各スキャンが前のスキャンにステッチされる必要があります。IE- スキャン 2 がスキャン 1 にステッチされ、4 つのデータ ポイントが追加されます。次に、スキャン 3 からの新しいデータがスキャン 1 とスキャン 2 の組み合わせに追加され、スキャン 4 からの新しいデータがスキャン 1、2、および 3 の組み合わせに追加されます。 、 等々....

何らかのパターン認識を使用して、Scan3 が Scan1 と約 8 単位異なり、Scan4 が約 11 単位ずれていることを特定して、スキャン間の大きなオーバーラップを照合することにより、すべてのデータを正規化する方法があると想定しています。ただし、私のデータには「ノイズ」があり、オーバーラップのパターンが完全に適合しないことに注意してください。

最終結果には、4 つのスキャンすべてを含む完全な地形プロファイルが含まれている必要があり、実際の数値が異なる場合は何らかの調整が行われます。次のようなもの：

5,6,7,8,15,16,18,20,25,23,20,16.5,15,10,10,9,8,9,11,10,13,15.5,17,19,19,25,28,29.5,29,30,32,35,38,37,35,34,32,30,24,22,19,18,17,19,20,21,24
4

1 に答える 1

1

配列アラインメントを調べたいと思うかもしれません.DNAアラインメントは基本的にこの問題ですが、数字ではなく塩基が関係しています.

もう一度簡単に説明すると、スキャンをスライドさせたときに値間の差の最小標準偏差を見つけることに基づいて、「最良の」シフトを見つけるための簡単な関数を次に示します。この関数は、指定された 2 つのシーケンスを取り、それらを指定されたシフト (デフォルトでは -15 から 15) と比較します。

aligner <- function(bestsequence, sequence2, shift = (-15):15){
  minsd <- sd(bestsequence[1:min(length(sequence2), length(bestsequence))] - sequence2[1:min(length(sequence2), length(bestsequence))])
  bestshift <- 0
  avgdiff <- mean(bestsequence[1:min(length(sequence2), length(bestsequence))] - sequence2[1:min(length(sequence2), length(bestsequence))])
  for(i in shift){
    if(i < 0){
      worksequence2 <- sequence2[abs(i):length(sequence2)]
      if(sd(bestsequence[1:min(length(worksequence2), length(bestsequence))]
            -  worksequence2[1:min(length(worksequence2), length(bestsequence))]) < minsd){
        minsd <- sd(bestsequence[1:min(length(worksequence2), length(bestsequence))]-
                      worksequence2[1:min(length(worksequence2), length(bestsequence))])
        bestshift <- i
        avgdiff <- mean(bestsequence[1:min(length(worksequence2), length(bestsequence))]-
                          worksequence2[1:min(length(worksequence2), length(bestsequence))])
      }
    }
    if(i > 0){
      workbest <- bestsequence[i:length(bestsequence)]
      if(sd(workbest[1:min(length(sequence2), length(workbest))]
            -sequence2[1:min(length(sequence2), length(workbest))]) < minsd){
        minsd <- sd(workbest[1:min(length(sequence2), length(workbest))]-
                      sequence2[1:min(length(sequence2), length(workbest))])
        bestshift <- i
        avgdiff <- mean(workbest[1:min(length(sequence2), length(workbest))]-
                        sequence2[1:min(length(sequence2), length(workbest))])
      }
    }
  }
  return(list(bestshift = bestshift, avgdiff = avgdiff, minsd = minsd))
}

だから、あなたのデータのために：

aligner(Scan1, Scan2)

$bestshift
[1] 5

$avgdiff
[1] 0.03846154

$minsd
[1] 0.1961161

したがって、Scan2 の 5 番目の要素は Scan1 の 1 番目の要素と同じです。ここから、サブセット化し、avgdiff で修正し、新しいデータ ポイントを追加して、再実行するのは簡単です。

編集：最終的なシーケンスを取得する方法は次のとおりです。まず、目的のシーケンスを出力するラッパーが必要です。基本的には前のコマンドを実行し、シフトが正か負かをチェックしてから、正しいシーケンスを出力します。

wrappedaligner <- function(bestseq, seq2){
  z <- aligner(bestseq, seq2)
  if(z$bestshift==0){
    if(length(bestseq) >= length(seq2)){
    return(bestseq)
    } else {return(c(bestseq, (seq2[(length(bestseq)+1):length(seq2)])-z$avgdiff))}
  }
  else if(z$bestshift > 0){
    if(length(bestseq)-z$bestshift >= length(seq2)){
      return(bestseq)
  } else {return(c(bestseq, seq2[(length(bestseq) - z$bestshift + 2):length(seq2)] - z$avgdiff))}
  }
  else if(z$bestshift <0){
    if((length(bestseq) - abs(z$bestshift))>= length(seq2)){
      return(bestseq)
    } else {return(c(seq2[1:abs(z$bestshift) - 1] - z$avgdiff, bestseq))}
  }
}

次に、データに対して再帰的に実行する必要があります-幸いなことに、次を使用できますReduce。

Reduce(wrappedaligner, list(Scan1, Scan2, Scan3, Scan4))

 [1]  5.00000  6.00000  7.00000  8.00000 15.00000 16.00000 18.00000 20.00000
 [9] 25.00000 23.00000 20.00000 17.00000 15.00000 10.00000 10.00000  9.00000
[17]  8.00000  9.00000 11.00000 10.00000 13.00000 16.00000 17.00000 19.00000
[25] 20.00000 25.00000 28.00000 30.00000 29.00000 30.00000 31.96154 34.96154
[33] 37.96154 36.96154 50.83974 49.83974 47.83974 45.83974 39.83974 37.83974
于 2015-09-14T20:23:18.347 に答える