Вроде как начало дискуссии :)

Спасибо за чат :)

А теперь вопрос: длинная последовательность N в скаффолде часто означает невозможность точного предсказания нуклеотиодного контента (повторы и т.п.). А нельзя здесь как-то сэкономить? не важно, в памяти или в скорости

Вот это как раз вопрос оптимизации. Короткий неизвестный участок в референсе (насколько короткий, - задаете уже вы) должен рассматриваться, как последовательность ЛЮБЫХ нуклеотидов. Длинный же участок должен, в принципе, изыматься из рассмотрения. Применительно к нашей ситуации, вы должны задуматься, какой % совпадения вас устраивает. Исходя из этого, и пляшете.

Но есть ещё один тонкий момент: ведь последовательность N может быть длинной, а может быть короткой. И это, по-хорошему, разные ситуации в геноме. Поэтому дерево должно строиться по коротким неизвестным последовательностям, но пропускать длинные

@ViktorTomilov Понял. Таким образом, я могу определить, например, точные координаты повторов.

@AlexeyKozlov Да, пожалуй: если один конец заходит в определенную область, а другой - в неопределенную, то да, вы уточните референс.

@ViktorTomilov Не понял в чем проблема комплиментарного поиска. зачем идти по дереву задом наперед. ведь можно заранее развернуть и перекодировать паттерн в обычный вариант. даже если паттерн поступает на вход потоком, можно накопить какое то количество букв и развернуть

И как предполагается при индексации писать N. Если N в паттерне то понятно, надо просматривать все 4 направления. Но в дереве тогда похоже понадобится пятая ветка

@Mike Можно. Собственно, в стандартном подходе и используют 2 вариант: 1. Сделать палиндром референса и искать комплементарные варианты по нему; 2. делать палиндром паттерна и возвращать результат поиска как бы "с конца"

все равно не понял. зачем вообще нужен палиндром референса. ведь из палиндромного паттерна можно сделать обычный и искать по обычному референсу

@Mike Если N - немного (насколько немного, тут решает сам бионформатик. У нас немного = не более 50% N на паттерн), то они входят во все children. Если же МНОГО, то просто не индексируются: есть обозначение "конец значимой цепочки и ссылка на первый значимый children. Таким образом экономится и место, и время

@Mike для нескольких раз - да. А теперь представьте, что ваша программа постоянно висит в памяти и обслуживает запросы поиска по разным паттернам на один референс. Экономичнее сделать один раз комплементарный референс и держать его в памяти, нежели, допустим, 100000 раз делать комплемент паттерну.

@ViktorTomilov Ага, т.е. запрашивающий дает просто паттрен, а мы должны искать как прямой вариант так и комплимент ?

и заранее что искать не известно

@Mike Да. В живом организме, фактически, не всегда понятно, какая из цепей ДНК будет активна, т.е. с какой стороны будет работать фермент. Разумеется, пользователь может сам задать два варианта, но такие случаи редки :)

интересна еще статистика текста. скажем если взять 20 буквенные последовательности, то велик ли разброс в их частоте. Есть ли такие, которые встречаются в разы чаще, чем другие

@Mike Если брать ДНК в целом, то да. Есть такие части, называемые повторы, роль которых мы только начинаем понимать, - они разбросаны по всему геному и очень похожи друг на друга. Например, система называемая "Alu repeats" включает в себя более миллиона похожих друг на друга паттернов, входящих в ДНК любого человека.

@Mike Если же брать только значимую часть ДНК, то там такая статистика маловстречаема

@ViktorTomilov Это я к тому, что если мы сделаем индекс на 20 букв, нам его не хватит что бы потом просто дочитать хвосты уже по оставшемуся списку вхождений. хотя можно конечно смотреть много ли раз встречается и если часто то делать именно в этом месте дерево глубже. Кроме того разбив входной паттерн на несколько 20 буквенных частей можно получить списки вхождений для каждого из них и сравнить их, добавляя требуемое смещение

И да, я же правильно понимаю, что цель дерева не только дать информацию "такое встречается", но и указать точное смещение в референсе ?

@Mike не могу точно ответить, просто не делал подобного. Дело в том, что повторяющиеся элементы в ДНК на данном этапе не интересны - мы почти не можем предсказывать и диагностировать заболевания на основании их анализа. Для диагностики мы анализируем значимые участки (гены, промоторы, CpG-островки), которые друг на друга, в общем-то, не похожи (кроме, может быть, небольших участков у некоторых из них)

@Mike Да, разумеется, в этом-то и вся цель

@ViktorTomilov Это я как раз к тому, что может быть можно хранить гораздо меньший индекс, если в нем будет меньше букв. там же все равно хранение идет так что бы можно было найти произвольный участок. в конце каждой цепочки у нас висят списки смещений (от этого похоже не уйти). списки отсортированы (сама индексация идет подряд, так что списки изначально отсортированы). Сравнение двух и более отсортированных списков проходит за один проход по ним всем одновременно.

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

@ViktorTomilov И где нибудь можно пример такого файла взять. хочу статистику посчитать, может на основе нее еще какие идеи придут ...

последний вопрос отпадает. нашел уже :)

формат Fasta, правда я не понял чем там отличаются маленькие буквы от больших. буду считать все в одном регистре

@Mike да, fasta подойдёт. Маленькие буквы значат, что "не 100% точно установлено, но скорее всего это так". У вас, судя по всему, старый релиз, но это не важно - переведите все в Uppercase, для общего представления подойдет

@Mike Но если вдруг захотите последний релиз, то вот он: ncbi.nlm.nih.gov/grc/human/data - кликаете на столбец GenBank accession напротив нужной хромосомы, на следующей странице - Send to, выбираете file, формат FASTA и кнопка Create file

@ViktorTomilov Я уже начал смотреть, любопытная картина. кстати буквы N в том что я взял встречаются редко. я пока эти места вообще пропускал. составил словрь частот 9и буквенных последовательностей. получил 262123 варианта, при том, что максимально возможно комбинаций из 9 букв - 262144. Т.е. пытаться строить структуры со ссылками друг на друга последовательно по 1 букве даже не стоит и пытаться, потому что покрыто все поле возможных вариантов.

Проще плоскую таблицу в которую можно обращаться сходу вычислив смещение по 9 символам входной строки.

@Mike Есть GC-богатые регионы, есть AT-богатые - думаю, термины пояснять не надо. Их наличие зависит от функционала и организма. То есть, у какого-нибудь гриба может быть перекос в одну сторону, а у клевера - в другую. :)

дошел до 11 буквенных ... их почти 4 миллиона вариантов. но на 266к меньше, чем теоретический предел. т.е. начали наблюдаться пропуски (6%). Даже в этом случае в дереве все еще нет смысла. Начиная с 12 символов оно может быть уже интересно. И это же значит что последовательности короче 11 букв вообще нет особого смысла искать, если конечно какой то конекретный вариант не особо уникален. но я так понял такие короткие и не ищут особо

@Mike Когда планируется эксперимент по нахождению определенного фрагмента ДНК, мы стараемся подобрать праймеры (уникальные короткие последовательности ДНК, которые "прилипнут" только к определенному участку и ни к какому другому) длиной не менее 18 букв, а что и 24 буквы.

*бывает, что

@ViktorTomilov Вы в update2 слишком сильно изменили задачу. Если ищем что то с небольшим алфавитом, но почти полным покрытием возможных перестановок букв на ограниченной длине ключа, то дерево не подходит совершенно, накладные расходы слишком большие. Если же алфавит велик, то покрытие возможных вариантов даже на очень коротких ключах порядка 5-6 букв очень мало и тогда нужны как раз деревья. какие конкретно деревья это уже другой вопрос.

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

@Mike Эмм... Я не менял задачу, - это просто мысль вслух, попытка объяснить, почему побитовые представления так и не появились в бесплатных реализациях

@ViktorTomilov Интересно, какова вероятность, что в паттерне при его длине в 100 букв вообще не будет редких 11 буквенных участков, встречающихся скажем менее 1000 раз. т.е. он весь будет состоять только из чего то очень распространненного. А то я все больше склоняюсь к мысли брать по 10-11 букв, переводить их непосредственно в смещение в таблице указателей на списки подходящих участков и сверять между собой для начала самые короткие из этих списков.

если встретится 1-3 редких участков списки из 100 значений сравнить и найти общее будет очень быстро. потом уже принимать решение, либо сразу проверять исходник по найденным нескольким адресам, либо сравнивать с более крупными списками

При этом индексацию я предполагаю выполнить простыми средствами за 2 прохода по исходнику.

@Mike честно говоря, скорее будет наоборот - будет более редкое сочетание, чем в среднем по хромосоме. Часто попадающие участки, за редким исключением, нужны в специфических исследованиях.

@Mike Каким образом? Хэшированием?

@ViktorTomilov Нет. С хешами нельзя даже 1 букву N обработать. Нам нужны 2 массива. один собственно индекс, скажем на 4 млн. элементов для 11 буквенного ключа. Элемент из 2 частей. Кол-во вхождений 4 байта и смещение в списке еще 4 байта. Второй массив собственно список, он большой, я так понимаю понадобится 4*размер исходной последовательности. берем 11 букв паттерна из них делаем адрес, буквально кодируем каждую букву 2 битами

(для других алфавитов представляем его как число в системе счисления с основанием в размер алфавита)

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

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

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

хотя думаю можно попробовать найти участки без букв N

при индексации думаю, что бы по 100 раз не перераспределять память выделяя кускам и и потом еще копируя если не влезает, делаем 2 прохода. Сначала считаем сколько раз встретилась та или иная последовательность. потом выделяем сразу один кусок памяти подо все. И за второй проход по исходнику наполняем списки. которые лежат в этом одном массиве просто друг за другом

при чтении исходника разумеется берем ключи со сдвигом на 1 букву

@Mike Если N много подряд, их можно вообще не брать в расчёт, просто выкидывать из рассмотрения, не забывая, однако, сдвигать смещение на их количество. "много" - это 70% от размера паттерна и более.

кстати я посмотрел исходник что вы дали. там одиночных N вообще нет, все идут блоками

@Mike последний релиз в этом цикле, видимо, исправили все недочеты. Я сам еще с последним не работал

ну вот вопрос, что если в раттерне через одну букву все N. тогда я хз как искать. вариантов слишком много

Да, и как вариант можно сделать такую таблицу на меньшую длину ключа например, но на ее концы повесить не списки, а корни деревьев продолжающих последовательнось.

тут главное не городить дерево с самого верха. кода бывают все возможные комбинации. а то только и будете что друг за другом считайте располагать по 4 структуры для каждого направления и тут же следующие уже 16 и все снабжать указателями друг на друга. И это при том, что на первых 10 уровнях дерева будет меньше тясячи пустых указателей, т.е. несуществующих последовательностей

@Mike Нет, такие паттерны - это моветон: либо прибор накосячил, либо руки кривые. А вот паттерн, где по краям значащие участки, а середина - NNNN, - нормальная ситуация

Кстати, если биты адресов для букв правильно подобрать можно получать палиндром буквально одной инверсией бит (xor'ом) и одним разворотом задом на перед. правда я что то не помню у x86 разворота задом наперед, по 1 биту выдвигать-вдвигать через перенос не хочется