根据已知的重叠将地形剖面拼接在一起

Question

我有一系列地形剖面扫描，我想将它们结合起来创建一个连续的剖面。唯一的问题是每次扫描可能会或可能不会从不同的高度进行，因此虽然不同的文件在覆盖区域方面有相当多的重叠，但不同的数据可能没有共同的参考点绝对高度。

以下是 4 种不同的扫描。每次扫描包含大约 30 个测量值，最后几个测量值代表新数据，其余的与前一次扫描重叠。第一次扫描包含唯一已知的绝对值，因此第一次扫描是 "gold standard"。第二次扫描恰好是从相同的高度进行的，因此重叠（几乎）完美匹配并且只为之前的扫描添加了 4 个新点。第三次和第四次扫描是从不同的高度拍摄的，所以虽然重叠覆盖了相同的区域（相对），但我不能简单地将它拼接到前两次扫描上。

   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- Scan2 被缝合到 Scan 1 上添加 4 个数据点，然后 Scan 3 的新数据被添加到 Scan1 和 Scan2 的组合中，然后 Scan4 的新数据被添加到 Scans 1,2 和 3 的组合中,等等....

我假设有一种方法可以通过匹配扫描之间的大重叠来标准化所有数据，使用某种模式识别来确定 Scan3 与 Scan1 大约有 8 个单位不同，而 Scan4 是大约 11 个单位。但请注意，我的数据中有一些 "noise"，重叠的模式并不完全合适。

最终结果应该包含一个完整的地形剖面，包含所有 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

Answer 1

您可能需要研究序列比对 - DNA 比对基本上就是这个问题，但使用的是碱基而不是数字。

快速浏览一下，这是一个快速编写的函数，用于查找 "best" 偏移，基于滑动扫描时查找值之间差异的最低标准差。该函数采用给定的两个序列，并将它们与给定的位移（默认为 -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

根据已知的重叠将地形剖面拼接在一起

Stiching together topographic profiles based on known overlap

r

pattern-matching

scanning

topography