[생물정보학] Biopython 기초

Biopython은 생물정보학에서 가장 널리 사용되는 파이썬 라이브러리입니다.

시퀀스 처리, 다양한 파일 형식 파싱, NCBI 데이터베이스 접근, BLAST 실행 등 생물정보학 분석에 필요한 기능을 폭넓게 제공합니다.

pip install biopython
# 또는
conda install -c bioconda biopython

---

Seq: 시퀀스 객체

Seq는 DNA, RNA, 단백질 시퀀스를 다루는 Biopython의 기본 객체입니다. 파이썬 문자열과 유사하게 사용할 수 있지만, 생물학적 연산 메서드가 추가되어 있습니다.

from Bio.Seq import Seq

# 시퀀스 생성
dna = Seq("ATCGATCGGGCC")
rna = Seq("AUCGAUCGGGCC")
protein = Seq("MADVKKLGEL")

# 기본 조작
print(len(dna))           # 12
print(dna[0:4])           # ATCG (슬라이싱)
print(dna.count("G"))     # 3
print(dna + Seq("TTTT"))  # 시퀀스 연결

# 대소문자
print(dna.upper())
print(dna.lower())

상보서열과 역상보서열

dna = Seq("ATCGATCGGGCC")

print(dna.complement())         # TAGCTAGCCCGG
print(dna.reverse_complement()) # GGCCCGATCGAT

# RNA로 전사
rna = dna.transcribe()
print(rna)    # AUCGAUCGGGCC

# 역전사
dna_back = rna.back_transcribe()
print(dna_back)    # ATCGATCGGGCC

번역 (Translation)

dna = Seq("ATGAAACCCGGGTTTTAA")

# DNA를 단백질로 번역
protein = dna.translate()
print(protein)    # MKPGf* (* = 종결 코돈)

# 종결 코돈까지만 번역
protein = dna.translate(to_stop=True)
print(protein)    # MKPGF

# 다른 유전 코드 사용 (예: 미토콘드리아)
protein = dna.translate(table=2, to_stop=True)

GC 함량 계산

from Bio.SeqUtils import gc_fraction

dna = Seq("ATCGATCGGGCC")
gc = gc_fraction(dna)
print(f"GC 함량: {gc:.2%}")    # GC 함량: 58.33%

---

SeqRecord: 시퀀스 + 메타데이터

SeqRecord는 Seq 객체에 ID, 이름, 설명, 어노테이션 등 메타데이터를 추가한 객체입니다.

파일에서 읽어온 레코드는 대부분 SeqRecord 형태입니다.

from Bio.SeqRecord import SeqRecord
from Bio.Seq import Seq

record = SeqRecord(
    Seq("ATCGATCGGGCC"),
    id="seq_001",
    name="test_sequence",
    description="Homo sapiens BRCA1 exon 1"
)

print(record.id)             # seq_001
print(record.description)   # Homo sapiens BRCA1 exon 1
print(record.seq)            # ATCGATCGGGCC
print(len(record))           # 12

---

SeqIO: 파일 읽기·쓰기

SeqIO는 다양한 시퀀스 파일 형식을 읽고 쓰는 모듈입니다. FASTA, FASTQ, GenBank, EMBL 등을 지원합니다.

파일 읽기

from Bio import SeqIO

# FASTA 읽기
for record in SeqIO.parse("sequences.fasta", "fasta"):
    print(record.id)
    print(len(record.seq))
    print(record.description)

# FASTQ 읽기
for record in SeqIO.parse("reads.fastq", "fastq"):
    print(record.id)
    print(str(record.seq))
    # 품질 점수 (Phred score)
    qualities = record.letter_annotations["phred_quality"]
    print(f"평균 품질: {sum(qualities)/len(qualities):.1f}")

# GenBank 읽기
for record in SeqIO.parse("sequence.gb", "genbank"):
    print(record.id)
    print(record.description)
    # 어노테이션된 특성 목록
    for feature in record.features:
        print(feature.type, feature.location)

 

파일에 레코드가 하나만 있다면 SeqIO.read()를 사용합니다.

record = SeqIO.read("single.fasta", "fasta")

파일 쓰기

from Bio.SeqRecord import SeqRecord
from Bio.Seq import Seq
from Bio import SeqIO

records = [
    SeqRecord(Seq("ATCGATCGGGCC"), id="seq_001", description=""),
    SeqRecord(Seq("GCTAGCTAGCTA"), id="seq_002", description=""),
]

# FASTA로 저장
SeqIO.write(records, "output.fasta", "fasta")

# FASTQ로 저장 (품질 점수 필요)
record = SeqRecord(
    Seq("ATCGATCG"),
    id="read_001",
    description="",
    letter_annotations={"phred_quality": [30, 35, 28, 32, 30, 35, 28, 32]}
)
SeqIO.write([record], "output.fastq", "fastq")

형식 변환

# FASTQ를 FASTA로 변환
count = SeqIO.convert("reads.fastq", "fastq", "reads.fasta", "fasta")
print(f"{count}개 레코드 변환 완료")

딕셔너리로 인덱싱

레코드 수가 적을 때 ID로 빠르게 접근하려면 딕셔너리로 변환합니다.

# 메모리에 전체 로드
seq_dict = SeqIO.to_dict(SeqIO.parse("sequences.fasta", "fasta"))
print(seq_dict["seq_001"].seq)

# 대용량 파일: 인덱스 파일 생성 (메모리 절약)
seq_index = SeqIO.index("large.fasta", "fasta")
print(seq_index["seq_001"].seq)
seq_index.close()

---

Entrez: NCBI 데이터베이스 접근

Entrez 모듈로 NCBI 데이터베이스에서 시퀀스, 논문, 유전체 정보를 검색하고 다운로드할 수 있습니다.

from Bio import Entrez

# 반드시 이메일 등록 (NCBI 정책)
Entrez.email = "your_email@example.com"

# 유전자 검색
handle = Entrez.esearch(db="nucleotide", term="BRCA1[Gene] AND Homo sapiens[Organism]")
record = Entrez.read(handle)
handle.close()

print(f"검색 결과 수: {record['Count']}")
print(f"ID 목록: {record['IdList'][:5]}")

# 시퀀스 다운로드
accession = "NM_007294"    # BRCA1 mRNA
handle = Entrez.efetch(db="nucleotide", id=accession,
                        rettype="fasta", retmode="text")
seq_record = SeqIO.read(handle, "fasta")
handle.close()

print(seq_record.id)
print(len(seq_record.seq))

NCBI는 초당 3회 이상의 요청을 제한합니다. 대량 요청 시 time.sleep(0.5)으로 요청 간격을 두는 것이 좋습니다.

---

BLAST: 시퀀스 유사성 검색

원격 BLAST (NCBI 서버)

from Bio.Blast import NCBIWWW, NCBIXML

query = Seq("ATGAAACCCGGGTTTTAA")

# BLAST 실행 (시간이 걸릴 수 있음)
result_handle = NCBIWWW.qblast("blastn", "nt", query)

# 결과 파싱
blast_records = NCBIXML.parse(result_handle)
blast_record = next(blast_records)

for alignment in blast_record.alignments[:3]:    # 상위 3개
    print(f"Subject: {alignment.title[:60]}")
    hsp = alignment.hsps[0]
    print(f"  Score: {hsp.score}, E-value: {hsp.expect:.2e}")
    print(f"  Identity: {hsp.identities}/{hsp.align_length}")

로컬 BLAST 결과 파싱

BLAST를 로컬에서 실행한 뒤 XML 결과를 파싱합니다.

from Bio.Blast import NCBIXML

with open("blast_result.xml") as f:
    blast_records = NCBIXML.parse(f)
    for record in blast_records:
        for alignment in record.alignments:
            hsp = alignment.hsps[0]
            if hsp.expect < 1e-5:    # E-value 필터링
                print(f"{alignment.title[:50]}: E={hsp.expect:.2e}")

---

Multiple Sequence Alignment

from Bio import AlignIO
from Bio.Align import MultipleSeqAlignment

# 정렬 파일 읽기
alignment = AlignIO.read("alignment.fasta", "fasta")
print(f"시퀀스 수: {len(alignment)}")
print(f"정렬 길이: {alignment.get_alignment_length()}")

for record in alignment:
    print(f"{record.id}: {str(record.seq)[:30]}...")

---

주요 모듈 정리표

모듈	주요 기능
Bio.Seq	시퀀스 객체, 상보서열, 번역
Bio.SeqRecord	시퀀스 + 메타데이터
Bio.SeqIO	시퀀스 파일 읽기·쓰기·변환
Bio.SeqUtils	GC 함량 등 시퀀스 유틸리티
Bio.Entrez	NCBI 데이터베이스 접근
Bio.Blast	BLAST 실행 및 결과 파싱
Bio.AlignIO	다중 서열 정렬 파일 처리
Bio.Phylo	계통수 분석 및 시각화