Écrivez une idée

candidates = [
    'GGGAGCAGGCAAGGACTCTG',
    'GCTCGGGCTTGTCCACAGGA',
    '...',
    # 被你看出来啦，这些其实人类基因的片段
]

bg_db = [
    'CTGCTGACGGGTGACACCCA',
    'AGGAACTGGTGCTTGATGGC',
    '...',
    # 这个更多，有十亿左右
]

Parce que vos données sont en réalité très uniques, elles peuvent être rationalisées ici.
Parce que toutes les chaînes comportent 20 caractères et se composent de ATCG quatre caractères. Ensuite, ils peuvent être convertis en nombres entiers à des fins de comparaison.
La représentation binaire est la suivante

A  ==>  00
T  ==>  01
C  ==>  10
G  ==>  11

Étant donné que la longueur d'une chaîne est fixe et que chaque caractère peut être représenté par 2 bits, chaque chaîne peut être représentée par un entier 40 bits. Il peut être exprimé sous la forme de 32+8, ou vous pouvez utiliser directement la mise en forme de bits 64. Il est recommandé d'utiliser le langage C pour ce faire.

Parlons de comparaison.
Parce que nous devons trouver 每一个candidate在bg_db中与之差异小于等于4的所有记录, donc tant que nous effectuons l'opération ^ XOR au niveau du bit sur deux entiers, le résultat 二进制中1不超过8个,且这不超过8个1最多只能分为4个组 peut répondre aux exigences (00^11 = 11 ,10^01=11).
Divisez les 40 bits du résultat en 20 groupes, c'est-à-dire qu'il y a au plus des 4 groupes avec les trois valeurs b01 b10 b11, et le reste est tous b00 .
Ensuite, l'algorithme de comparaison est facile à écrire.
Vous pouvez obtenir plusieurs groupes de trois valeurs non nulles pour chaque octet (quatre groupes) pour obtenir brièvement le résultat global de la comparaison.
Parce qu'il n'y a que des 256 valeurs possibles pour chaque octet, et qu'il n'y a que des 3^4=81 types de valeurs qualifiées , le résultat peut d'abord être stocké puis récupéré.
Voici une fonction pour obtenir combien de groupes non nuls il y a dans le résultat.

/*****************下面table中值的生成******//**
  int i;
  for( i=0;i<256;++i){
    int t =0;
    t += (i&0x01 || i&0x02)?1:0;
    t += (i&0x04 || i&0x08)?1:0;
    t += (i&0x10 || i&0x20)?1:0;
    t += (i&0x40 || i&0x80)?1:0;
    printf("%d,",t);
    if(i%10 ==9){putchar('\n');}
  }
********************************************//

int table[] = {
0,1,1,1,1,2,2,2,1,2,
2,2,1,2,2,2,1,2,2,2,
2,3,3,3,2,3,3,3,2,3,
3,3,1,2,2,2,2,3,3,3,
2,3,3,3,2,3,3,3,1,2,
2,2,2,3,3,3,2,3,3,3,
2,3,3,3,1,2,2,2,2,3,
3,3,2,3,3,3,2,3,3,3,
2,3,3,3,3,4,4,4,3,4,
4,4,3,4,4,4,2,3,3,3,
3,4,4,4,3,4,4,4,3,4,
4,4,2,3,3,3,3,4,4,4,
3,4,4,4,3,4,4,4,1,2,
2,2,2,3,3,3,2,3,3,3,
2,3,3,3,2,3,3,3,3,4,
4,4,3,4,4,4,3,4,4,4,
2,3,3,3,3,4,4,4,3,4,
4,4,3,4,4,4,2,3,3,3,
3,4,4,4,3,4,4,4,3,4,
4,4,1,2,2,2,2,3,3,3,
2,3,3,3,2,3,3,3,2,3,
3,3,3,4,4,4,3,4,4,4,
3,4,4,4,2,3,3,3,3,4,
4,4,3,4,4,4,3,4,4,4,
2,3,3,3,3,4,4,4,3,4,
4,4,3,4,4,4
};

int getCount(uint64_t cmpresult)
{
    uint8_t* pb = &cmpresult;    // 这里假设是小段模式，且之前比较结果是存在低40位
    return table[pb[0]]+table[pb[1]]+table[pb[2]]+table[pb[3]]+table[pb[4]];
}

Tout d'abord, votre estimation du temps est complètement fausse. Ce type de traitement de données à grande échelle doit exécuter des dizaines de milliers d'éléments et durer plus de dix secondes avant de pouvoir être utilisé pour la multiplication afin de calculer le temps total. Si un seul élément est compté, cela Le temps représente presque toute la surcharge du processus d'initialisation, plutôt que la surcharge critique d'E/S et de CPU

Le texte suivant

Les quatre possibilités d'ACTG sont équivalentes à 2 bits. C'est trop inutile d'utiliser un caractère pour représenter une position de gène. Un caractère fait 8 bits et peut contenir 4 positions de gène

Même si vous n'utilisez aucun algorithme et écrivez simplement vos 20 gènes sous forme binaire, vous pouvez gagner 5 fois plus de temps

De plus, après avoir bouclé 20 fois, le nombre d'instructions CPU est de 20*n, et n est estimé à au moins 3. Mais pour le binaire, l'opération XOR de comparaison est directement une instruction CPU, et le nombre de les instructions sont 1

Je ne connais pas grand-chose à l'algorithme, mais par expérience, les algorithmes complexes prennent plus de temps et ne sont pas aussi rapides que cet algorithme simple et grossier

Vous pouvez envisager le multithreading et le clustering pour traiter les données

Au fait, il semble que la distance de Hamming puisse également être calculée

Aussi aucun algorithme n'est utilisé, solution de force brute, écrite en C

Sur ma machine (CPU : Core 2 Duo E7500, RAM : 4G, OS : Fedora 19), résultats des tests

candidates    bg.db        cost
10000    1000000    318758165微秒
500      1000000    14950302微秒 

Si vous passez au processeur à 24 cœurs du sujet, les performances seront améliorées de 20 fois, puis ajoutez 48 machines pour fonctionner ensemble, les opérations de 5 000 W prendront 15 secondes et le temps sera de
1 000 000 * 1 000 000 000. / 500 / 1000000 * 15 / 20 / 48 / 3600 / 24 = 3,616898 jours

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdbool.h>
#include <sys/time.h>

#define START_CC(flag)  \
    struct timeval st_##flag##_beg; \
    gettimeofday(&st_##flag##_beg, 0)

#define STOP_CC(flag)   \
    struct timeval st_##flag##_end; \
    gettimeofday(&st_##flag##_end, 0)

#define PRINT_CC(flag)  \
    double cost_##flag = 0.0L; \
    cost_##flag = (double)(st_##flag##_end.tv_sec - st_##flag##_beg.tv_sec); \
    cost_##flag = cost_##flag * 1000000L + (double)(st_##flag##_end.tv_usec - st_##flag##_beg.tv_usec);    \
    printf(#flag" cost time %.6lf microsecond.\n", cost_##flag);

#define GENEORDER_CODE_LENGTH 20 + 1

typedef struct GeneOrder
{
    char code[GENEORDER_CODE_LENGTH];
}GeneOrder, *GeneOrderPtr;

typedef struct GOArray
{
    size_t          capacity;
    size_t          length;
    GeneOrderPtr    data;
}GOArray;

GOArray createGOAarray(size_t capacity)
{
    GOArray goa;

    goa.capacity    = capacity;
    goa.length      = 0;
    goa.data        = (GeneOrderPtr)malloc(capacity * sizeof(GeneOrder));

    return goa;
}

void destroyGOArray(GOArray* goa)
{
    if (goa->capacity > 0) {
        free(goa->data);
    }
}

bool readGOFile(char const *file, GOArray *goarray)
{
    FILE* fp = NULL;

    if ((fp = fopen(file, "r+")) == NULL) {
        return false;
    }

    char buff[64];

    while (fgets(buff, 64, fp) != NULL) {
        if (goarray->length < goarray->capacity) {
            memcpy(goarray->data[goarray->length].code,
                buff,
                GENEORDER_CODE_LENGTH * sizeof(char)
            );
            goarray->data[goarray->length].code[GENEORDER_CODE_LENGTH - 1] = '
gcc -Wall -Wextra -o ccs main.c -std=c99 -Os && ./ccs candidate.list bg.db
';
            goarray->length ++;
        } else {
            fclose(fp);
            return true;
        }
    }

    fclose(fp);
    return true;
}

int main(int argc, char* argv[])
{
    (void)argc;

    GOArray condgo  = createGOAarray(10000);
    GOArray bggo    = createGOAarray(1000000);

    printf("loading ...\n");

    START_CC(loading);
    if (!readGOFile(argv[1], &condgo) || !readGOFile(argv[2], &bggo)) {
        destroyGOArray(&condgo);
        destroyGOArray(&bggo);
        return -1;
    }
    STOP_CC(loading);
    PRINT_CC(loading);


    int count = 0;

    START_CC(compare);
    for (size_t i = 0;i < 500;i ++) {
        const GeneOrderPtr gop = condgo.data + i;
        for (size_t j = 0;j < bggo.length;j ++) {
            const GeneOrderPtr inner_gop = bggo.data + j;
            int inner_count = 0;

            for (size_t k = 0;k < 20;k ++) {
                if (gop->code[k] != inner_gop->code[k]) {
                    if (++inner_count > 4) {
                        break;
                    }
                }
            }

            if (inner_count <= 4) {
            #ifdef DBGPRINT
                printf("%d %s - %s\n", i, gop->code, inner_gop->code);
            #endif
                count++;
            }
        }
    }
    STOP_CC(compare);
    PRINT_CC(compare);

    printf("result = %d\n", count);

    destroyGOArray(&condgo);
    destroyGOArray(&bggo);

    return 0;
}

Compiler les paramètres et exécuter

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <stdbool.h>
#include <sys/time.h>
#include <pthread.h>

#define START_CC(flag)              \
    struct timeval st_##flag##_beg; \
    gettimeofday(&st_##flag##_beg, 0)

#define STOP_CC(flag)               \
    struct timeval st_##flag##_end; \
    gettimeofday(&st_##flag##_end, 0)

#define PRINT_CC(flag)                                                                                  \
    double cost_##flag = 0.0L;                                                                          \
    cost_##flag = (double)(st_##flag##_end.tv_sec - st_##flag##_beg.tv_sec);                            \
    cost_##flag = cost_##flag * 1000000L + (double)(st_##flag##_end.tv_usec - st_##flag##_beg.tv_usec); \
    printf(#flag " cost time %.6lf microsecond.\n", cost_##flag);

#define GENEORDER_CODE_LENGTH 20 + 1

typedef struct GeneOrder {
    char code[GENEORDER_CODE_LENGTH];
} GeneOrder, *GeneOrderPtr;

typedef struct GOArray {
    size_t capacity;
    size_t length;
    GeneOrderPtr data;
} GOArray;

GOArray createGOAarray(size_t capacity)
{
    GOArray goa;

    goa.capacity = capacity;
    goa.length = 0;
    goa.data = (GeneOrderPtr)malloc(capacity * sizeof(GeneOrder));

    return goa;
}

void destroyGOArray(GOArray* goa)
{
    if (goa->capacity > 0) {
        free(goa->data);
    }
}

bool readGOFile(char const* file, GOArray* goarray)
{
    FILE* fp = NULL;

    if ((fp = fopen(file, "r+")) == NULL) {
        return false;
    }

    char buff[64];

    while (fgets(buff, 64, fp) != NULL) {
        if (goarray->length < goarray->capacity) {
            memcpy(goarray->data[goarray->length].code, buff,
                GENEORDER_CODE_LENGTH * sizeof(char));
            goarray->data[goarray->length].code[GENEORDER_CODE_LENGTH - 1] = '
gcc -Wall -Wextra -o ccs main.c -std=c99 -O3 -lpthread && ./ccs candidate.list bg.db
';
            goarray->length++;
        }
        else {
            fclose(fp);
            return true;
        }
    }

    fclose(fp);
    return true;
}

typedef struct ProcessST {
    GOArray* pcond;
    GOArray* pbg;
    size_t beg;
    size_t end; // [beg, end)
} ProcessST;

void* processThread(void* parg)
{
    ProcessST* ppst = (ProcessST*)parg;
    GOArray* pcond = ppst->pcond;
    GOArray* pbg = ppst->pbg;
    int count = 0;

    for (size_t i = ppst->beg; i < ppst->end; i++) {
        const GeneOrderPtr gop = pcond->data + i;
        for (size_t j = 0; j < pbg->length; j++) {
            const GeneOrderPtr inner_gop = pbg->data + j;
            int inner_count = 0;

            for (size_t k = 0; k < 20; k++) {
                if (gop->code[k] != inner_gop->code[k]) {
                    if (++inner_count > 4) {
                        break;
                    }
                }
            }

            if (inner_count <= 4) {
#ifdef DBGPRINT
                printf("%d %s - %s\n", i, gop->code, inner_gop->code);
#endif
                count++;
            }
        }
    }

    return (void*)count;
}

int main(int argc, char* argv[])
{
    (void)argc;

    GOArray condgo = createGOAarray(10000);
    GOArray bggo = createGOAarray(1000000);

    printf("loading ...\n");

    START_CC(loading);
    if (!readGOFile(argv[1], &condgo) || !readGOFile(argv[2], &bggo)) {
        destroyGOArray(&condgo);
        destroyGOArray(&bggo);
        return -1;
    }
    STOP_CC(loading);
    PRINT_CC(loading);

    size_t range[] = { 0, 250, 500 };
    pthread_t thr[2] = { 0 };
    ProcessST pst[2];

    START_CC(compare);

    for (size_t i = 0; i < 2; i++) {
        pst[i].pcond = &condgo;
        pst[i].pbg = &bggo;
        pst[i].beg = range[i];
        pst[i].end = range[i + 1];
        pthread_create(&thr[i], NULL, processThread, &pst[i]);
    }

    int count = 0;
    int* ret = NULL;

    for (size_t i = 0; i < 2; i++) {
        pthread_join(thr[i], (void**)&ret);
        count += (int)ret;
    }

    STOP_CC(compare);
    PRINT_CC(compare);

    printf("result = %d\n", count);

    destroyGOArray(&condgo);
    destroyGOArray(&bggo);

    return 0;
}

Si vous le changez en multi-thread, la vitesse sera plus rapide. Sur ma machine, je l'ai changé en 2 threads et j'ai simplement utilisé 500条candidates pour tester. La vitesse peut être augmentée à 9040257微秒. le nombre de threads est augmenté à 4, les performances s'amélioreront. Ce n'est pas très gros, mais les processeurs les plus récents ont la technologie hyper-threading et la vitesse devrait être meilleure. . .

Compiler et tester

Désolé, j'ai vu quelqu'un d'autre répondre aujourd'hui.
En regardant attentivement la question, j'ai réalisé que je pensais que c'était juste un match.
J'ai donc proposé d'utiliser une machine automatique à courant alternatif.

Mais le but de la question est de trouver les différences dans les séquences.
Il s'agit de trouver la distance d'édition entre les deux.
wiki : Modifier la distance
wiki : Distance de Ravenstein

Quand je brossais OJ, j'utilisais DP (programmation dynamique) pour trouver le nombre minimum de modifications pour convertir une chaîne en une autre chaîne.

for(i:0->len)
    for(j:0->len)
        str1[i]==str2[j] ? cost=0 : cost=1
        dp[i,j]=min(
            dp[i-1, j  ] + 1,     // 刪除
            dp[i  , j-1] + 1,     // 插入
            dp[i-1, j-1] + cost   // 替換
        )

Par exemple :

str1    "abcdef"
str2    "acddff"

str2 est converti en str1

Insérer b Compter une fois
Supprimer d Compter une fois
Modifier f Compter une fois

Pour la séquence du gène ATCG du sujet, faut-il seulement retrouver celle modifiée ?
Cependant, comment doit-il être calculé ainsi
ATCGATCG
TCGATCGA
 ?

Si vous ne trouvez que des modifications, comparez simplement str1[i] et str2[i] directement.

for(i:0->len)
if(str1[i]!=str2[i] )++count;

Inspiré par @rockford.
Nous pouvons prétraiter les données brutes.

Chaîne en candidats
GGGAGCAGGCAAGGACTCTG
A5 T2 C4 G9

Données supplémentaires après traitement A5 T2 C4 G9

Chaîne en bg_db
CTGCTGACGGGTGACACCCA
A4 T3 C7 G6

Données supplémentaires après traitement A4 T3 C7 G6

A5 -> A4 est enregistré comme -1
T2 -> -3

Évidemment, A ne peut devenir TCG que s'il est modifié.

pour savoir au moins combien de différences ils ont.
Tout ce qui est supérieur à 4 n’a pas besoin d’être comparé.

En comparant d'abord les données supplémentaires prétraitées, puis en effectuant la comparaison via un algorithme de comparaison unique.

Vos tâches individuelles sont déterminées. Ce dont vous avez besoin est d'envoyer ces tâches aux travailleurs. De tels calculs ne sont pas effectués de manière synchrone en un seul processus.
En fait, cela équivaut à avoir [a, b] et [c, d] à comparer, votre tâche est

1. [a,c]

2. [a, d]

3. [b,c]

4. [b,d]

Si vous sérialisez de manière synchrone, le temps dont vous avez besoin est de 4 * temps unique
Si vous avez 4 processeurs ou 4 machines en parallèle, votre temps sera presque temps unique

Les calculs comme les génomes sont donc essentiellement effectués à l'aide de tâches parallèles multicœurs sur des machines à grande échelle. Fondamentalement, les principes de référence sont les principes de l'article Google MapReduce

Je ne suis pas doué en algorithmes, mais avec une grande quantité de données comme la vôtre, il n'est certainement pas possible de les comparer avec un seul ordinateur. Pour les tâches de données gourmandes en CPU comme la vôtre, je suis d'accord avec ce que d'autres ont dit sur l'utilisation. une méthode de cluster ou multi-processus pour calculer. Autrement dit, nous utilisons le modèle map-reduce pour calculer
map is mapping. Vous mappez d'abord chacun de vos candidats à bg_db un par un pour former une structure de données comme celle-ci (candidates,bg_db) dans une file d'attente puis transmettez-le à différents serveurs. , chaque serveur utilise plusieurs processus pour calculer, c'est le seul moyen, mais votre volume de données est trop important, trouvez un moyen d'attribuer vos tâches et calculez en parallèle. ,

Vous pouvez essayer d'utiliser une arborescence de dictionnaire pour enregistrer toutes les chaînes. Ensuite, vous pouvez utiliser la méthode de parcours de l'arborescence du dictionnaire lors de l'interrogation.
Lorsque vous parcourez l'arborescence, vous pouvez maintenir une séquence de nœuds actuels. Cette séquence stocke le nombre de discordances entre le nœud actuellement parcouru et le nœud correspondant.
Lorsque vous traversez le nœud suivant, essayez de descendre tous les nœuds de la séquence actuelle et de former une nouvelle séquence de nœuds.
L'avantage est que les bits actuels de plusieurs chaînes peuvent être comparés entre eux, ce qui peut faire gagner du temps. Puisqu’il n’y a pas beaucoup de choix pour chaque position et que le décalage n’est pas important, il ne devrait pas y avoir d’expansion excessive de la séquence de nœuds actuelle. (C'est une supposition... Je ne l'ai pas vérifié trop attentivement...)

def match(candidate, root):
  nset = [[],[]]
  currentId = 0
  nset[currentId].append((root, 0))
  for ch in candidate:
    nextId = 1 - currentId
    for item in nset[currentId]:
      node, mismatch = item
      for nx in node.child:
        if nx.key == ch:
          nset[nextId].append((nx, mismatch))
        else:
          if mismatch:
            nset[nextId].append((nx, mismatch - 1))
    nset[currentId].clear()
    currentId = 1 - currentId
  return nset[currentId]

Le code ci-dessus est une idée approximative. Ce serait beaucoup plus rapide s’il était écrit en C++.
L'ensemble du processus peut être effectué à l'aide de clusters pour l'informatique distribuée.

Visuellement, le questionneur n'a pas plusieurs machines à calculer.
J'ai une idée simple pour calculer la somme des nombres alphabétiques de chaque chaîne (A:0,T:1,C:2,G : 3), Calculez d'abord la différence entre la somme ordinale des deux chaînes. Le maximum ne peut pas dépasser 12. Quatre A deviennent quatre G. Si la différence est inférieure à 12, alors traitez-la
n'est qu'une idée approximative. le poids spécifique peut être Avec des réglages supplémentaires, cela peut théoriquement être beaucoup plus rapide.
Il existe un autre algorithme qui calcule la distance d'édition de chaîne (le nombre minimum de modifications pour ajouter, supprimer et modifier une chaîne par une autre chaîne. Je ne m'en souviens pas pour le moment. Vous pouvez le vérifier vous-même.

J'utilise blast blastn-short

:)