문제 해설

Next Generation Sequencing(NGS가) 개발되고 나서, 저렴한 비용으로 DNA의 sequencing을 할 수 있게 되었습니다. NGS는 전체 DNA를 읽지는 못하고 DNA를 작게 쪼개 염기서열을 읽어들입니다. 그렇게 읽은 작은 DNA조각을 read라고 합니다.

그런데 NGS는 만능이 아니기에 주의를 기울여야 합니다. NGS는 무시 못할 확률로 리드를 읽을 때 에러를 발생시키기 때문입니다. 심지어 에러는 어떤 파트에 발생할지 예측할 수 없습니다. 그러므로 Genome Assembly를 수행 할 때는 반드시 에러를 수정하는 절차를 거쳐야 합니다.

문제 해석

원문

As is the case with point mutations, the most common type of sequencing error occurs when a single nucleotide from a read is interpreted incorrectly.

Given : collection of up to 1000 reads of equal length (at most 50 bp) in FASTA format. Some of these reads were generated with a single-nucleotide error. For each read $s$ in the dataset, one of the following applies:

• $s$ was correctly sequenced and appears in the dataset at least twice (possibly as a reverse complement);
• $s$ is incorrect, it appears in the dataset exactly once, and its Hamming distance is 1 with respect to exactly one correct read in the dataset (or its reverse complement).

Return : A list of all corrections in the form "[old read]->[new read]". (Each correction must be a single symbol substitution, and you may return the corrections in any order.)

해석

Point Mutation처럼 오직 read에서 한개의 nucleotide에서만 에러가 발생한다고 합니다. 주어지는 건 최대 1000개의 리드가 주어집니다. 이 중에서 리드는 맞을 수도 있고 틀릴 수도 있습니다.

• 맞는 리드의 경우 자신과 같은 리드(reverse compliment을 포함하여)가 2개 이상 있을 수 있습니다.
• 틀린 리드의 경우 오직 하나의 옳바른 리드에 대해 hamming distance가 1만큼 떨어져 있습니다.

그리고 Return으로 [old read]->[new read]형식으로 출력하라고 합니다.

(여기서 reverse compliment란 $TCATC = reverseCompliment(GATGA)$ 를 의미합니다. )
문제가 이상하게 여겨질 수 있습니다. 동일한 리드가 2개 이상 있을 수 있는데 어떻게 단 하나의 옳은 리드와만 매칭이 될 수 있는지 의문이 생기실지도 모르겠습니다.

예를 들어 $TCGTC$가 $TCATC$와 error -> correct pair를 맺을 수 있다고 합시다. 이때, $TCATC$가 2개가 있던 3개가 있던 상관 없습니다. 왜냐하면 correct는 모두 동일함으로 error->correct pair는 한개입니다.

$TCGTC$ 와 $reverseCompliment(x)$가 있어 페어를 이룰 수 있다고 해도, 결국 $x$를 $reverseCompliemnt$한 값과 pair를 맺어야 하기 때문에 결국 error->correct 페어는 한개가 됩니다.

만일 리드에서 서로 다른 위치에 에러가 있다면, 여러개와 페어가 맺어 질 가능성이 있지만, 문제 자체에서 그런 가능성을 봉쇄해 놓았습니다. 그러므로 우리는 좀더 편하게 문제에 접근 할 수 있습니다.

문제 풀이

문제의 데이터가 지나가는 파이프 라인은 다음과 같습니다.

find_correct: seqs -> correct_set, seqs

find_error : correct_set, seqs -> correct_set, error_set

correction : correct_set, error_set -> pairs

run : seqs -> correction(find_error(find_correct(x)))

파이프라인의 순서대로 프로그램은 3가지 순서대로 진행이 됩니다.

1. 옳바른 리드를 찾기
2. 에러 찾기
3. 페어 맞추기

---

1. 옳바른 리드 찾기
   옯바른 리드의 조건은 자신과 같은 리드를 찾는 것입니다. read들을 2중으로 순회하며, 자기 자신을 제외한 다른 리드에서 자신과 같은(reverse compliment 포함) 리드가 있을 경우 저장해 둡니다. 그리고 원본 read set과 옳바른 read들의 set을 반환합니다.

def find_correct(seqs):
    correct = []
    for i, seq_i in enumerate(seqs):
        for seq_j in seqs[:i] + seqs[i+1: ]:
            if seq_i == seq_j or seq_i == rev_seq(seq_j):
                correct.extend([seq_i, seq_j])
    return set(correct), seqs

2. 에러 찾기
   옳바른 read의 set을 얻었으니, 전체 리드에서 옳바른 리드를 제외한 나머지 리드가 error가 됩니다. $error\_set = set(reads) / correct\_set$ 으로 error_set을 얻어 낼 수가 있습니다.

def find_error(datapacakge):
    correct, seqs = datapacakge
    errors = set(seqs).difference(correct)
    return correct, errors

3. 페어 맞추기
   마지막으로 pair를 맞춰주면 됩니다. pair의 조건은 hamming distance가 1만큼 떨어진 리드를 옳바른 리드로부터 찾아내는 것입니다. 이때, hamming distance 1의 조건은 straight와 reverse compliment 케이스 두개를 모두 조사해야 합니다. 그러면 문제에서 알려 준 것과 같이 딱 한개의 데이터만 correct로 나오게 됩니다. 이제 해야 할 일은 옳은 read 셋에서 straight가 hamming distance가 1이면 straight 대로, reverse compiment와 hamming distance가 1이면 reverse compliment로 변환한 read를 error read와 매칭해주기만 하면 됩니다.

def correction(datapacakge):
    correct, errors = datapacakge
    pairs = []
    for error in errors:
        cor = first([cor for cor in correct if hamming_distance(cor, error) == 1 or hamming_distance(rev_seq(cor), error) == 1])
        if hamming_distance(error, cor) == 1:
            pairs.append((error, cor))
        elif hamming_distance(error, rev_seq(cor)) == 1:
            pairs.append((error, rev_seq(cor)))
    return pairs

전체 코드는 github에 올려놨습니다.