Qualification Round. Problem 3 Introns Detection

DNA원본과 splicing후 조각들의 read로 exon을 재조합하는 문제이다.

exon1 - intron1 - exon2 - intron2 - exon3 - intron3 이고

exon1 - exon2

exon2 - exon3

이런식으로 주어지면

exon1 - exon2 - exon3을 구하면 되는것이다.

 

간단하게 말하자면 원본 DNA seq에서 임의의 seq를 제거해 주어진 exon들을 모두 포함하는 새로운 sequence를 만들면 된다.

 

가장 단순하게 생각해보면 DNA seq에서 모든 base들을 하나씩 포함하거나 제외하면서 만든 seq가 exon들을 전부 포함하는지 보면 된다. base별로 포함/미포함 조합을 가지고 있으므로 길이 l인 sequence의 모든 조합만 구하는것만으로도 2^l의 복잡도를 가진다. exon의 갯수를 n, 길이를 k라고 하면 단순히 비교를 한다고 할때 n * k * (2^l)의 복잡도를 가지게 된다.

for spliced_dna in get_all_combi(dna_seq):
    for exon in exons:
        if exon not in spliced_dna:
            return False
    return True

 

다른 방법을 생각해보자. DNA seq는 백만의 길이가 넘어갈수도 있지만 exon의 갯수는 10번째 케이스에서도 4000개이다. 그렇다면 갯수가 적은 exon의 조합을 이용해보자. exon의 조합의 갯수는 n!이다. 각 조합별로 현재 DNA seq에서 splicing으로 얻을 수 있는지 여부를 검사하는덴 seq의 길이인 l의 복잡도가 필요하다. 그러므로 l * n!의 복잡도이다.

for spliced_dna in get_all_combi(exons):
    if can_be_spliced(dna_seq, spliced_dna):
        return True

위 두가지 방법으로는 도저히 현실적인 시간 안에 답을 구할수 없다. 여러 문자열을 substring으로 포함하는지 알아내는데는 suffix tree를 이용한 방법이 좋지만 문제는 splicing은 임의의 문자를 삭제하고 난 뒤에 substring 매칭을 해야한다는것이다.

 

조금만 생각해보면 원 DNA sequence는 4가지 염기로 이루어져있다. 그리고 seq가 랜덤하게되어있든 어떤 패턴을 가지고 있든 A가 없다거나 G로만 이루어져있는 극단적인 경우가 아니라면 충분한 길이가 있다면 어떤 subsequence든 만들어낼수있다. 그래서 다음과 같은 greedy한 방법으로 풀었다.

spliced_dna = ''
whlie True:
	min_last_idx_exon = min(map(lambda exon: process_and_get_last_idx(dna_seq, exon)))
    spliced_dna = weld_exon(spliced_dna, min_last_idx_exon)
    if reached_end():
    	break
return spliced_dna

위 pseudo코드를 실제로 구현한 코드는 아래와 같다.

import time

def remove_overlapped(exons):
    result = set()
    for e1 in exons:
        overlapped = False
        for e2 in exons:
            if e1 != e2 and e1 in e2:
                overlapped = True
        if overlapped is False:
            result.add(e1)
    return list(result)

def weld(exon1, exon2):
    max_common = 0
    for idx in range(len(exon2)):
        if exon1[-idx:] == exon2[0:idx]:
            max_common = idx
    return exon1[0:len(exon1)-max_common] + exon2

def process_dna(sequence, exon):
    seq_idx = -1
    for exon_base in exon:
        seq_idx += 1
        if seq_idx == len(sequence):
            return len(sequence) + 1
        while sequence[seq_idx] != exon_base:
            seq_idx += 1
            if seq_idx == len(sequence):
                return len(sequence)+1
    return seq_idx+1


if __name__ == '__main__':
    input_file = 'Q3'
    with open(input_file) as f:
        data = f.read().strip().split('\n')

    seq = data[0]
    exons = data[2:]

    exons = sorted(remove_overlapped(exons))

    cur_exon = ''
    cur_seq_end = 0
    old_time = time.time()
    while len(exons) != 0:
        end_idxs = []
        for exon in exons:
            weld_exon = weld(cur_exon, exon)
            new_exon_part = exon[len(cur_exon) + len(exon) - len(weld_exon):]
            end_idx = process_dna(seq[cur_seq_end:], new_exon_part) + cur_seq_end
            end_idxs.append(end_idx)
        min_exon_idx, seq_idx = min(enumerate(end_idxs), key= lambda x: x[1])
        if seq_idx <= len(seq):
            cur_exon = weld(cur_exon, exons.pop(min_exon_idx))
            cur_seq_end = seq_idx
            old_time = time.time()
        else:
            break
    with open('Q3_'+input_file, 'wt') as f:
        f.write(cur_exon + '\n')
        for exon in data[2:]:
            idx = cur_exon.find(exon)
            if idx != -1:
                idx += 1
            f.write(str(idx) + '\n')

완벽한 방법은 아니지만 꽤 적당한 성능을 내는 코드이다.

시간복잡도를 줄이기 위해 DNA seq를 exon과 비교해 splicing하는 부분을 새로 추가되는 exon부분만 splicing하도록 했다.

위 방법으로 seq길이를 l이라고 하고 exon갯수를 n, 길이를 s라고 할때 n * l * k의 복잡도를 가지게 구현할수있었다. (내 계산이 맞겠지?)

9번 케이스의 경우 1분정도 걸리고 200만의 seq길이와 2000길이의 4000개의 exon인 10번째 케이스는 1시간 반정도 걸려서 4000개중 3997의 exon을 매칭할 수 있었다.