Нахождение минимального набора вырожденных последовательностей ДНК

#алгоритм #строки #биоинформатика #sof

Как всегда, небольшое упрощенное вступление:

В биологии в большинстве случаев для описания последовательности ДНК достаточно 4х
литер - A,G,T и C (это обозначения нуклеотидов: аденин, гуанин, тимин и цитозин). Для
специфических целей удобно также пользоваться литерами, обозначающими вырождение, например:

R = A или G, т.е. на данной позиции присутствует либо A, либо G, других вариантов нет
Y = T или C
W = A или T
V = A или C или G
N = A или C или G или T 


существуют обозначения для всех вариантов вырождений (вот ссылка на текстовый файл
со всеми вариантами)

Соответственно, к примеру, когда мы говорим о последовательности RAGCTY, это значит,
что мы имеем в виду варианты AAGCTT или AAGCTC или GAGCTT или GAGCTC. Это прямая задача
- нахождение всех вариантов вырожденной строки. Она проста.

Но есть обратная задача. Очень часто для исследований/диагностики требуется набор
коротких ДНК-последовательностей (праймеров) (обычно, от 8 до 30 нуклеотидов). Заказ
каждой последовательности стоит времени/денег, поэтому я могу, к примеру, рассчитывать,
что мою заявку на 50 вариантов удовлетворят, а вот больше - уже откажут. При этом,
если я вижу, что две или более последовательности составляют вырожденность (как в вышеуказанном
примере), я стараюсь подать в заявке именно вырожденный вариант. Синтез по времени
и деньгам занимает практически столько же, зато я имею 2 или более последовательностей
на "одноразовый заказ" для своего исследования.

Это просто, когда последовательностей мало и они короткие. Но порой (часто такое
бывает при исследовании нового вируса или сильно измененного штамма уже известного) имеем:

Набор большого количества последовательностей (строк) одной длины, состоящих из литер
A, G, C, T.
Задача: найти минимальное множество вырожденных последовательностей, описывающее
именно этот набор. 

Т.е., допустим, на входе  AAGCTT + AAGCTC + GAGCTT. Решение: RAGCTT + AAGCTC
Неправильные решения: (NAGCTT + AAGCTC) или NAGCTN , поскольку в первом случае мы
имеем 5 последовательностей, а во втором - 16, тогда как у нас их всего 3, и в дальнейшей
работе (как теоретической, так и практической мы будем иметь неправильные последовательности).

В качестве примера хорошей (возможно, идеальной) работы: набор реальных невырожденных
праймеров (отсортированы, в количестве 94 штук) и результат после поиска вырождений
(14 вырожденных праймеров)

Это просто, когда последовательностей хотя бы несколько десятков. Но когда речь идёт
о сотнях тысяч (дальше будет отбор по другим критериям, которые, опять-таки, лучше
работают с вырожденностями), то решение неочевидно.
Существующие пакеты предлагают решение для относительно небольших наборов. Алгоритм
в них достаточно прост: первая строка попарно сравнивается со всеми остальными, при
нахождении разницы в одну литеру обе строки удаляются из набора, заменяясь одной с
вырождением по этой литере. Алгоритм неоптимален, и часто потом удаётся вручную слегка
улучшить результат. Плюс, он не работает с большими наборами.

Можно ускорить этот алгоритм путём перевода строк в полубайтовые массивы: 
A = 0001, C = 0010, G = 0100, T = 1000. Тогда, к примеру, R = 0101, V = 0111, N=
1111 и т.п. Это должно весьма сильно ускорить сравнение при переборе, но, думаю, серьезного
ускорения не добавит.

Есть идеи других алгоритмов решения обратной задачи?

Если кому интересен реальный материал - 2 миллиона 23-буквенных праймеров (в будущем
проекте их больше -  более, чем в 2 раза, но и так файл весит 14 мегабайтов).
    

Ответы

Ответ 1

Допустим у нас есть набор:

AAA
AAT
ATA
CAA
GAA


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

AAA
CAA
GAA
ATA
AAT


То, найдя такую серию, можно будет её смерджить в более короткую

VAA
ATA
AAT


И если проделать такую операцию сначала для всех серий для первого символа, потом
для всех для второго и т.д., то в итоге получим что то сжатое. Корректность алгоритма
я навряд ли докажу, но пока не смог найти контрпримеров. 

Итак, идем к коду. 

Для начала я добавил перечисление для типов

[Flags]
public enum Types
{
    None = 0,
    A = 1,
    C = 1 << 1,
    G = 1 << 2,
    T = 1 << 3,
    N = A | C | G | T,
    V = A | C | G,
    H = A | C | T,
    D = A | G | T,
    B = C | G | T,
    Y = T | C,
    M = A | C,
    K = G | T,
    S = G | C,
    W = A | T,
    R = A | G
}


Потом добавил методы для сопоставления типов и символов

internal static class Exts
{


В одну сторону

    public static char GetChar(this Types t)
    {
        switch (t)
        {
            case Types.A: return 'A';
            case Types.C: return 'C';
            case Types.G: return 'G';
            case Types.T: return 'T';
            case Types.N: return 'N';
            case Types.V: return 'V';
            case Types.H: return 'H';
            case Types.D: return 'D';
            case Types.B: return 'B';
            case Types.Y: return 'Y';
            case Types.M: return 'M';
            case Types.K: return 'K';
            case Types.S: return 'S';
            case Types.W: return 'W';
            case Types.R: return 'R';
            default: return '?';
        }
    }


В обратную сторону        

    public static Types GetTypes(this char c)
    {
        switch (c)
        {
            case 'A': return Types.A;
            case 'C': return Types.C;
            case 'G': return Types.G;
            case 'T': return Types.T;
            case 'N': return Types.N;
            case 'V': return Types.V;
            case 'H': return Types.H;
            case 'D': return Types.D;
            case 'B': return Types.B;
            case 'Y': return Types.Y;
            case 'M': return Types.M;
            case 'K': return Types.K;
            case 'S': return Types.S;
            case 'W': return Types.W;
            case 'R': return Types.R;
            default: return Types.None;
        }
    }
}


После я добавил слегка модифицированную сортировку подсчетом

public class CountingSort
{
    public void Sort(string[] a, int w, int shift, int n, string[] aux)
    {
        int R = 21;

        for (int d = 0; d < w; d++)
        {
            var ind = (shift + d) % w;

            int[] count = new int[R + 1];
            for (int i = 0; i < n; i++)
                count[GetInd(a[i], ind) + 1]++;

            for (int r = 0; r < R; r++)
                count[r + 1] += count[r];

            for (int i = 0; i < n; i++)
                aux[count[GetInd(a[i], ind)]++] = a[i];

            for (int i = 0; i < n; i++)
                a[i] = aux[i];
        }
    }


Ну и получение значения по индексу я особо не продумывал, сделал как сделал. Знаю,
что можно было переиспользовать код выше, но меня это сильно не заботило

    private static int GetInd(string s, int p)
    {
        if (s == null) return 20;
        var c = s[p];

        switch (c)
        {
            case 'A': return 1;
            case 'C': return 2;
            case 'G': return 3;
            case 'T': return 4;

            case 'N': return 15;
            case 'V': return 14;
            case 'H': return 13;
            case 'D': return 12;
            case 'B': return 11;
            case 'Y': return 10;
            case 'M': return 9;
            case 'K': return 8;
            case 'S': return 7;
            case 'W': return 6;
            case 'R': return 5;
        }

        return 0;
    }
}


В итоге, добавил класс для декомпрессии

public class DnaCompressor
{


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

    public void Compress(string[] lines, IProgress progress)
    {
        var tmp = new string[lines.Length];
        int n = lines.Length;
        var maxpos = lines[0].Length;
        var sorter = new CountingSort();

        for (var pos = 0; pos < maxpos; pos++)
        {
            sorter.Sort(lines, lines[0].Length, pos, n, tmp);
            var ind = 0;

            while (ind < lines.Length && lines[ind] != null)
            {
                int start = ind;
                while (ind < (lines.Length - 1) && lines[ind + 1] != null && IsSame(lines[ind],
lines[ind + 1], pos)) ind++;
                int end = ind;

                Types t = Types.None;
                for (int i = start; i <= end; i++)
                    t |= lines[i][pos].GetTypes();

                var added = 0;
                foreach (var tt in new[]
                {
                    Types.N, Types.V, Types.H, Types.D, Types.B, Types.Y, Types.M,
Types.K, Types.S, Types.W,
                    Types.R, Types.A, Types.C, Types.G, Types.T
                })
                {
                    if (t == Types.None) break;

                    if ((t & tt) == tt)
                    {
                        t &= ~tt;
                        var sb = new StringBuilder(lines[start]) { [pos] = tt.GetChar() };
                        lines[start] = sb.ToString();
                        added++;
                    }
                }

                ind++;
                for (var i = start + added; i <= end; i++)
                    lines[i] = null;
            }

            n = Scroll(lines, n);

            progress.Report((100 * pos) / maxpos);
        }

        progress.Report(100);
    }


Этот метод просто двигает пустые строки в конец массива

    internal static int Scroll(string[] lines, int n)
    {
        int startInd = 0; int endInd = n - 1;
        while (startInd < endInd)
        {
            while (startInd < endInd && lines[startInd] != null) startInd++;
            while (startInd < endInd && lines[endInd] == null) endInd--;
            if (startInd < endInd)
            {

                lines[startInd] = lines[endInd];
                lines[endInd] = null;
            }

        }
        n = endInd + 1;
        return n;
    }


Проверяем, все ли символы в строке одинаковые, за исключением указанной позиции

    internal static bool IsSame(string s1, string s2, int except)
    {
        for (var i = 0; i < s1.Length; i++)
        {
            if (i == except) continue;
            if (s1[i] != s2[i]) return false;
        }
        return true;
    }
}


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

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



Все исходники доступны тут

Бинарник потыкать можно скачать тут