풀이 날짜 : 2023-06-19(R)/ 2023-07-05(python)
| 문제
consensus 서열이란 여러 리드에서 제일 빈도수가 많은 염기가 표기된 서열이다.
이는 시퀀싱 단계에서 일어날 수 있는 에러를 보정하기 위한 방법으로, 가장 빈도수가 높은 염기를 그 자리의 consensus 서열로 지정한다.
| 데이터와 결과
Sample Dataset
>Rosalind_1
ATCCAGCT
>Rosalind_2
GGGCAACT
>Rosalind_3
ATGGATCT
>Rosalind_4
AAGCAACC
>Rosalind_5
TTGGAACT
>Rosalind_6
ATGCCATT
>Rosalind_7
ATGGCACT
Sample Output
ATGCAACT
A: 5 1 0 0 5 5 0 0
C: 0 0 1 4 2 0 6 1
G: 1 1 6 3 0 1 0 0
T: 1 5 0 0 0 1 1 6
아니 output이 서열 하나가 아니라 각 자리의 빈도수까지 출력해야 한다.
사실 각 base마다 자르고 -> 개수 새고 -> 개수가 가장 큰 염기를 consensus서열에 추가하면 된다.
생각보다 간단하겠다 하고 파일을 열어보았는데, 아래와 같이 개행문자로 나누어져 있는 파일이었다.
>Rosalind_3920
GTAATCTCAAATTGCCGCCGGCACAATGATCGAGCTCGCCCAGCGCACCGGTGGTGCTAC
TCTCCGCGCATCCCGACCGCCCGGTGATGCGGTGCGAAGTGGAGGACACGTGGGTTTTGC
CATGGACCGCATTCTGCAGAACATTGCTGTAACTGGATCTAAGATGCTGGAACACTAGTA
CATTAACACAAGCTCCACTCGGCTCCTTGAACCTGATTACCTTGACAACCGTATTTGTGG
CTTTCCCTTGTATTGGTGCTGCCGGCTAGGCCGCCGTAGTTGAAATGTGGCGATGTTTTC
GCTGGTGCGTGCCTCAAGCGAGATTAGGGGCGTTGTACTTGGCTTAGAATTCCCGCAGCA
GACAGTTATCAACGGGAGGGCAACAATGGTGGCCCTAAAGGTACCGGTACTCTGCTGACG
GGAGAGGCCTAGCCGTATATTGCCGAACTTAGCTCCATCGAGGTAAGCTGCCCAGTAGGC
GTACGCATACACGCCTGGCTCATTCAGAAATAATAGTGTAACCCGGCGACTAAATGCTAT
AATACTCCAGTATTACGGGAAATAATATAGTCACCCAATCTATGCCCACCCACCAAATAA
TGATTAAAGCCAACGACGTATCAGGTCGCCTAAGACTGATGATTTTTACGTTTGCTAAAC
TTTATACATGTGAAGTGACGACCTCACAATAATGCCGAGCGAAATGCATGCCACTTTCCG
GAGCGAGACCGATATGGAGTACTTGATGGACGTGGGCTCTTTGGTCAGCTAGCTGGGGTC
TCGCTTACAAAATGAGCAAGGTTTCTTGCCTAGTAATTTAAAATGGGCTCGTAGGGTTAC
CTGTTAGGGTGTAGACACCCGAGGGCAACTTGGTAAAAGGGGAGAGCATGGTGATCCGGT
AGGGCAGTTAAGAAT
>Rosalind_6360
GGGGGAGAATAATGGGCTTTACATCATGTCAAACACGCCCTAGACGATCATCGGGCGTAG
TTGCTCTAGGTCTCGCGCTATTATGTTTGGCCAGCTGTGGATAGTCTATCGTAAACATCT
ATCGATAGGGGTACGGGCTCGGAGACCCATCGTCACTGGAGAGATCTTCGACCGGCCGAA
CGAACGGTTCCTGCTTGCACGGTTTTGGTCACAAGTGGTAGATGTTTTCTATTTCAGGTC
GCTAACATTATAGTTCATTATGGGGCTCGGAATGTGACGCTGCTAGACAGTAGAGCAACT
CATCTAGCGTCACCCTCCCACGCTAGAGACACCACCCTTCGCGTCTGTGGTGGGATTGAG
CTTATAGATCCCTACTGTGGAATCTCCAACCCTTTCGTCTAGTAAATCGATAATAGGGTA
CCAGTTACTCTACGGTCTGATTTGTCGGTCTCACCAATCTGCACGTGCCGGTTTAGTCTG
GAAGCCGGATGACTCGTGGGCAAGTAAGTTGAGCATATGCACTAGGAACCCTGTTACTCC
AACTACGAGATTTGTAAAGACGTTGTATCAGCCAAGAGTACTACCGCGAACGGGGCATAT
TAGCAACCCACCGGGTTAGAAGATATCGAGCTCCCGTATCTGCCTCATTCCAATAGAAGC
TCGAAGAGCGACTTGCAAGCCCGAACTCCCAGACGAAATCCTGGGAGGCTCCCATTGTTG
CGGAGAAAATAAACGGTCGCCGGCGTAGAGACTCTCCGGGGGTTGAGTGGACTCGACGTG
GTGTTTCACACGCCGTGCAAGGCATTACTGTTTTTCGCGTAGTCAGACACCAGCAAAGTG
AGATGCAGATAGTCGAATCCTCTCACTCTCAAGCTTATCAGCGATCCCGTCCAGCGTCGT
| Python
import re, sys, os
import numpy as np
# 파일 읽어오기
with open('./rosalind_cons.txt', 'r') as f :
lines = f.readlines()
# fasta id의 패턴을 인식하는 정규표현식
id_pattern = re.compile('>(\S+)')
# fasta를 dictionary 형태로 저장
Seq_d = {}
for line in lines:
outh = id_pattern.search(line)
if outh:
Seq_id=outh.group(1)
else:
outl=line
if(Seq_id in Seq_d.keys()):
Seq_d[Seq_id] += line.replace('\n', '')
else:
Seq_d[Seq_id] =line.replace('\n', '')
# 서열 길이 파악
seq_len = len(list(Seq_d.values())[0])
print(seq_len) # 976
key_len = len(Seq_d.keys())
print(key_len) # 976
# Numpy를 사용해서 array만들기
value_list = list(Seq_d.values())
array_list = []
for line_string in value_list:
line_character= list(line_string)
array_list.append(line_character)
seq_array = np.array(array_list).T
# 각 열마다 염기 빈도수 확인하고 결과 출력하기
from collections import Counter
element = np.zeros([4, seq_len], dtype=int)
consensus = []
for i in list(range(seq_len)) :
a = Counter(seq_array[i, :])
element[0, i] = a['A']
element[1, i] = a['C']
element[2, i] = a['G']
element[3, i] = a['T']
consensus.append(max(a, key = a.get))
print("".join(consensus))
print("A:", *element[0,:])
print("C:", *element[1,:])
print("G:", *element[2,:])
print("T:", *element[3,:])R이 더 익숙해서, R을 이용해 문제를 푼 후 이 방식을 활용하여 어떻게 Python에서 풀 수 있을지 고민하였다.
일단 정규식을 사용하는 것이 편하고, 관련 예제도 있기에 re 패키지를 조금 공부하고, R의 matrix 같은 array데이터를 다루는 Numpy를 알면 훨씬 편할 것 같아서 이에 대한 강의를 본 후에 풀어보니 이틀 정도 걸렸던 거 같다.
물론 1년 뒤에 푼다면 더 좋은 풀이를 쓸 수도 있을 것이다. 하지만 지금 실력으로는 이 정도가 최선이었다.
| R
library(dplyr)
library(stringr)
# 파일 읽어오기
content <- paste(readLines("./rosalind_cons.txt"), collapse = "")
# >Rosalind_XXXX로 구분하기
pattern <- ">Rosalind_[0-9]{4}"
Read_name <- str_extract_all(content, pattern)[[1]]
Read <- str_split(content, ">Rosalind_[0-9]{4}")[[1]]
Read <- Read[-1]
# 각 Read 쪼개주기
seqs <- Read %>% strsplit("") %>% unlist %>% matrix( nrow = length(Read_name), byrow = TRUE)
read_length <- nchar(Read[1])
# Cosensus 서열 담을 vector
output <- vector(length = read_length)
# 빈도수를 담을 matrix만들기
element <- c("A", "C", "G", "T")
elemet_c <- matrix(nrow = 4, ncol = read_length)
rownames(elemet_c) <- element
matrix <- seqs
# 반복 반복
for (i in 1:read_length) {
base <- matrix[,i]
elemet_c["A", i] <- sum(base == "A")
elemet_c["C", i] <- sum(base == "C")
elemet_c["G", i] <- sum(base == "G")
elemet_c["T", i] <- sum(base == "T")
elemet_c
for (j in element) {
if (elemet_c[j,i] == max(elemet_c[,i])){
out_leter <- names(elemet_c[j,i] )
}
}
output[i] <- out_leter
}
concensus <- base::paste0(output, collapse="")
elemet_c
# 텍스트 파일에 문자열 작성하기
Lines <- c(concensus,
paste0(c("A:", elemet_c[1,]), collapse=" "),
paste0(c("C:", elemet_c[2,]), collapse=" "),
paste0(c("G:", elemet_c[3,]), collapse=" "),
paste0(c("T:", elemet_c[4,]), collapse=" "))
file <- file("./output_CAP.txt", "w")
writeLines(Lines, file)
close(file)
# 아니면 R에서 그냥 출력
cat(Lines, sep = "\n")
- ">Rosalind_[0-9]{4}"은 정규화 패턴으로 문자열에서 여러 패턴을 찾기 위해 사용된다.
- [0-9]는 숫자를 의미하며, {4}는 이 숫자가 4번 반복된다는 의미이다.
- 리눅스, R, Python 모두에서 사용가능하며, 실제 데이터 분석에서 활용도가 높다 (어려워서 안 쓸 뿐이다)
- R은 Python처럼 출력되지 않기 때문에 WriteLines를 통해서 결과를 작성하였다.
Python의 len은 R의 nchar이고
R의 length는 Python의 string.count이다.
각 함수가 정리된 비교표를 찾아봐야겠다.
반응형
