Я изучаю некоторые данные о геномном покрытии, которые в основном представляют собой длинный список (несколько миллионов значений) целых чисел, каждый из которых говорит о том, насколько хорошо (или «глубоко») охвачена эта позиция в геноме.
Я хотел бы найти «долины» в этих данных, то есть регионы, которые значительно «ниже», чем их окружение.
Обратите внимание, что размер долин, которые я ищу, может варьироваться от 50 оснований до нескольких тысяч.
Какую парадигму вы бы порекомендовали использовать для поиска этих долин?
ОБНОВИТЬ
Некоторые графические примеры для данных:
ОБНОВЛЕНИЕ 2
Определение долины - это, конечно, один из вопросов, с которыми я борюсь. Это очевидные для меня:
но есть несколько более сложных ситуаций. В целом, есть три критерия, которые я рассматриваю: 1. (Среднее? Максимальное?) Покрытие в окне по отношению к глобальному среднему. 2. Покрытие (...) в окне относительно его непосредственного окружения. 3. Насколько велико окно: если я вижу очень низкий охват для короткого промежутка, это интересно, если я вижу очень низкий охват для длинного промежутка, это также интересно, если я вижу слегка низкое покрытие для короткого промежутка, это не очень интересно , но если я вижу слегка низкий охват для длинного промежутка - это .. Так что это сочетание длины sapn и его покрытия. Чем дольше, тем выше я оставляю покрытие и все равно считаю его долиной.
Спасибо,
Дейв
Ответы:
Вы можете использовать своего рода подход Монте-Карло, например, используя скользящее среднее ваших данных.
Возьмите скользящее среднее данных, используя окно разумного размера (я думаю, что вам решать, насколько широко).
Пропуски в ваших данных будут (конечно) характеризоваться более низким средним, поэтому теперь вам нужно найти какой-то «порог», чтобы определить «низкий».
Для этого вы случайным образом меняете значения ваших данных (например, используя
sample()
) и пересчитываете скользящее среднее для ваших замененных данных.Повторите этот последний отрывок достаточно много раз (> 5000) и сохраните все средние значения этих испытаний. По сути, у вас будет матрица из 5000 строк, по одной на пробу, каждая из которых содержит скользящую среднюю для этого испытания.
В этот момент для каждого столбца вы выбираете квантиль 5% (или 1% или что вы хотите), то есть значение, под которым лежит только 5% от среднего значения рандомизированных данных.
Теперь у вас есть «доверительный предел» (я не уверен, что это правильный статистический термин) для сравнения ваших исходных данных. Если вы обнаружите, что часть ваших данных ниже этого предела, вы можете назвать это сквозным.
Конечно, следует помнить, что ни этот, ни какой-либо другой математический метод никогда не могли бы дать вам никаких указаний на биологическую значимость, хотя я уверен, что вы хорошо это знаете.
РЕДАКТИРОВАТЬ - пример
Это только позволит вам графически найти регионы, но вы можете легко найти их, используя что-то на линии
which(values>limits.5)
.источник
Я совершенно не знаю этих данных, но, предполагая, что данные упорядочены (не по времени, а по положению?), Имеет смысл использовать методы временных рядов. Существует множество методов идентификации временных кластеров в данных. Обычно они используются для поиска высоких значений, но могут использоваться для низких значений, сгруппированных вместе. Я имею в виду статистику сканирования, статистику суммарной суммы (и другие), используемые для обнаружения вспышек заболеваний в данных подсчета. Примеры этих методов находятся в пакете наблюдения и пакете DCluster.
источник
surveillance
иDCluster
, но не могли бы вы быть более конкретным? Они оба относительно большие пакеты, и их цель кажется довольно конкретной. Я не уверен, с чего начать.Есть много вариантов для этого, но один хороший: вы можете использовать
msExtrema
функцию вmsProcess
пакете .Редактировать:
В анализе финансовых показателей этот вид анализа часто выполняется с использованием концепции «просадки».
PerformanceAnalytics
Пакет имеет некоторые полезные функции , чтобы найти эти долины . Вы можете использовать тот же алгоритм здесь, если вы рассматриваете свои наблюдения как временные ряды.Вот несколько примеров того, как вы могли бы применить это к своим данным (где «даты» не имеют значения, но используются только для упорядочивания), но первыми элементами в
zoo
объекте будут ваши данные:источник
Некоторые из пакетов Bioconductor (например, ShortRead , Biostrings , BSgenome , IRanges , genomeIntervals ) предлагают средства для работы с позициями генома или векторами покрытия, например, для ChIP-seq и идентификации обогащенных областей. Что касается других ответов, я согласен, что любой метод, основанный на упорядоченных наблюдениях с некоторым пороговым фильтром, позволил бы изолировать низкий сигнал в пределах определенной полосы пропускания.
Может быть, вы также можете посмотреть на методы, используемые для выявления так называемых "островов"
источник