Блог

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

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

Ниже приведен код, который написал:

 int main()
{
    char dna[1000];
    char read[3] = {0};
    char most_freq[3];

    printf("Enter the DNA sequencen");
    fgets(dna, sizeof(dna), stdin);
    int i, j;
    for(i=0; i<strlen(dna)-3; i  )
    {

        read[0] = dna[i];
        read[1] = dna[i 1];
        read[2] = dna[i 2];
        int count=0, maxcount=1;
        for(j = 0; j < strlen(dna); j  )
        {
            if(dna[j] == read[0] amp;amp; dna[j 1] == read[1] amp;amp; dna[j 2] == read[2])
            {
                count  ;
            }
            else
            {
                continue;
            }
        }
        if(count > maxcount)
        {
            maxcount = count;
            
           printf("%s %dn", read, maxcount);
        }

    }
}
  

Это то, что я получаю, если я ввожу :

CGCCTAAATAGCCTCGCGGAGCCTTATGTCATACTCGTCCT

 CGC 2
GCC 3
CCT 4
ATA 2
AGC 2
GCC 3
CCT 4
CTC 2
TCG 2
CGC 2
AGC 2
GCC 3
CCT 4
GTC 2
ATA 2
CTC 2
TCG 2
GTC 2
CCT 4
  

Понятно, что ответ — CCT, но я не хочу избыточности в выводе. Как мне решить эту проблему?

Вот достаточно быстрый способ сделать это в C

 #include <stdlib.h>
#include <string.h>
#include <stdio.h>

typedef struct {
    char n_base[4];
    int count;
} NMer_3;

typedef struct {
    int count;
    NMer_3 trimer[4 * 4 * 4];
} dict;

int cmp(const void* a, const void* b) {
    return strncmp((char*)a, (char*)b, 3);
}

void insertTrimer(dict* d, char c[3]) {
    NMer_3* ptr = (NMer_3*)bsearch((void*)c, (void*)d->trimer, d->count,
                                   sizeof(NMer_3), cmp);
    if (ptr == NULL) {
        int offset = d->count;
        strncpy(d->trimer[offset].n_base, c, 3);
        d->trimer[offset].count = 1;
        d->count  ;
        qsort(d->trimer, d->count, sizeof(NMer_3), cmp);
    } else {
        ptr->count  ;
    }
}

int main() {
    char dna[1000];
    dict d;
    printf("Enter the DNA sequencen");
    char* res = fgets(dna, sizeof(dna), stdin);
    if (res == NULL)
        return 1;
    char* ptr = amp;dna[0];
    for (int i = 0; i < strlen(dna) - 2; i  )
        insertTrimer(amp;d, ptr  );
    for (int i = 0; i < d.count; i  )
        printf("%s : %dn", d.trimer[i].n_base, d.trimer[i].count);
    return 0;
}
  

По сути, каждый 3-мер является записью в более крупной структуре. Большая структура выполняется с двоичным поиском и q-сортируется каждый раз, когда обнаруживается новый 3-мер. В противном случае, если найдено повторение, их запись увеличивается.

Вот результат, используемый с вашим вводом

 AAA : 1
AAT : 1
ACT : 1
AGC : 2
ATA : 2
ATG : 1
CAT : 1
CCT : 4
CGC : 2
CGG : 1
CGT : 1
CTA : 1
CTC : 2
CTT : 1
GAG : 1
GCC : 3
GCG : 1
GGA : 1
GTC : 2
TAA : 1
TAC : 1
TAG : 1
TAT : 1
TCA : 1
TCC : 1
TCG : 2
TGT : 1
TTA : 1
  

Способы повышения скорости:

• Используйте программу, подобную Jellyfish
• Используйте hashmap. Для хэш-карт / таблиц нет стандартной C-библиотеки. Вы в основном собираетесь сделать что-то очень похожее на то, что я сделал здесь. Управление памятью может быть проблемой. Однако вы собираетесь выполнять поиск O (1) для каждого 3-мер в последовательности вместо O (log (n)), кроме того, добавление будет только O (1) вместо необходимости сортировки O (n * log (n)) .

Если вы делаете это на C , вы получаете много преимуществ, первое из которых — гораздо более простой код:

 #include <string>
#include <iostream>
#include <map>

int main() {
    std::string dna;
    printf("Enter the DNA sequencen");
    std::getline(std::cin, dna);
    auto d = std::map<std::string,int>{};
    for (int i = 0; i < dna.size() - 2; i  ){
        auto mer3 = dna.substr(i,3);
        auto itr = d.find(mer3);
        if (itr == d.end()){
            d[mer3] = 1;
        } else {
            itr->second  = 1;
        }
    }
    for (auto i : d) std::cout << i.first << ':' << i.second << 'n';
    std::cout <<std::endl;
    return 0;
}
  

Это фактически то же самое, что и в примере C.

map Однако, если вы замените unordered_map его, он станет намного быстрее, вывод не будет отсортирован.