如何按升序对最常见的 3 聚体列表进行排序?

How to sort my list of most frequent 3-mers in ascending order?

我正在编写代码以找出 DNA 序列中最常见的 3 聚体。我写了一个代码来计算 3-mer 的出现次数,如果它大于 1,那么代码会同时记录字符串和出现次数。

这给了我一个本质上多余的列表。我想对列表进行排序,以便我只在列表中看到每个 3-mer 一次。

下面是写的代码:

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

    printf("Enter the DNA sequence\n");
    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] && dna[j+1] == read[1] && dna[j+2] == read[2])
            {
                count++;
            }
            else
            {
                continue;
            }
        }
        if(count > maxcount)
        {
            maxcount = count;
            
           printf("%s %d\n", 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 sequence\n");
    char* res = fgets(dna, sizeof(dna), stdin);
    if (res == NULL)
        return 1;
    char* ptr = &dna[0];
    for (int i = 0; i < strlen(dna) - 2; i++)
        insertTrimer(&d, ptr++);
    for (int i = 0; i < d.count; i++)
        printf("%s : %d\n", d.trimer[i].n_base, d.trimer[i].count);
    return 0;
}

基本上,每个 3 聚体都是更大结构中的一个条目。较大的结构是 binary-searched,每次发现新的 3-mer 时 q-sorted。否则,如果发现重复,则他们的条目会递增。

这是您输入的结果

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 这样的程序
  • 使用哈希图。 hashmaps/tables 没有标准 C-library。你基本上要做的事情与我在这里所做的非常相似。内存管理可能是一个挑战。但是,您将对序列中的每个 3-mer 进行 O(1) 搜索,而不是 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 sequence\n");
    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,它会变得更快,但是,输出不会被排序。