바이오피슨의 바이오.SeqIO에서 gbff 파일을 처리하는 방법

개시하다

첫 투고.생체정보학과 기사가 거의 보이지 않아 심상치 않게 느껴졌지만, 큐타보다 UI가 더 멋져 젠에 기고하기로 했다.
이 글은biopythonBio.SeqIO을 사용하여Genbank 파일(gbff)을 조작하는 방법에 대해 자신이 자주 사용하는 부분만 설명한다.2022년 들어 한 달 간격으로 이런 일이 벌어진다자기를 위한 해설 기사.

읽다

기본적으로 gz로 압축되어 있기 때문에 해동의 절차를 표시합니다.gzip -d 등으로 미리 해동해도 괜찮아요.

path = "./sample.gbff.gz"

import gzip
from Bio import SeqIO

with gzip.open(path, mode="rt", encoding="utf-8") as f:
    records = list(SeqIO.parse(f, "genbank"))

!
SeqIO.parse () 에서 gbff 파일을 읽을 수 있지만, 이 방법은 교체기를 되돌릴 수 있다는 것을 주의해야 한다.
지정된 gbff 파일이 별도의 항목으로 구성되면 겉으로는 문제가 없지만 MAG(Metagenome assembed genome)를 비롯한 콘텐츠에 대한 gbff에 넥스트() 등을 추가하지 않으면 데이터가 돌아오지 않는다.
여기서 간단하게 보기 위해list()를 통해 모든 데이터를yeld로 변환하지만 메모리의 데이터량을 먼저 먹고 이동전화의 장점을 이용하여 하나하나next를 사용하는 것이 좋다고 생각합니다.
20222.03.17 보충
엔트를 제네레이션이라고 부르는 글이어서 오해가 생길 수 있어 살짝 수정했다.
또한 SeqIO.parse()해 봤는데 with에서 파일을 여는 기간이 아니면 넥스트()를 할 수 없기 때문에 메모리를 관리하기 위해 하나하나 읽는 것이 더욱 주의가 필요하다.

SeqRecord 클래스 객체

읽은 record는 SeqRecord 클래스의 인스턴스입니다.

print(f"{records[0]=}")
print(f"{type(records[0])=}")

records[0]=SeqRecord(seq=Seq('GCGTTAAAACTTAAAGTGGAGGCGCGACCCGGAGTCGAACCGAGATCGACGGAT...ACA'), id='NZ_MCBT01000001.1', name='NZ_MCBT01000001', description='Shewanella colwelliana strain CSB03KR, whole genome shotgun sequence', dbxrefs=['BioProject:PRJNA224116', 'BioSample:SAMN05384437', 'Assembly:GCF_001735525.1'])
type(records[0])=
읽은 SeqRecord 대상은 개인적으로 처리하기 어려워서 나중에 처리할 때 종종 헷갈릴 수 있습니다.
따라서 아래 부분에서 설명한 내용은 총결의 본질일 수 있다.

인스턴스 변수

다음 인스턴스 변수를 사용하여 레코드의 프로파일을 볼 수 있습니다.굳이 말하자면, 마지막 annotations는 매우 중요하고, 그 외에 거의 사용하지 않을 수도 있다.
value
설명
예제
dbxrefs
세션 번호 확인 등
['Bio Project: PRJNA 224116', (생략)]
id
연결 ID
'NZ_MCBT01000001.1'
name
시퀀스 확인 세션
'NZ_MCBT01000001'
description
조합 설명
'Shewanella colwelliana strain CSB03KR, whole genome shotgun sequence'
annotations
상세 시퀀스 정보
{'molecule type':'DNA','topology':'linear',...(생략)

annotations의 하층

방문record.annotation 후 dict 형식으로 매우 많은 정보를 얻을 수 있습니다.
NCBI의 항목이기 때문에 메타 정보가 다른 샘플과 비교되지 않는 경우가 종종 있는데, 어쨌든 현재 샘플 데이터로 사용되고 있는 것들에 포함된 메타 정보를 나열해 보겠습니다.
value
예제
molecule_type
dna
topology
linear
data_file_division
CON
date
05-JUL-2021
accessions
['NZ_MCBT01000001', 'NZ_MCBT01000000']\
sequence_version
1
keywords
['WGS', 'RefSeq']
source
Shewanella colwelliana
taxonomy
['Bacteria','Proteobacteria','Gammaproteobacteria','Alteromonadales','Shewanellaceae','Shewanella']
references
생략하다
comment
생략하다
contig
'join(MCBT01000001.1:1..325457)'
stductured_comment
내포된 객체
예를 들어 taxonomy의 정보를 보고 싶을 때 dict의 키로 부른다.

records[0].annotations["taxonomy"]

SeqFeature 클래스 객체

SeqFeature류는 유전인자 정보를 기록한 것이다.염색체 1의start,end의 영역(strand는 +/-), 야옹이 기능을 가진 유전자인 등의 정보가 존재한다.
아까SeqRecord류에 대해SeqRecord.features하면 리스트 형식으로 빠르게 데이터가 나온다.
예를 들어 피처에서 rRNA만 꺼내려면 이렇게 하세요.

res = []
for record in records:
    minires = []
    for feature in record.features:
        if feature.type == "rRNA":
            minires.append(feature)
        else:
            pass
    if minires:
        res.append(minires)
    else:
        pass

res

[[SeqFeature(FeatureLocation(ExactPosition(149124), ExactPosition(150667), strand=1), type='rRNA')],
[SeqFeature(FeatureLocation(ExactPosition(241270), ExactPosition(241386), strand=-1), type='rRNA')],
[SeqFeature(FeatureLocation(ExactPosition(11), ExactPosition(1554), strand=1), type='rRNA'),
SeqFeature(FeatureLocation(ExactPosition(1909), ExactPosition(4803), strand=1), type='rRNA')],
[SeqFeature(FeatureLocation(ExactPosition(11), ExactPosition(1554), strand=1), type='rRNA'),
SeqFeature(FeatureLocation(ExactPosition(2403), ExactPosition(5297), strand=1), type='rRNA')],
[SeqFeature(FeatureLocation(BeforePosition(0), AfterPosition(516), strand=-1), type='rRNA')],
[SeqFeature(FeatureLocation(ExactPosition(319), AfterPosition(637), strand=-1), type='rRNA')],
[SeqFeature(FeatureLocation(ExactPosition(216), ExactPosition(332), strand=1), type='rRNA'),
SeqFeature(FeatureLocation(ExactPosition(497), ExactPosition(613), strand=1), type='rRNA')]]

SeqFeature가 유지하는 유전적 인자 정보

여기서 얻은 SeqFeature류의 대상은 언뜻 보기에는 단순한 데이터형을 구성하지만, print()는 딱 봐도 다양한 정보가 담겨 있다.

print(res[0][0])

type: rRNA
location: 149124:150667
qualifiers:
Key: db_xref, Value: ['RFAM:RF00177']
Key: inference, Value: ['COORDINATES: nucleotide motif:Rfam:12.0:RF00177', 'COORDINATES: profile:INFERNAL:1.1.1']
Key: locus_tag, Value: ['BEL05_RS03635']
Key: note, Value: ['Derived by automated computational analysis using gene prediction method: cmsearch.']
Key: old_locus_tag, Value: ['BEL05_09200']
Key: product, Value: ['16S ribosomal RNA']
sequence의 위치와 위치 (방향) 는 다음과 같습니다.

print(res[0][0].location.start)
print(res[0][0].location.end)
print(res[0][0].location.strand)

149124
150667
1
qualifiers는 dict로 정보를 저장합니다.
나는 자주'제품'을 방문한다.다른 모티프 정보 등을 활용할 수 있을지 연구해 본 적이 없어 모르겠다.
자세한 내용은 공식 문서를 보십시오.

res[0][0].qualifiers['product']

['16S ribosomal RNA']

Seq 클래스 객체

seq 대상은 수조 정보 자체이지만, 그 행위가str류와 같도록 특수한 방법을 정의합니다.

records[0].seq

Seq('GCGTTAAAACTTAAAGTGGAGGCGCGACCCGGAGTCGAACCGAGATCGACGGAT...ACA')

>>> len(records[0].seq)
325457
>>> str(records[0].seq)[:10]
'GCGTTAAAAC'
>>> records[0].seq[:10]
'GCGTTAAAAC'
>>> records[0].seq[:10] + "AAAAAAAAAAAAAAAA"
Seq('GCGTTAAAACAAAAAAAAAAAAAAAA')
>>> records[0].seq.count("AAG")
5692

행위와str류는 기본적으로 같은 것이 좋다. 그러나 이 정도의 조작만 하면 Seq류는 좋은 점이 없다. 나는 바로str류로 전환했다.
아마 아미노산 번역 등을 하면 Seq반의 좋은 점이 있을 거라고 생각했는데, 지금은 그런 일이 없어서 테스트를 안 했어요.

최후

아주 간단한 메모 정도의 정보만 있을 뿐이지만 적어도 미래의 나는 얻을 수 있을 것이다.
여기까지 읽은 너를 위해서라면, 나는 매우 기쁠 것이다!