Biopython
Bio.SeqIO
#
Find similar titles
- 최초 작성자
- 최근 업데이트
Structured data
- Category
- Programming
Bioptython을 활용한 서열 시퀀스를 읽고 쓰는 법(서열 입출력)을 알아보고자 한다. 서열을 읽고 쓰는 핵심은 바로 from Bio import SeqIO 이다. 이는 서열을 읽고 쓰는 여러 함수를 제공한다. 기본 사용은 다음과 같다.
>>> from Bio import SeqIO
>>> help(SeqIO)
...
서열을 읽는 법은 기본적으로 SeqIO 객체를 이용하여 SeqIO.parse()를 이용하면 된다. 이는 크게 기본적으로 두 가지 인자를 필요로 하는데 파일이름과 포맷이다. 포맷은 익히 알려진대로 fasta와 genbank가 있다. 두 가지 인수를 이용하여 우리가 자주 사용하는 파일을 읽을 수 있다.
>>> from Bio import SeqIO
>>> for seq_record in SeqIO.parse("ls_orchid.fasta", "fasta"):
... print(seq_record.id)
... print(repr(seq_record.seq))
... print(len(seq_record))
위의 예제에서는 SeqIO.parse("ls_orchid.fasta", "fasta")를 이용하여 ls_orchid.fasta 라는 파일을 읽어 들이며 이 파일의 기본 포맷은 fasta라고 알려주는 것이다. 이렇게 파일을 읽으면 기본적으로 루프 구조를 가지기 때문에 for 구문을 이용하여 순환하며 각 서열 하나하나를 출력할 수가 있다. 출력에는 for 구문에서 받은 변수에 .xxx 를 이용하여 출력할 수가 있다.
>>> from Bio import SeqIO
>>> for seq_record in SeqIO.parse("ls_orchid.gbk", "genbank"):
... print(seq_record.id)
... print(seq_record.seq)
... print(len(seq_record))
위의 예제에서는 동일한 방법으로 genbank로 된 ls_orchid.gbk 파일을 읽어들이는 방법이다. fasta와 같이 읽어 들였던 방법으로 동일하게 지정하여 파일명과 포맷 형식을 바꿔주면 된다.
>>> from Bio import SeqIO
>>> identifiers = [seq_record.id for seq_record in SeqIO.parse("ls_orchid.gbk", "genbank")]
>>> identifiers
['Z78533.1', 'Z78532.1', 'Z78531.1', 'Z78530.1', 'Z78529.1', 'Z78527.1', ..., 'Z78439.1']
위의 예제에서는 리스트를 이용하여 각 서열을 파싱한 다음 리스트에 담아 저장하는 예제를 보여준다. 자료구조론에 입각하여 다양한 방법을 통해 각 서열을 응용할 수 있다.
>>> from Bio import SeqIO
>>> record_iterator = SeqIO.parse("ls_orchid.fasta", "fasta")
>>> first_record = next(record_iterator)
>>> print(first_record.id)
>>> print(first_record.description)
>>> second_record = next(record_iterator)
>>> print(second_record.id)
>>> print(second_record.description)
위의 예제에서는 서열을 순환하는 또 다른 방법을 보여준다. 기본적으로 SeqIO와 for 구문을 이용하여 순환하여 데이터를 출력할 수 있지만 변수에 값을 입력한 후 next()를 이용하여 다음 서열로 이동할 수 있게 해준다. 이 예제에서는 next()만 잘 기억하고 있으면 된다.
>>> from Bio import SeqIO
>>> records = list(SeqIO.parse("ls_orchid.gbk", "genbank"))
>>> print("Found %i records" % len(records))
다음 예제는 전체 서열을 리스트에 담아 서열의 사이즈를 구하는 예제이다. list() 함수를 이용하여 서열의 각 레코드 객체를 저장하고 리스트에 len()함수를 이용하여 서열의 전체 갯수를 구하는 예제이다. 다른 많은 서열의 사이즈를 구하는 방법이 있지만 위의 예제에서는 리스트를 이용한 방법이다. 갯수가 많아질 경우에는 메모리를 많이 점유하는 단점이 존재한다.
>>> print("The last record")
>>> last_record = records[-1] # using Python's list tricks
>>> print(last_record.id)
>>> print(repr(last_record.seq))
>>> print(len(last_record))
>>> print("The first record")
>>> first_record = records[0] # remember, Python counts from zero
>>> print(first_record.id)
>>> print(repr(first_record.seq))
>>> print(len(first_record))
다음 예제는 파이썬의 리스트를 이용하여 각 원하는 순서의 서열을 출력하는 예제이다. 위의 예제에서는 리스트[0], 리스트1과 같이 파이썬의 자료 구조론을 이용하여 출력하는 방법이다. 필요에 따라서 원하는 형태로 출력하면 된다. 다만 파이썬을 비롯한 많은 언어에서는 리스트는 0이 첫 번째를 지칭한다는 것만 이해하고 있으면 된다.
>>> from Bio import SeqIO
>>> record_iterator = SeqIO.parse("ls_orchid.gbk", "genbank")
>>> first_record = next(record_iterator)
>>> print(first_record)
앞에서 언급하였듯이 리스트의 트릭을 사용하지 않고 next()를 이용하여 출력할 수도 있다. 출력 결과는 다음과 같다.
ID: Z78533.1
Name: Z78533
Description: C.irapeanum 5.8S rRNA gene and ITS1 and ITS2 DNA.
Number of features: 5
/sequence_version=1
/source=Cypripedium irapeanum
/taxonomy=['Eukaryota', 'Viridiplantae', 'Streptophyta', ..., 'Cypripedium']
/keywords=['5.8S ribosomal RNA', '5.8S rRNA gene', ..., 'ITS1', 'ITS2']
/references=[...]
/accessions=['Z78533']
/data_file_division=PLN
/date=30-NOV-2006
/organism=Cypripedium irapeanum
/gi=2765658
Seq('CGTAACAAGGTTTCCGTAGGTGAACCTGCGGAAGGATCATTGATGAGACCGTGG...CGC', IUPACAmbiguousDNA())
다음의 예제에서는 서열의 종 정보를 추출/출력하기 위한 방법에 대해 설명한다. 먼저 종 정보를 담기 위해 하나의 리스트를 먼저 선언해주고 그 다음 각 서열의 종 정보를 추가하여 출력하는 예제이다.
>>> from Bio import SeqIO
>>> all_species = []
>>> for seq_record in SeqIO.parse("ls_orchid.gbk", "genbank"):
... all_species.append(seq_record.annotations["organism"])
>>> print(all_species)
즉, all_species=[] 와 같이 먼저 종 정보를 넣을 빈 리스트를 생성한다. 그 다음 SeqIO 객체를 이용하여 각 서열을 순환하여 처리한다. 여기서는 seq_record.annotations["organism"] 와 같이 이미 정해진 annotations라는 사전에 값이 들어 있으므로 annotations["organism"] 와 같이 종 정보를 추출하여 all_species.append()를 이용하여 빈 리스트에 추가하면 된다. for의 순환이 끝나면 all_species에는 전체 서열의 종 정보가 들어 있게 되고 이를 출력하면 된다.
>>> from Bio import SeqIO
>>> all_species = [seq_record.annotations["organism"] for seq_record in SeqIO.parse("ls_orchid.gbk", "genbank")]
>>> print(all_species)
물론 위의 예제와 같이 변수=[ 수식]을 이용하여 한번에 처리할 수도 있다. 위의 두 종 정보를 출력하면 다음과 같이 출력되면 정상이다. 이는 다시 for 구문 순환을 통해서 별도의 데이터로 활용할 수 있다.
['Cypripedium irapeanum', 'Cypripedium californicum', ..., 'Paphiopedilum barbatum']
위의 경우에는 genbank를 이용한 방법이다. 그러나 fasta 포맷에서는 좀 다르다. 아는 바와 같이 fasta에서는 종 정보가 >로 시작하는 설명라인에 포함되어 있으므로 설명라인을 확인하여 종 정보가 있는 위치를 파싱하여 추가하여여한다. 물론 아래의 예제에서는 좀 더 파싱 라인이 추가로 필요할 수 있다.
>>> from Bio import SeqIO
>>> all_species = []
>>> for seq_record in SeqIO.parse("ls_orchid.fasta", "fasta"):
... all_species.append(seq_record.description.split()[1])
>>> print(all_species)
아래 예제는 파일에서 서열을 불러 들이는 것이 아니라 Entrez 기능을 이용하여 서열을 검색 후 서열을 읽어 오는 방법이다. 여기서는 추가적으로 하나의 객체 요소가 더 필요한 데 from Bio import Entrez 이다. 이는 NCBI에서 제공하는 다양한 데이터를 추출하기 위한 메서드를 제공한다. 또한 예전 버전과 달리 메일 주소는 반드시 적어줘야 한다. 물론 임의의 메일 주소를 적어줘도 된다. 다만 메일의 형식은 유지해야 한다. Entrez를 이용하는 방법은 아래와 같이 Entrez.efetch(db="서열 데이터베이스 이름", rettype="포맷", retmode="데이터 타입", id="검색어 아이디")를 사용한다.
>>> from Bio import Entrez
>>> from Bio import SeqIO
>>> Entrez.email = "A.N.Other@example.com"
>>> handle = Entrez.efetch(db="nucleotide", rettype="fasta", retmode="text", id="6273291")
>>> seq_record = SeqIO.read(handle, "fasta")
>>> handle.close()
>>> print("%s with %i features" % (seq_record.id, len(seq_record.features)))
그런 후에 SeqIO.read(handle, "fasta") 와 같이 서열의 입출력 객체인 SeqIO를 이용하여 서열을 정보를 파싱할 수 있다. 다만 여기서 read를 사용한 것은 단 하나의 서열이기 때문이다. 단 하나의 서열인 경우는 read, 둘 이상의 서열은 parse를 이용하면 된다.
>>> from Bio import Entrez
>>> from Bio import SeqIO
>>> Entrez.email = "A.N.Other@example.com"
>>> handle = Entrez.efetch(db="nucleotide", rettype="gb", retmode="text", id="6273291")
>>> seq_record = SeqIO.read(handle, "gb") #using "gb" as an alias for "genbank"
>>> handle.close()
>>> print("%s with %i features" % (seq_record.id, len(seq_record.features)))
위의 예제는 마찬가지로 서열의 포맷 형태를 genbank로 가져오는 것만 다를 뿐 동일한 코드 형태를 취하고 있다.
이처럼 SeqIO 객체를 이용하면 다양한 서열 파일들을 손쉽게 다룰 수 있고 Entrez를 이용하여 NCBI에서 제공하는 많은 정보를 추가적으로 검색하고 다운로드하여 이용할 수 있다. 이를 응용한다면 우리가 다루고자 하는 많은 서열을 좀 더 손쉽게 이용할 수 있을 것이다.
Tip #
너무 서열이 긴 FASTA 포맷의 서열을 SeqIO를 이용하여 읽을 경우 Memory error가 발생한다. FASTA 포맷을 읽는 모듈에서 공백이나 캐리지 리턴('\r') 같은 문자를 제거 할 때 서열을 모두 한 줄로 모은 다음 replace문을 이용하여 변경하는데,이 때 서열이 너무 길면 replace문이 실행될 때 아래와 같은 Memory error가 발생하게 된다.
File "C:\Python27\lib\site-packages\Bio\SeqIO_init_.py", line 600, in parse
for r in i:
File "C:\Python27\lib\site-packages\Bio\SeqIO\FastaIO.py", line 122, in FastaIterator
for title, sequence in SimpleFastaParser(handle):
File "C:\Python27\lib\site-packages\Bio\SeqIO\FastaIO.py", line 67, in SimpleFastaParser
yield title, "".join(lines).replace(" ", "").replace("\r", "")
MemoryError
이 경우 @BIOPYTHON경로@/Bio/SeqIO/FastaIO.py에서 각 서열 라인마다 replace문을 실행하고 모으도록 변경하면 해결된다.
-
파일경로 확인
$ python Python 2.6.6 (r266:84292, Jan 22 2014, 09:42:36) [GCC 4.4.7 20120313 (Red Hat 4.4.7-4)] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> import Bio.SeqIO >>> Bio.SeqIO.__file__ '/usr/lib64/python2.6/site-packages/Bio/SeqIO/__init__.pyc' >>> exit() -
Fasta 포맷 처리 모듈(@BIOPYTHON경로@/Bio/SeqIO/FastaIO.py) 변경
- 이 경우 FastaIO.py의 경로는 '/usr/lib64/python2.6/site-packages/Bio/SeqIO/FastaIO.py'
- 대부분의 경우 일반유저에게는 파일 변경 권한이 없으므로 관리자 계정으로 수정
-
버전에 따라 조금씩 다르지만 50~70줄 사이에 아래 코드 확인
while True: if line[0] != ">": raise ValueError( "Records in Fasta files should start with '>' character") title = line[1:].rstrip() lines = [] line = handle.readline() while True: if not line: break if line[0] == ">": break lines.append(line.rstrip()) line = handle.readline() # Remove trailing whitespace, and any internal spaces # (and any embedded \r which are possible in mangled files # when not opened in universal read lines mode) yield title, "".join(lines).replace(" ", "").replace("\r", "") -
이 부분을 아래와 같이 변경
while True: if line[0] != ">": raise ValueError( "Records in Fasta files should start with '>' character") title = line[1:].rstrip() lines = [] line = handle.readline() while True: if not line: break if line[0] == ">": break #lines.append(line.rstrip()) lines.append(line.rstrip().replace(" ", "").replace("\r", "")) line = handle.readline() # Remove trailing whitespace, and any internal spaces # (and any embedded \r which are possible in mangled files # when not opened in universal read lines mode) #yield title, "".join(lines).replace(" ", "").replace("\r", "") yield title, "".join(lines)
