検索

ホームページ  >  に質問  >  本文

如何提高大量的字符串比较速度

补充数据链接:

candidates: https://pan.baidu.com/s/1nvGWbrV
bg_db: https://pan.baidu.com/s/1sllFLAd

每个字符串长度都是23,只要前面20个字符就行,由于数据太大,我只传了五分之一,大神们可以挑战一下,有速度快的可以贴一下代码,让小弟拜读一下,谢谢!

下面是正式的问题:

我现在有两个字符串数组,姑且称为candidates和bg_db,全部都是长度为20的短字符串,并且每个字符串的每个字符只有ATCG四种可能(没错!就是基因组序列啦!):

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

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

我的任务是对candidates的每一个candidate,找到bg_db中所有与其小于等于4个差异的记录,举个例子来说:

# 上面一条为candidate,即candidates的一个记录
# 中间|代表相同,*代表不相同
# 下面一条代表bg_db的一条记录

A T C G A T C G A T C G A T C G A T C G
| | | | | | | | | | | | | | | | | | | |    # 差异为0
A T C G A T C G A T C G A T C G A T C G

A T C G A T C G A T C G A T C G A T C G
* | | | | | | | | | | | | | | | | | | |    # 差异为1
T T C G A T C G A T C G A T C G A T C G

A T C G A T C G A T C G A T C G A T C G
* | | | * | | | | | | | | | | | | | | |    # 差异为2
T T C G T T C G A T C G A T C G A T C G

A T C G A T C G A T C G A T C G A T C G
* | | | * | | | | | | * | | | | | | | |    # 差异为3
T T C G T T C G A T C C A T C G A T C G

A T C G A T C G A T C G A T C G A T C G
* | | | * | | | | | | * | | | * | | | |    # 差异为4
T T C G T T C G A T C C A T C A A T C G

我的问题是如果快速地找到:每一个candidate在bg_db中与之差异小于等于4的所有记录,如果采用暴力遍历的话,以Python为例:

def align(candidate, record_from_bg_db):
    mismatches = 0
    for i in range(20):
        if candidate[i] != record_from_bg_db[i]:
            mismatches += 1
            if mismatches >= 4:
                return False
    return True

candidate = 'GGGAGCAGGCAAGGACTCTG'
record_from_bg_db = 'CTGCTGACGGGTGACACCCA'

align(candidate, record_from_bg_db) # 1.24微秒左右

# 总时间:

10000000 * 1000000000 * 1.24 / 1000 / 1000 / 60 / 60 / 24 / 365
# = 393
# 1千万个candidates,10亿条bg_db记录
# 耗时大约393年
# 完全无法忍受啊

我的想法是,bg_db是高度有序的字符串(长度固定,每个字符的可能只有四种),有没有什么算法,可以让candidate快速比较完所有的bg_db,各位大神,求赐教。

高洛峰高洛峰2941日前873

全員に返信(1)返信します

  • 三叔

    三叔2016-10-27 09:29:44

    写一个思路

    candidates = [
        'GGGAGCAGGCAAGGACTCTG',
        'GCTCGGGCTTGTCCACAGGA',
        '...',
        # 被你看出来啦,这些其实人类基因的片段
    ]
    
    bg_db = [
        'CTGCTGACGGGTGACACCCA',
        'AGGAACTGGTGCTTGATGGC',
        '...',
        # 这个更多,有十亿左右
    ]

    因为你的数据其实是很有特点的,这里可以进行精简。
    因为所有的字符串都是20个字符长度,且都由ATCG四个字符组成。那么可以把它们变换为整数来进行比较。
    二进制表现形式如下

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

    因为一个字符串长度固定,每个字符可以由2个比特位表示,所以每个字符串可以表示为一个40位的整数。可以表示为32+8的形式,也可以直接使用64位整形。建议使用C语言来做。

    再来说说比较。
    因为要找到每一个candidate在bg_db中与之差异小于等于4的所有记录,所以只要两个整数一做^按位异或操作,结果中二进制中1不超过8个,且这不超过8个1最多只能分为4个组的才有可能是符合要求的(00^11=11,10^01=11)。
    把结果的40个比特位分作20个组,那么就是说最多只有4个组为b01 b10 b11这三个值,其余的全部为b00。
    那么比较算法就很好写了。
    可以对每个字节(四个组)获取其中有几个组是为三个非零值的,来简介获取整体的比较结果。
    因为每个字节只有256种可能的值,而符合条件的值只有3^4=81种,所以可以先将结果存储起来,然后进行获取。
    这里给出一个函数,来获取结果中有几个是非零组。

    /*****************下面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]];
    }


    返事
    0
  • キャンセル返事