Как искать долины на графике?

10

Я изучаю некоторые данные о геномном покрытии, которые в основном представляют собой длинный список (несколько миллионов значений) целых чисел, каждый из которых говорит о том, насколько хорошо (или «глубоко») охвачена эта позиция в геноме.

Я хотел бы найти «долины» в этих данных, то есть регионы, которые значительно «ниже», чем их окружение.

Обратите внимание, что размер долин, которые я ищу, может варьироваться от 50 оснований до нескольких тысяч.

Какую парадигму вы бы порекомендовали использовать для поиска этих долин?

ОБНОВИТЬ

Некоторые графические примеры для данных: альтернативный текст альтернативный текст

ОБНОВЛЕНИЕ 2

Определение долины - это, конечно, один из вопросов, с которыми я борюсь. Это очевидные для меня: альтернативный текст альтернативный текст

но есть несколько более сложных ситуаций. В целом, есть три критерия, которые я рассматриваю: 1. (Среднее? Максимальное?) Покрытие в окне по отношению к глобальному среднему. 2. Покрытие (...) в окне относительно его непосредственного окружения. 3. Насколько велико окно: если я вижу очень низкий охват для короткого промежутка, это интересно, если я вижу очень низкий охват для длинного промежутка, это также интересно, если я вижу слегка низкое покрытие для короткого промежутка, это не очень интересно , но если я вижу слегка низкий охват для длинного промежутка - это .. Так что это сочетание длины sapn и его покрытия. Чем дольше, тем выше я оставляю покрытие и все равно считаю его долиной.

Спасибо,

Дейв

Дэвид Б
источник
Не могли бы вы предоставить небольшой образец данных?
Шейн
@Shane см. Обновление
Дэвид Б.
@ Дэвид Спасибо. Поскольку оба ответа подразумевают, анализ временных рядов может быть применен здесь, так как вы заказали наблюдения.
Шейн
Это довольно сложно ответить, не зная точно, что вы ищете. Можете ли вы обвести точки на графиках, которые вы хотите захватить? Что вы считаете "долиной"? как низко это должно идти и что вы хотите вернуть? Трудно сформулировать решение, не зная вопроса, то есть порогов и тому подобное.
Фалмарри,
@ Шейн ♦ Спасибо. Поскольку у меня нет опыта и в анализе временных рядов, не могли бы вы оставить несколько указаний, с чего мне начать?
Дэвид Б

Ответы:

5

Вы можете использовать своего рода подход Монте-Карло, например, используя скользящее среднее ваших данных.

Возьмите скользящее среднее данных, используя окно разумного размера (я думаю, что вам решать, насколько широко).

Пропуски в ваших данных будут (конечно) характеризоваться более низким средним, поэтому теперь вам нужно найти какой-то «порог», чтобы определить «низкий».

Для этого вы случайным образом меняете значения ваших данных (например, используя sample() ) и пересчитываете скользящее среднее для ваших замененных данных.

Повторите этот последний отрывок достаточно много раз (> 5000) и сохраните все средние значения этих испытаний. По сути, у вас будет матрица из 5000 строк, по одной на пробу, каждая из которых содержит скользящую среднюю для этого испытания.

В этот момент для каждого столбца вы выбираете квантиль 5% (или 1% или что вы хотите), то есть значение, под которым лежит только 5% от среднего значения рандомизированных данных.

Теперь у вас есть «доверительный предел» (я не уверен, что это правильный статистический термин) для сравнения ваших исходных данных. Если вы обнаружите, что часть ваших данных ниже этого предела, вы можете назвать это сквозным.

Конечно, следует помнить, что ни этот, ни какой-либо другой математический метод никогда не могли бы дать вам никаких указаний на биологическую значимость, хотя я уверен, что вы хорошо это знаете.

РЕДАКТИРОВАТЬ - пример

require(ares) # for the ma (moving average) function

# Some data with peaks and throughs 
values <- cos(0.12 * 1:100) + 0.3 * rnorm(100) 
plot(values, t="l")

# Calculate the moving average with a window of 10 points 
mov.avg <- ma(values, 1, 10, FALSE)

numSwaps <- 1000    
mov.avg.swp <- matrix(0, nrow=numSwaps, ncol=length(mov.avg))

# The swapping may take a while, so we display a progress bar 
prog <- txtProgressBar(0, numSwaps, style=3)

for (i in 1:numSwaps)
{
# Swap the data
val.swp <- sample(values)
# Calculate the moving average
mov.avg.swp[i,] <- ma(val.swp, 1, 10, FALSE)
setTxtProgressBar(prog, i)
}

# Now find the 1% and 5% quantiles for each column
limits.1 <- apply(mov.avg.swp, 2, quantile, 0.01, na.rm=T)
limits.5 <- apply(mov.avg.swp, 2, quantile, 0.05, na.rm=T)

# Plot the limits
points(limits.5, t="l", col="orange", lwd=2)
points(limits.1, t="l", col="red", lwd=2)

Это только позволит вам графически найти регионы, но вы можете легко найти их, используя что-то на линии which(values>limits.5).

Nico
источник
Очевидно, что вы можете применить тот же подход, используя нечто иное, чем скользящее среднее, это просто для того, чтобы дать представление.
Нико
+1 Большое спасибо, Нико. Позвольте мне увидеть, правильно ли я вас понял: в конце концов, это в основном похоже на установление некоторого глобального порога и определение любой точки со значением <порог как части долины. Выборка и т. Д. Просто используется, чтобы получить значимую меру (квантиль) для установки порога. Почему мы не можем использовать один порог для всех точек, я имею в виду, что если мы выполним достаточное количество симуляций, мы получим прямые (чётные и желтые) линии. Кроме того, поправьте меня, если я ошибаюсь, но это не учитывает окружающую среду, а проверяет абсолютную ценность каждой точки.
Дэвид Б
@ Дэвид Б: конечно, вы можете использовать глобальный порог, и это, вероятно, сэкономит вам некоторое время на вычисления. Я предполагаю, что выбор чего-то вроде 1/3 глобального среднего значения может быть началом. Этот процесс обмена, вероятно, более полезен, если вы используете какую-то другую статистику, а не скользящую среднюю, это было в основном дать представление. В любом случае скользящая средняя будет учитывать окружение, в примере будет учитываться окно из 10 пунктов.
Нико
4

Я совершенно не знаю этих данных, но, предполагая, что данные упорядочены (не по времени, а по положению?), Имеет смысл использовать методы временных рядов. Существует множество методов идентификации временных кластеров в данных. Обычно они используются для поиска высоких значений, но могут использоваться для низких значений, сгруппированных вместе. Я имею в виду статистику сканирования, статистику суммарной суммы (и другие), используемые для обнаружения вспышек заболеваний в данных подсчета. Примеры этих методов находятся в пакете наблюдения и пакете DCluster.


источник
@cxr Спасибо за ваш ответ. Я взглянул на surveillanceи DCluster , но не могли бы вы быть более конкретным? Они оба относительно большие пакеты, и их цель кажется довольно конкретной. Я не уверен, с чего начать.
Дэвид Б
2

Есть много вариантов для этого, но один хороший: вы можете использовать msExtremaфункцию в msProcessпакете .

Редактировать:

В анализе финансовых показателей этот вид анализа часто выполняется с использованием концепции «просадки». PerformanceAnalyticsПакет имеет некоторые полезные функции , чтобы найти эти долины . Вы можете использовать тот же алгоритм здесь, если вы рассматриваете свои наблюдения как временные ряды.

Вот несколько примеров того, как вы могли бы применить это к своим данным (где «даты» не имеют значения, но используются только для упорядочивания), но первыми элементами в zooобъекте будут ваши данные:

library(PerformanceAnalytics)
x <- zoo(cumsum(rnorm(50)), as.Date(1:50))
findDrawdowns(x)
table.Drawdowns(x)
chart.Drawdown(x)
Шейн
источник
Спасибо, Шейн, но, похоже, здесь можно найти локальные минимумы (или максимумы) - то есть одну точку в регионе. Мои данные (как и любые биологические данные) являются ШУМНЫМИ.> Мне не важны сами точечные минимумы, а большие регионы с низким уровнем.
Дэвид Б
Если у вас есть локальные максимальные и минимальные баллы, вы можете легко рассчитать разницу. Итак, вы хотите знать случаи, когда различия велики как по величине, так и по «продолжительности»? Это данные временного ряда?
Шейн
@david Возможно, вы можете использовать эту функцию итеративно. Используйте функцию для определения минимумов. Отбросьте эту точку и окружающие точки (скажем, x точек в пределах некоторого уровня допуска). Вы можете выбрать уровень допуска (например, + - 10 отсчетов), который будет определять плоскую область для вашего приложения. Найдите новые минимумы в новом наборе данных. Будет ли это работать?
@shane Аналогия, которая приходит на ум, - это долины в горном регионе. Я думаю, что цель состоит в том, чтобы идентифицировать все долины, и проблема в том, что некоторые долины «глубже», а некоторые «неглубокие» по отношению к горам.
@Shane Это не временные ряды, это координаты вдоль генома (хромосомы).
Дэвид Б.
2

Некоторые из пакетов Bioconductor (например, ShortRead , Biostrings , BSgenome , IRanges , genomeIntervals ) предлагают средства для работы с позициями генома или векторами покрытия, например, для ChIP-seq и идентификации обогащенных областей. Что касается других ответов, я согласен, что любой метод, основанный на упорядоченных наблюдениях с некоторым пороговым фильтром, позволил бы изолировать низкий сигнал в пределах определенной полосы пропускания.

Может быть, вы также можете посмотреть на методы, используемые для выявления так называемых "островов"

Zang, C, Schones, DE, Zeng, C, Cui, K, Zhao, K и Peng, W (2009). Кластерный подход для идентификации обогащенных доменов по данным модификации гистонов ChIP-Seq . Биоинформатика, 25 (15) , 1952-1958.

хл
источник