[REVP] Complementing a Strand of DNA
안녕하세요. 일을 쉬면서 감을 다시 찾기 위해 Rosalind를 다시 사용하고 있습니다.
그동안 제가 사용한 코드를 다시 보는 것도 좋지만 아무래도 굳어진 뇌를 다시 활성화하고 python이라는 언어의 감을 되찾으려면 백준이나 Rosalind를 이용하는 것이 더 괜찮다고 생각해 몇 문제를 풀어보았는데요.
아무래도 오랜만에 풀다 보니 스크립트가 참고용으로 보기에도 난잡해보일 수 있어 양해 부탁드립니다.
최근 업무가 epigenome 관련 업무였다 보니 이런 여러 restriction site가 생각나 매우 익숙하게 느껴졌습니다.
Problem
A DNA string is a reverse palindrome if it is equal to its reverse complement. For instance, GCATGC is a reverse palindrome because its reverse complement is GCATGC. See Figure 2.
Given: A DNA string of length at most 1 kbp in FASTA format.
Return: The position and length of every reverse palindrome in the string having length between 4 and 12. You may return these pairs in any order.
Sample Dataset
>Rosalind_24
TCAATGCATGCGGGTCTATATGCAT
Sample Output
4 6
5 4
6 6
7 4
17 4
18 4
20 6
21 4
이 문제를 풀기 위한 스크립트를 짜는데 있어서 고려해야할 점은 다음과 같다.
- 기존 염기서열과 reverse complement된 염기서열을 비교하기 위해 기존 염기서열을 바꿔준 뒤 비교하는 작업
- 결과값을 도출할 때 염기서열의 위치가 첫 번째 기준으로 정렬되고 reverse complement로 판단된 길이가 두 번째로 정렬되도록 하는 것
- 마지막 부분을 비교할 때 염기서열의 수보다 많은 부분을 읽으려 하는 것을 방지하는 것
(마지막 남은 (n - 12)bp의 염기서열에서 비교작업 진행 후 (n - 11)bp에 대하여 다시 4, 6, 8, 10bp의 reverse complement 비교 시 12bp에 대한 비교는 예외처리를 하는 것.) - 가능하다면 방대한 양을 처리한다고 했을 때 메모리 및 CPU 사용량을 최소화하는 것
Solution
물론 >Rosalind 라는 header 부분을 제외하고 진행하는 과정 역시 필요하지만 이 부분은 if문과 startswith를 사용한 방법이나 읽을 때 1번째 줄부터 읽는 방법 들이 있으며 fasta 파일의 특징인 줄바꿈은 strip()을 사용해 '\n'을 제거할 수 있다.
염기서열의 경우 "" 자리에 open을 통해 문제의 파일을 읽은 뒤 위 방법 중 하나를 사용하면 된다.
현 문제에서 reverse complement의 길이를 4, 6, 8, 10, 12bp로 제한했기 때문에 이에 대한 변수로 rc_len을 설정했다.
dna = ""
"""
f = open("/path/rosalind_revp.txt")
s = "".join(f.read().strip().split("\n")[1:]).strip()
"""
rc_len = [4,6,8,10,12]
가장 먼저 Reverse complement에 대한 비교를 위해 염기서열을 바꿔주는 작업이 필요하다.
해당 작업은 계속해서 사용될 것이기 때문에 def문을 사용하여 처리했다.
이는 기존의 문제 중 하나였던 [REVC] Complementing a Strand of DNA 에서의 작업과 동일하게 진행하였다.
(이 과정을 dictionary를 사용하여 진행한 사람도 확인했다.)
# get rev_compl
def cal_rev(bp):
bp_rev = bp.replace('A', 't').replace('T', 'a').replace('G', 'c').replace("C", 'g').upper()[::-1]
return bp_rev
두 번째로 기존 서열과 reverse complement 서열을 비교하는 작업이 필요하다.
이 때 결과값을 도출 했을 때 bp(염기서열의 위치)를 기준으로 먼저 정렬될 수 있도록 작성해야 한다.
# calculate rc
def cal_rc(dna, rc_len, pos):
fst_dna = dna[:int(rc_len/2)] #ATGC
sec_dna = dna[int(rc_len/2):rc_len] #GCAT
rev_dna = cal_rev(fst_dna) #GCAT
if sec_dna == rev_dna:
print(pos, rc_len)
이후 위 2개의 def문을 사용하여 for문 내에서 사용하여 진행한다.
첫 번째 for문에서는 전체 서열 중 마지막 12개의 서열 전까지만 비교를 진행한다.
두 번째 for문에서 마지막 12개의 서열에 대하여 진행한다.
이렇게 작성한 이유는 단순하게 마지막 12개의 서열에 대해서 동일하게 4~12개의 bp에 대한 검사를 할 시 오류가 발생할 수 있으며 이를 위해 if문을 설정할 경우 코드가 지저분해지기도 하지만 이를 위해 if문의 계산이 계속 진행되기 때문에 불필요한 계산을 줄이기 위해 넣었다.
for문 내에서 진행 과정은 slide window 방식을 참고하여 변수의 저장 및 사용에 대한 메모리를 줄이려고 노력했다.
for i in range(0, len(dna)-12):
read_dna = dna[i:i+12]
for j in rc_len:
cal_rc(read_dna, j, i+1)
for i in range(0,9):
read_dna = dna[len(dna)-12+i:]
pos = len(dna)-11+i
if i == 0:
rc_len = [4,6,8,10,12]
for j in rc_len:
cal_rc(read_dna, j, pos)
elif i >= 1 and i < 3:
rc_len = [4,6,8,10]
for j in rc_len:
cal_rc(read_dna, j, pos)
elif i >=3 and i < 5:
rc_len = [4,6,8]
for j in rc_len:
cal_rc(read_dna, j, pos)
elif i >= 5 and i < 7:
rc_len = [4,6]
for j in rc_len:
cal_rc(read_dna, j, pos)
else:
rc_len = [4]
for j in rc_len:
cal_rc(read_dna, j, pos)