Labels

Tuesday, February 10, 2015

Repeated DNA Sequences


Repeated DNA Sequences



 


All DNA is composed of a series of nucleotides abbreviated as A, C, G, and T, for example: "ACGAATTCCG". When studying DNA, it is sometimes useful to identify repeated sequences within the DNA.
Write a function to find all the 10-letter-long sequences (substrings) that occur more than once in a DNA molecule.
For example,
Given s = "AAAAACCCCCAAAAACCCCCCAAAAAGGGTTT",

Return:
["AAAAACCCCC", "CCCCCAAAAA"]. 
 
 
Naive Way: 一开始试了两种方法都不行。一种是用一个set去存
已经遍历过的substring,然后每取得一个新的就看set是否contain。 
这样的做法的算法复杂度是O(n),但是需要O(n)的space,
会出现内存不够。第二种方法是不用set,每次取得一个substring
都和之前开头到当前位置的每一个substring比较,
这样的算法复杂度是O(n^2),space O(1),会出现超时。
 
 


然后我想到了这个在time complexity 和 
memory之间折中的办法。Map存的是对应前m个字符一致时,
各个(m字符串)出现的位置,1<m<=10,这样m越大,space 
需求越大,time complexity越小。 但是还是不行,
在m=3的时候会报超时,m=4的时候会报空间不足。


public class Solution {
    public List<String> findRepeatedDnaSequences(String s) {
        List<String> rlst = new ArrayList<String>();
        Map<String, List<Integer>> map = new HashMap<String, List<Integer>>();
        int k = 10;
        int offset = 4;
        for(int i = 0;i < s.length()-k+1;i++){
            String ss = s.substring(i,i+offset);
            if(!map.containsKey(ss)){
                List<Integer> list = new ArrayList<Integer>();
                list.add(i);
                map.put(ss, list);
            }else{
                List<Integer> list = map.get(ss);
                boolean found = false;
                for(int j = 0;j < list.size();j++){
                    int t = offset;
                    for(;t < k;t++){
                        if(s.charAt(i+t)!=s.charAt(list.get(j)+t))
                            break;
                    }
                    if(t==k){
                        rlst.add(s.substring(i,i+k));
                        found = true;
                        break;
                    }
                }
                if(!found)
                    list.add(i);
            }
        }
        return rlst;
    }
}


然后还试了一种树形的存贮方式,也是memory limit exceed。看了tag标记为bit manipulate,似乎有点头绪了。四种字符可以用两位二进制完全表示。10个连续字符是20位二进制,一个int型是32位。是否可以用一个Set<Integer>来存遍历过的字符串,这样每一次存贮可以节省16*10-32=128bit 空间。这种方法还是会内存不足。










 最后的最后!



  Improved Way:



终于让这道题Accepted了,已经是第五天了,这回我也算是竭尽所能的减少空间使用率。首先使用了一个byte array来存10个连续单位中A,C,G,T出现次数 ,又因为正好4个byte是一个int,就可以用一个int来表示出现次数,使用一个Map<Integer, Set<String>>的map来将出现次数和同一种出现次数下的对应String Set做映射。然后由于10个单位的两比特数是20比特,又可以用一个int表示,就将所有String换成int来存,这样最后的map就变成了Map<Intger, Set<Integer>>。每次只要看该String对应的Integer是否在Set中就可以判断是否出现过了。







public class Solution {
    public List<String> findRepeatedDnaSequences(String s) {
        List<String> list = new ArrayList<String>();
        Map<Integer, Set<Integer>> map = new HashMap<Integer, Set<Integer>>();
        Set<String> visited = new HashSet<String>();
        int k = 10;
        byte[] appear = new byte[4];
        Arrays.fill(appear, (byte)0);
        for(int i = 0;i < s.length();i++){
            appear[dna2int(s.charAt(i))]++;
            if(i >= k-1){
                String sample = s.substring(i-k+1,i+1);
                int key = byte2int(appear);
                int intSample = seq2int(sample);
                if(!map.containsKey(key)){
                    Set<Integer> set = new HashSet<Integer>();
                    set.add(intSample);
                    map.put(key,set);
                }else{
                    Set<Integer> set = map.get(key);
                    if(set.contains(intSample)){
                        if(!visited.contains(sample))
                            visited.add(sample);
                    }else{
                        set.add(intSample);
                        //map.put(key,set);
                    }
                }
                appear[dna2int(s.charAt(i-k+1))]--;
            }
        }
        list.addAll(visited);
        return list;
    }
    
    private int dna2int(char a){
        switch(a){
            case 'A':return 0;
            case 'C':return 1;
            case 'G':return 2;
            case 'T':return 3;
            default: return 0;
        }
    }
    
    private int byte2int(byte[] b){
        int x = b[0] << 24 | b[1] << 16 | b[2] << 8 | b[3];
        return x;
    }
    
    private int seq2int(String s){
        int x = 0;
        for(int i = 0;i < s.length();i++){
            x <<= 2;
            switch(s.charAt(i)){
                case 'A':
                    x |= 0;
                    break;
                case 'C':
                    x |= 1;
                    break;
                case 'G':
                    x |= 2;
                    break;
                case 'T':
                    x |= 3;
                    break;
                default:
                    break;
            }
        }
        return x;
    }
    
}



 



 



最后在discuss里看到了一个最厉害的方法,叫做Rabin-Karp,可用来解strStr()。它采用Rolling Hash的方法,每次固定大小的窗移动一位时,只需减去对应的hash量和加上新的hash量就可以形成新的hash码,它的理论基础大概是这样的。



 



对于这道题,一共有ACGT四个变量,为他们赋值0123,那么以4为底的话,就能让连续k个组合的序列具有唯一的hash码。比如k=3时,




ACG = 0 * 4^2 + 1 * 4^1 + 2 * 4^0;



TTC  = 3 * 4^2 + 3 * 4^1 + 1 * 4^0;



就像10进制,只要每一位上的数不是完全相同,两个数就不相同。



使用这个的前提是4^(k+1) <= Integer.MAX_VALUE。



k=10的时候 4^11 = 2^22 < 2^31。






维基百科中说,即使变量的个数很大时,也可以用小于变量个数的大质数做底,因为质数幂次方后很少出现两个不同数的hash码一致。









public class Solution {
    private static final int base = 4;
    public List<String> findRepeatedDnaSequences(String s) {
        // Rabin-Karp rolling hash
        List<String> list = new ArrayList<String>();
        Set<Integer> set = new HashSet<Integer>();
        Set<String> rslt = new HashSet<String>();
        int k = 10;
        int head_factor = (int)Math.pow(base,k-1);
        int hash = 0;
        for(int i = 0;i < s.length();i++){
            hash -= i>=k?head_factor*dna2int(s.charAt(i-k)):0;
            hash = hash*base+dna2int(s.charAt(i));
            if(i >= k-1 && !set.add(hash))
                rslt.add(s.substring(i-k+1,i+1));
        }
        list.addAll(rslt);
        return list;
    }
    
    private int dna2int(Character c){
        switch(c){
            case 'A':return 0;
            case 'C':return 1;
            case 'G':return 2;
            case 'T':return 3;
            default:return 0;
        }
    }
}





No comments:

Post a Comment