2

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

以下是 4 种不同的扫描。每次扫描包含大约 30 次测量，最后几次测量代表新数据，其余的与前一次扫描重叠。第一次扫描包含唯一已知的绝对值，因此第一次扫描是“黄金标准”。第二次扫描恰好是从相同的高度进行的，因此重叠（几乎）完美匹配，并且仅在前一次扫描中增加了 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 个单位。但请注意，我的数据中有一些“噪音”，重叠的模式并不完美。

最终结果应包含包含所有 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 比对基本上就是这个问题，但使用的是碱基而不是数字。

快速回顾一下，这里有一个快速编写的函数来找到“最佳”偏移，基于在滑动扫描时找到值之间差异的最低标准偏差。该函数采用给定的两个序列，并将它们与给定的移位（默认为 -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

因此，您的 Scan2s 第 5 个元素等于 Scan1 的第一个元素。从这里它应该很容易进行子集化，通过 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 回答