专栏首页生信了序列比对(26)精准匹配之KMP算法、Trie树以及AC自动机

序列比对(26)精准匹配之KMP算法、Trie树以及AC自动机

前文已经介绍过KMP算法和Trie树,本文将在此基础上介绍AC自动机。

之前的序列比对文章大都在利用动态规划算法解决字符串的非精准匹配(允许错配、插入和缺失),比如全局比对和局部比对问题。当然,后来我们还介绍了模序发现和中间字符串问题,并初次学习了如何用分支定界法解决这一类问题。

既然有非精准匹配问题,就有精准匹配问题。所谓精准匹配,就是两个字符串在比对时,不允许错配、插入和缺失。其实,笔者之前已经有过介绍:

KMP算法

算法(四)(转载)KMP算法》一文介绍了KMP算法,该算法常用来解决如何在一个长字符串(模板串)中寻找一个短字符串(模式串)的问题。其特点就是要预先对模式串进行处理,对其构建一个next数组,该数组中存放着模式串每一个位置所对应的跳转指针,这样做的效果是模板串的指针不用回溯,一直在递增,从而增加效率。

KMP算法的完整代码:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define MAXLEN 1000

void getNext(const char *pat, int *next, const int n) {
/* assume n >= 2 */
    int i, j;
    for (next[0] = -1, next[1] = 0, i = 1, j = 0; i < n - 1;) {
        while (j != -1 && pat[i] != pat[j])
            j = next[j];
        next[++i] = ++j;
    }
}

int KMP(const char *temp, const int m, const char *pat, const int n) {
    int i, j;
    int *next = (int*) malloc(sizeof(int) * n);
    getNext(pat, next, n);
    for (i = 0, j = 0; i < m && j < n;) {
        if (j == -1 || temp[i] == pat[j]) {
            i++;
            j++;
        } else {
            j = next[j];
        }
    }
    free(next);
    if (j < n)
        return -1;
    return i - j;
}

int main(void) {
    char t[MAXLEN], p[MAXLEN];
    int index;
    printf("the temp str: ");
    scanf("%s", t);
    printf("the pattern str: ");
    scanf("%s", p);
    if ((index = KMP(t, strlen(t), p, strlen(p))) == -1)
        printf("Not found!\n");
    else
        printf("Index is %d.\n", index);
    return 0;
} 

其效果如下:

Trie树(字典树)

算法(五)字典树算法快速查找单词前缀》一文介绍了Trie树,该算法常用来解决如何在多个模板串中寻找一个模式串的问题。其特点就是将多个模板串装填进一个树结构中,虽然使用了较多的空间,但是查找效率得到了提升。

Trie树的完整代码:

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#define N 26
#define MAXLEN 1000

struct Node {
    struct Node *a[N];
    int is_end;
};
typedef struct Node *tnode;

tnode init() {
    tnode p;
    int i;
    if ((p = (tnode) malloc(sizeof(struct Node))) == NULL) {
        fputs("Error: out of space!\n", stderr);
        exit(1);
    }
    for (i = 0; i < N; i++)
        p->a[i] = NULL;
    p->is_end = 0;
    return p;
}

void delete(tnode p) {
    int i;
    for (i = 0; i < N; i++)
        if (p->a[i] != NULL)
            delete(p->a[i]);
    free(p);
}

void insert(const char *s, tnode p) {
    int i;
    while (*s != '\n' && *s != '\0') {
        i = *s - 'a';
        if (p->a[i] == NULL) {
            p->a[i] = init();
        }
        p = p->a[i];
        s++;
    }
    p->is_end = 1;
}

int find(const char *s, tnode p) {
    int i;
    while (*s != '\n' && *s != '\0') {
        i = *s - 'a';
        if (p->a[i] == NULL)
            return 0;
        p = p->a[i];
        s++;
    }
    if (p->is_end)
        return 1;
    return 0;
}

void strLower(char *s) {
    while (*s != '\0') {
        if (*s >= 'A' && *s <= 'Z') {
            *s += 32;
        }
        s++;
    }
}

int main(int argc, char *argv[]) {
    FILE *fpp, *fpw;
    char buf[MAXLEN];
    tnode root = init();
    if (argc < 3) {
        fputs("Usage: trie.exe <pattern.txt> <words.txt>\n", stderr);
        exit(2);
    }
    if ((fpw = fopen(argv[2], "r")) == NULL) {
        fprintf(stderr, "Error: cannot open %s\n", argv[2]);
        exit(4);
    }
    while (fgets(buf, MAXLEN, fpw) != NULL) {
        strLower(buf);
        insert(buf, root);
    }
    if (! feof(fpw)) {
        fprintf(stderr, "Error: error while reading %s\n", argv[2]);
        exit(8);
    }
    fclose(fpw);
    if ((fpp = fopen(argv[1], "r")) == NULL) {
        fprintf(stderr, "Error: cannot open %s\n", argv[1]);
        exit(4);
    }
    while (fgets(buf, MAXLEN, fpp) != NULL) {
        strLower(buf);
        if (find(buf, root)) 
            printf("%s", buf);
    }
    if (! feof(fpp)) {
        fprintf(stderr, "Error: error while reading %s\n", argv[1]);
        exit(8);
    }    
    fclose(fpp);
    delete(root);
}

其效果如下:

其中,pat.txt文件内容(里面是多个模式串)如下:

how are you

words.txt文件内容(里面是多个模板串)如下:

lily we are friends

AC自动机

本文要介绍的AC自动机算法,是用来解决如何有效地在一个模板串中查找多个模式串的问题。该算法是建立在KMP算法和Trie树基础上的。其特点是:

  • 预先将多个模式串装填进Trie树中。
  • 在进行查找(匹配)之前,对Trie树中的每个有效节点构建fail指针,该指针的作用类似KMP算法中next数组的作用:如果到该节点时匹配失败,就可以利用fail指针跳转到下一个可利用节点继续进行匹配,从而避免了模板串的回溯而提升查找效率。

如同KMP算法的next数组充分利用了模式串的内部信息,AC自动机中的fail指针也是充分利用了多个模式串的内部信息,每一次跳转都是跳到“最大可利用后缀子字符串”的节点。

笔者在学习该算法时,主要参考了这一篇文章(https://www.jianshu.com/p/641142c6d477),看里面的图非常有助于理解AC自动机原理。

关于AC自动机的实现有几点说明:

  • 在对Trie树里所有有效节点构建fail指针的过程中,采用了BFS(广度优先搜索)的策略,具体的实现是利用了队列这种数据结构。采用BFS的原因是Trie树中某一层节点的fail指针指向的节点一定是来自于上面的层(可能是上面第一层,也有可能是上面第二层等等,甚至是根节点)。
  • 无论是Trie树还是队列,都是利用一种链表实现的,不得不感叹,对链表进行操作真是很容易出错,经常要填坑。

AC自动机的完整代码如下:

#include <stdio.h>
#include <stdlib.h>
#define ALPHA_NUM 26
#define QUEUE_EMPTY NULL
#define BUFSIZE 64
#define MAXTEMP 10000

/* Data structure for trie */
struct TNode {
    struct TNode *a[ALPHA_NUM];
    struct TNode *fail;
    struct TNode *parent;
    int c;
    int is_end;
};
typedef struct TNode *tnode;

tnode initTrie(int c, tnode parent) {
    tnode p;
    int i;
    if ((p = (tnode) malloc(sizeof(struct TNode))) == NULL) {
        fputs("Error: out of space!\n", stderr);
        exit(1);
    }
    for (i = 0; i < ALPHA_NUM; i++)
        p->a[i] = NULL;
    p->fail = NULL;
    p->parent = parent;
    p->c = c;
    p->is_end = 0;
    return p;
}

void deleteTrie(tnode p) {
    int i;
    for (i = 0; i < ALPHA_NUM; i++)
        if (p->a[i] != NULL)
            deleteTrie(p->a[i]);
    free(p);
}

void insertTrie(const char *s, tnode p) {
    int i;
    while (*s != '\n' && *s != '\0') {
        i = *s - 'a';
        if (p->a[i] == NULL) {
            p->a[i] = initTrie(*s, p);
        }
        p = p->a[i];
        s++;
    }
    p->is_end = 1;
}

/*
void printTrie(tnode p) {
    int i;
    printf("%c ", p->c);
    for (i = 0; i < ALPHA_NUM; i++) {
        if (p->a[i] != NULL)
            printTrie(p->a[i]);
    }
}
*/

void backtrack(char *s, const int n, tnode p) {
    int i, j, c;
    for (i = 0; i < n - 1 && p != NULL && p->parent != NULL; i++) {
        s[i] = p->c;
        p = p->parent;
    }
    if (i >= n - 1) {
        fputs("Error: pattern word is too long!\n", stderr);
        exit(2);
    }
    for (j = 0; j < i / 2; j++) {
        c = s[j];
        s[j] = s[i - 1 - j];
        s[i - 1 - j] = c;
    }
    s[i] = '\0';
}

/* Data structure for queue */
typedef tnode QNodeType;

struct QNode {
    QNodeType value;
    struct QNode *next;
};
typedef struct QNode *qnode;

qnode initQueue(QNodeType v) {
    qnode p;
    if ((p = (qnode) malloc(sizeof(struct QNode))) == NULL) {
        fputs("Error: out of space!\n", stderr);
        exit(1);
    }
    p->value = v;
    p->next = NULL;
    return p;
}

void addQueue(QNodeType v, qnode head) {
    qnode p;
    for (p = head; p->next != NULL; p = p->next)
        ;
    p->next = initQueue(v);
}

int isEmptyQueue(qnode head) {
    return head->next == NULL ? 1 : 0;
}

QNodeType popQueue(qnode head) {
    qnode p;
    QNodeType q;
    if (isEmptyQueue(head))
        return QUEUE_EMPTY;
    else {
        p = head->next;
        head->next = p->next;
        q = p->value;
        free(p);
        return q;
    }
}

void destroyQueue(qnode head) {
    while (popQueue(head) != QUEUE_EMPTY)
        ;
    free(head);
}

void strLower(char *s) {
    while (*s != '\0') {
        if (*s >= 'A' && *s <= 'Z') {
            *s += 32;
        }
        s++;
    }
}

void getFailPointer(tnode root) {
    int i;
    qnode head = initQueue(NULL);
    tnode this, tmp;
    addQueue(root, head);
    while (! isEmptyQueue(head)) {
        this = popQueue(head);
        for (i = 0; i < ALPHA_NUM; i++) {
            if (this->a[i] != NULL) {
                addQueue(this->a[i], head);
                tmp = this;
                while (tmp->fail != NULL) {
                    if (tmp->fail->a[i] != NULL) {
                        this->a[i]->fail = tmp->fail->a[i];
                        break;
                    } else {
                        tmp = tmp->fail;
                    }
                }
                if (tmp->fail == NULL)
                    this->a[i]->fail = root;
            }
        }       
    }
    destroyQueue(head);
}

void ACauto(char *buf, const int n, char *temp, tnode root) {
    char *s = temp;
    tnode p, tmp;
    int i;
    for (p = root; *s != '\0'; s++) {
        i = *s - 'a';
        while (p->a[i] == NULL && p->fail != NULL)  /* 这里必须先判断p->a[i],否则有可能错过root->a[i] != NULL的情况 */
            p = p->fail;
        if ((p = p->a[i]) == NULL)
            p = root;
        for (tmp = p; tmp->fail != NULL; tmp = tmp->fail)
            if (tmp->is_end) {
                backtrack(buf, BUFSIZE, tmp);
                printf("%s\n", buf);
            }
    }
}

int main(int argc, char *argv[]) {
    FILE *fpp, *fpt;
    char buf[BUFSIZE];
    char temp[MAXTEMP];
    tnode root = initTrie('\0', NULL);
    int n;
    if (argc < 3) {
        fputs("Usage: trie.exe <pattern.txt> <template.txt>\n", stderr);
        exit(2);
    }
    if ((fpp = fopen(argv[1], "r")) == NULL) {
        fprintf(stderr, "Error: cannot open %s\n", argv[1]);
        exit(4);
    }
    while (fgets(buf, BUFSIZE, fpp) != NULL) {
        strLower(buf);
        insertTrie(buf, root);
    }
    if (! feof(fpp)) {
        fprintf(stderr, "Error: error while reading %s\n", argv[1]);
        exit(8);
    }
    fclose(fpp);
    if ((fpt = fopen(argv[2], "r")) == NULL) {
        fprintf(stderr, "Error: cannot open %s\n", argv[2]);
        exit(4);
    }
    n = fread(temp, 1, MAXTEMP - 1, fpt);
    temp[n] = '\0';
    strLower(temp);
    if (! feof(fpt)) {
        fprintf(stderr, "Error: error while reading %s\n", argv[2]);
        exit(8);
    }    
    fclose(fpt);
    getFailPointer(root);
    ACauto(buf, BUFSIZE, temp, root);
    deleteTrie(root);
}

其效果如下:

其中pat.txt文件内容(里面是多个模式串)如下:

how are you

temp.txt文件内容(里面是模板串)如下:

iloveautumnhowaboutyou

至此,本文就全部结束了,谢谢大家看此长文。如果有发现错误或者有任何建议和意见,欢迎指正!

本文分享自微信公众号 - 生信了(gh_ed36a29a9a9d),作者:hxj7

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-09-26

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 算法(四)(转载)KMP算法

    字符串匹配是一个既古老又现代的问题,历久弥新。生信领域中字符串处理更是daily work。诸如bwa这般神一样的软件,本质上也是在解决字符串非精准匹配的问题。...

    一只羊
  • 序列比对(14)viterbi算法和后验解码的比较

    前文《序列比对(十)viterbi算法求解最可能路径》介绍了用viterbi算法求解最可能路径:在符号序列已知而状态序列未知时,最可能路径是:

    一只羊
  • 算法(二)蓄水池抽样算法快速随机抽取reads

    fastq文件往往都很大,出于测试目的,我们经常要从fastq文件中随机抽取reads,生成一个小一点的fastq文件,以加快测试效率。假设我们要从一个包含大约...

    一只羊
  • python基础教程:包,对,没错,绝对不是双肩包!

    包,Package,是一种Python模块的集合,从文件组织形式上看,包就是一个文件夹,里面放着各种模块(.py文件),也可以有子文件夹(子包)。包名构建了一个...

    一墨编程学习
  • RedHat 开源企业镜像项目 Quay

    Quay 是一个registry,存储,构建和部署容器的镜像仓库。它分析您镜像中的安全漏洞,可帮助您减轻潜在的安全风险问题。此外,它提供地理复制和BitTorr...

    YP小站
  • 朝阳群众,你关注的问题,答案都在这里!

    近期,后台很多粉丝留言想要咨询举报相关的问题,为此小助手特意整理了一份官方版举报攻略,希望大家大伙儿和小助手共同携手打造一个干净的网络世界。

    用户6966869
  • 独爱 Vim 的Linux老司机理由竟然是这个!!

    Vim是一个类似于Vi的著名的功能强大、高度可定制的文本编辑器,在Vi的基础上改进和增加了很多特性。VIM是自由软件。 Vim普遍被推崇为类Vi编辑器中最好的一...

    小小科
  • Spring support in Geronimo,看来It行业新的就是好的

    There was an EJB 3 talk at the "licensee" day, which happens before JavaOne star...

    田春峰-JCJC错别字检测
  • Blazor带我重玩前端(三)

    需要升级VS2019以及.NET Core到最新版(具体的最低支持,我已经忘了,总是越新支持的就越好),以更好的支持自己开发Blazor项目。

    Edison.Ma
  • poj 2251 Dungeon Master(广搜)

    题意:三维空间,可以走上下左右前后六个方向,求最短路径,BFS #include<stdio.h> #include<queue> #include<strin...

    用户1624346

扫码关注云+社区

领取腾讯云代金券