Соединение топографических профилей на основе известного перекрытия
У меня есть серия сканирований топографических профилей, которые я хотел бы объединить, чтобы создать единый непрерывный профиль. Единственная проблема заключается в том, что каждое сканирование может или не может быть взято с другой высоты, поэтому, хотя разные файлы имеют достаточное перекрытие с точки зрения охватываемой области, разные данные могут не иметь общую контрольную точку с точки зрения абсолютная высота
Ниже представлены 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 накладывается на Scan1, добавляя 4 точки данных, затем новые данные из Scan 3 добавляются в комбинацию Scan 1 и Scan2, затем новые данные из Scan4 добавляются в комбинацию Scan1,2 и 3., и так далее....
Я предполагаю, что есть способ нормализовать все данные путем сопоставления большого перекрытия между сканированиями, используя какое-то распознавание образов, чтобы выяснить, что Scan3 примерно на 8 единиц отличается от Scan1, а 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
1 ответ
Возможно, вы захотите взглянуть на выравнивание последовательностей - выравнивание ДНК в основном это проблема, но с основами, а не числами.
Еще раз, вот быстрая функция для поиска "лучшего" сдвига, основанная на поиске наименьшего стандартного отклонения различий между значениями при скольжении сканов. Функция берет две заданные последовательности и сравнивает их с заданными сдвигами (по умолчанию от -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
Итак, ваш 5-й элемент Scan2s равен 1-му элементу 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