EMBOSS에서의 서열정렬

EMBOSS의 needle과 water를 이용한 서열 정렬

EMBOSS(European Molecular Biology Open Software Suite)는 다양한 생물정보학 분석 기능을 제공하는 소프트웨어 모음입니다. 이 가운데 needle과 water는 두 서열을 비교하는 대표적인 정렬 도구입니다.

needle은 Needleman–Wunsch 알고리즘을 이용하여 전역 정렬(Global Alignment)을 수행합니다. 두 서열의 전체 길이를 기준으로 최적의 정렬을 찾기 때문에, 길이와 구조가 유사한 서열을 비교할 때 적합합니다.

반면 water는 Smith–Waterman 알고리즘을 이용하여 지역 정렬(Local Alignment)을 수행합니다. 두 서열 전체가 아니라 가장 유사한 부분을 찾아 정렬하므로, 일부 영역만 보존된 서열을 비교할 때 유용합니다.

두 프로그램은 명령줄 인터페이스가 유사하므로 하나의 사용법을 익히면 다른 프로그램에도 쉽게 적용할 수 있습니다.

입력 서열 준비

다음은 사람의 알파 글로빈 단백질 서열을 FASTA 형식으로 저장한 예입니다.

>HBA_HUMAN
MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLSHGSAQVKGHG
KKVADALTNAVAHVDDMPNALSALSDLHAHKLRVDPVNFKLLSHCLLVTLAAHLPAEFTP
AVHASLDKFLASVSTVLTSKYR

이 서열을 alpha.faa 파일에 저장합니다.

비교할 두 번째 서열은 사람의 베타 글로빈 단백질 서열입니다.

>HBB_HUMAN
MVHLTPEEKSAVTALWGKVNVDEVGGEALGRLLVVYPWTQRFFESFGDLSTPDAVMGNPK
VKAHGKKVLGAFSDGLAHLDNLKGTFATLSELHCDKLHVDPENFRLLGNVLVCVLAHHFG
KEFTPPVQAAYQKVVAGVANALAHKYH

이 서열은 beta.faa 파일에 저장합니다.

needle 명령 객체 생성

Biopython에서는 Bio.Emboss.Applications 모듈의 NeedleCommandline 클래스를 이용하여 EMBOSS의 needle 프로그램을 실행할 수 있습니다.

from Bio.Emboss.Applications import NeedleCommandline

needle_cline = NeedleCommandline(
    asequence="alpha.faa",
    bsequence="beta.faa",
    gapopen=10,
    gapextend=0.5,
    outfile="needle.txt"
)

print(needle_cline)

출력되는 명령은 다음과 같습니다.

needle -outfile=needle.txt -asequence=alpha.faa -bsequence=beta.faa -gapopen=10 -gapextend=0.5

각 옵션의 의미는 다음과 같습니다.

• asequence: 첫 번째 입력 서열
• bsequence: 두 번째 입력 서열
• gapopen: 새로운 공백을 생성할 때 적용하는 페널티
• gapextend: 이미 생성된 공백을 연장할 때 적용하는 페널티
• outfile: 정렬 결과를 저장할 파일

공백 열기 페널티는 새로운 삽입이나 결실을 허용할 때 부과되는 비용이며, 공백 연장 페널티는 이미 존재하는 공백을 한 칸 더 늘릴 때 적용되는 비용입니다. 일반적으로 공백 열기 페널티를 공백 연장 페널티보다 크게 설정합니다.

명령줄에서 직접 실행하기

생성된 명령은 명령 프롬프트나 터미널에서도 직접 실행할 수 있습니다.

needle -outfile=needle.txt -asequence=alpha.faa -bsequence=beta.faa -gapopen=10 -gapextend=0.5

정상적으로 실행되면 alpha.faa와 beta.faa의 전역 정렬 결과가 needle.txt 파일에 저장됩니다.

프로그램이 설치되어 있음에도 다음과 같은 메시지가 나타날 수 있습니다.

command not found

Windows에서는 다음과 같은 메시지가 출력될 수 있습니다.

'needle' is not recognized as an internal or external command

이 오류는 EMBOSS 실행 파일이 운영체제의 PATH 환경 변수에 등록되어 있지 않을 때 발생합니다. 이 경우 EMBOSS 설치 경로를 PATH에 추가하거나, 실행 파일의 전체 경로를 직접 지정해야 합니다.

from Bio.Emboss.Applications import NeedleCommandline

needle_cline = NeedleCommandline(
    r"C:\EMBOSS\needle.exe",
    asequence="alpha.faa",
    bsequence="beta.faa",
    gapopen=10,
    gapextend=0.5,
    outfile="needle.txt"
)

Windows 경로에는 역슬래시(\)가 포함되므로 문자열 앞에 r을 붙여 원시 문자열(Raw String)로 지정하는 것이 안전합니다.

속성을 개별적으로 설정하기

명령 객체를 먼저 생성한 다음 각 옵션을 속성으로 지정할 수도 있습니다.

from Bio.Emboss.Applications import NeedleCommandline

needle_cline = NeedleCommandline()

needle_cline.asequence = "alpha.faa"
needle_cline.bsequence = "beta.faa"
needle_cline.gapopen = 10
needle_cline.gapextend = 0.5
needle_cline.outfile = "needle.txt"

print(needle_cline)

생성되는 명령은 앞의 예제와 동일합니다.

needle -outfile=needle.txt -asequence=alpha.faa -bsequence=beta.faa -gapopen=10 -gapextend=0.5

설정된 속성은 다음과 같이 확인할 수 있습니다.

print(needle_cline.outfile)

실행 결과는 다음과 같습니다.

needle.txt

이 방식은 프로그램 실행 조건을 단계적으로 설정하거나, 특정 옵션만 동적으로 변경해야 할 때 유용합니다.

Python에서 needle 실행하기

명령 객체를 실행하려면 다음과 같이 호출합니다.

stdout, stderr = needle_cline()

실행이 정상적으로 완료되면 stdout에는 프로그램의 일반적인 실행 메시지가, stderr에는 경고나 오류 메시지가 저장됩니다.

print(stdout + stderr)

출력 예시는 다음과 같습니다.

Needleman-Wunsch global alignment of two sequences

정렬 결과 자체는 outfile로 지정한 needle.txt 파일에 저장됩니다.

AlignIO로 정렬 결과 읽기

EMBOSS 형식으로 저장된 결과는 Bio.AlignIO 모듈을 이용하여 읽을 수 있습니다.

from Bio import AlignIO

alignment = AlignIO.read(
    "needle.txt",
    "emboss"
)

print(alignment)

출력 예시는 다음과 같습니다.

Alignment with 2 rows and 149 columns
MV-LSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTY...KYR HBA_HUMAN
MVHLTPEEKSAVTALWGKV--NVDEVGGEALGRLLVVYPWTQRF...KYH HBB_HUMAN

정렬 결과에서 하이픈(-)은 삽입 또는 결실에 의해 생긴 공백을 나타냅니다. 정렬 길이가 원래 두 서열보다 길어진 이유는 최적의 대응 관계를 찾는 과정에서 공백이 추가되었기 때문입니다.

결과를 표준 출력으로 받기

정렬 결과를 파일에 저장하지 않고 표준 출력(stdout)으로 직접 받을 수도 있습니다. 임시 파일을 생성하지 않고 Python 메모리에서 결과를 처리하려는 경우 유용합니다.

이때 outfile 대신 stdout=True를 사용합니다.

from Bio.Emboss.Applications import NeedleCommandline

needle_cline = NeedleCommandline(
    asequence="alpha.faa",
    bsequence="beta.faa",
    gapopen=10,
    gapextend=0.5,
    stdout=True
)

stdout, stderr = needle_cline()

반환된 결과는 StringIO를 이용하여 AlignIO에서 읽을 수 있습니다.

from io import StringIO
from Bio import AlignIO

alignment = AlignIO.read(
    StringIO(stdout),
    "emboss"
)

print(alignment)

이 방식은 소규모 정렬이나 자동화된 분석 과정에서 편리합니다. 다만 결과가 매우 큰 경우에는 전체 문자열이 메모리에 저장되므로 파일 기반 처리 방식이 더 안정적일 수 있습니다.

표준 입력 사용하기

EMBOSS의 needle과 water는 입력 서열 중 하나를 표준 입력(stdin)으로 받을 수도 있습니다.

예를 들어 첫 번째 서열을 표준 입력으로 전달하려면 다음과 같이 설정합니다.

needle_cline = NeedleCommandline(
    asequence="stdin",
    bsequence="beta.faa",
    gapopen=10,
    gapextend=0.5,
    stdout=True
)

이후 FASTA 형식의 문자열을 stdin 매개변수로 전달할 수 있습니다.

fasta_data = """>HBA_HUMAN
MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTY
"""

stdout, stderr = needle_cline(
    stdin=fasta_data
)

파일을 생성하지 않고 메모리에 있는 서열을 직접 정렬할 수 있다는 장점이 있지만, 데이터 크기가 크면 메모리 사용량이 증가할 수 있습니다.

하나의 서열과 여러 서열 비교하기

두 번째 입력 파일에 여러 개의 서열이 포함되어 있으면 EMBOSS는 첫 번째 서열과 각 서열을 차례로 비교합니다.

예를 들어 beta.faa 파일에 다섯 개의 서열이 들어 있다면, needle은 총 다섯 개의 서열 쌍 정렬을 수행합니다.

이 기능은 기준 서열 하나를 여러 후보 서열과 비교하거나, 특정 단백질과 유사한 서열을 일괄적으로 분석할 때 유용합니다.

water를 이용한 지역 정렬

water는 needle과 거의 동일한 방식으로 사용할 수 있습니다. 차이점은 water가 Smith–Waterman 알고리즘에 기반한 지역 정렬을 수행한다는 점입니다.

Biopython에서는 WaterCommandline 클래스를 사용합니다.

from Bio.Emboss.Applications import WaterCommandline

water_cline = WaterCommandline(
    asequence="alpha.faa",
    bsequence="beta.faa",
    gapopen=10,
    gapextend=0.5,
    outfile="water.txt"
)

print(water_cline)

출력되는 명령은 다음과 같습니다.

water -outfile=water.txt -asequence=alpha.faa -bsequence=beta.faa -gapopen=10 -gapextend=0.5

실행 방법과 결과를 읽는 방법은 needle과 동일합니다.

stdout, stderr = water_cline()

from Bio import AlignIO

alignment = AlignIO.read(
    "water.txt",
    "emboss"
)

print(alignment)

needle은 두 서열 전체를 대상으로 정렬하지만, water는 두 서열에서 가장 유사한 구간을 찾아 정렬합니다. 따라서 전체적으로는 차이가 크지만 특정 기능 영역이나 도메인만 유사한 단백질을 비교할 때는 water가 더 적합합니다.

서열 쌍 정렬

서열 쌍 정렬(Pairwise Alignment)은 두 개의 생물학적 서열을 비교하여 문자가 가장 잘 대응되는 배치를 찾는 과정입니다. DNA, RNA 또는 단백질 서열 사이의 유사성을 조사하고, 보존 영역이나 삽입·결실·치환 위치를 확인하는 데 사용됩니다.

Biopython은 서열 쌍 정렬을 위해 다음과 같은 기능을 제공합니다.

• 기존의 Bio.pairwise2 모듈
• Bio.Align 모듈의 PairwiseAligner 클래스

두 방식 모두 전역 정렬과 지역 정렬을 지원하며, 일치 점수, 불일치 점수, 공백 열기 페널티, 공백 연장 페널티 등을 조정할 수 있습니다.

pairwise2는 비교적 간단한 인터페이스와 다양한 사용자 정의 기능을 제공하지만, 최근 Biopython에서는 PairwiseAligner 사용을 권장합니다. PairwiseAligner는 일반적으로 처리 속도가 빠르고, 긴 서열을 다룰 때 더 효율적입니다.

다만 기존 코드와 예제에서는 pairwise2가 여전히 많이 사용되므로 두 방식의 기본 원리를 함께 이해할 필요가 있습니다.

pairwise2를 이용한 서열 정렬

Bio.pairwise2 모듈은 EMBOSS의 needle 및 water와 유사한 전역·지역 정렬 기능을 Python 내부에서 제공합니다.

짧거나 중간 길이의 서열에서는 별도의 외부 프로그램을 실행하지 않고도 정렬을 수행할 수 있다는 장점이 있습니다.

기본적인 전역 정렬

다음 예제에서는 앞에서 사용한 알파 글로빈과 베타 글로빈 서열을 FASTA 파일에서 읽어 전역 정렬을 수행합니다.

from Bio import pairwise2
from Bio import SeqIO

seq1 = SeqIO.read(
    "alpha.faa",
    "fasta"
)

seq2 = SeqIO.read(
    "beta.faa",
    "fasta"
)

alignments = pairwise2.align.globalxx(
    seq1.seq,
    seq2.seq
)

globalxx는 두 서열 전체를 대상으로 전역 정렬을 수행합니다.

함수 이름의 마지막 두 글자는 점수 계산 방식을 나타냅니다.

첫 번째 글자는 문자 일치와 불일치에 대한 점수 방식을 의미하고, 두 번째 글자는 공백 페널티 방식을 의미합니다.

globalxx에서 첫 번째 x는 일치하는 문자에 1점을 부여하고, 불일치하는 문자에는 0점을 부여한다는 의미입니다. 두 번째 x는 공백에 별도의 페널티를 부여하지 않는다는 의미입니다.

따라서 globalxx는 가장 단순한 방식으로 두 서열에서 동일한 문자의 수를 최대화하는 정렬을 찾습니다.

정렬 결과 개수 확인하기

하나의 점수에 대해 여러 개의 최적 정렬이 존재할 수 있습니다. 반환값은 이러한 정렬 결과를 담은 리스트입니다.

print(len(alignments))

예제 데이터에서는 동일한 최고 점수를 갖는 여러 정렬이 생성될 수 있습니다.

첫 번째 정렬 결과는 다음과 같이 확인합니다.

print(alignments[0])

출력되는 정렬 객체에는 다음 정보가 포함됩니다.

• 첫 번째 정렬 서열
• 두 번째 정렬 서열
• 정렬 점수
• 정렬 시작 위치
• 정렬 종료 위치

전역 정렬에서는 일반적으로 시작 위치가 0이고, 종료 위치는 정렬된 전체 서열의 길이에 해당합니다.

정렬 결과를 보기 좋게 출력하기

pairwise2는 정렬 결과를 보기 쉬운 형태로 출력하기 위한 format_alignment() 함수를 제공합니다.

print(
    pairwise2.format_alignment(
        *alignments[0]
    )
)

출력 예시는 다음과 같습니다.

MV-LSPADKTNV---K-A--A-WGKVGAHAG---EY-GA-EALE-RMFLSF----PTTK-TY--F...YR-
|| | | | | | |||| | | ||| | | | | |...|
MVHL-----T--PEEKSAVTALWGKV----NVDE-VG-GEAL-GR--L--LVVYP---WT-QRF...Y-H
Score=72

두 서열 사이의 세로선(|)은 동일한 문자가 정렬된 위치를 나타냅니다. 하이픈(-)은 삽입 또는 결실로 인해 추가된 공백입니다.

이 예제에서는 불일치와 공백에 페널티를 부여하지 않았기 때문에 생물학적으로 지나치게 많은 공백이 포함된 정렬이 만들어질 수 있습니다. 실제 단백질 서열 분석에서는 대체 행렬과 공백 페널티를 함께 사용하는 것이 일반적입니다.

키워드 인수 사용하기

최근 Biopython 버전에서는 함수 인수를 명시적인 키워드 형태로 지정할 수 있습니다.

alignments = pairwise2.align.globalxx(
    sequenceA=seq1.seq,
    sequenceB=seq2.seq
)

키워드 인수를 사용하면 각 값의 의미가 명확해지고, 코드의 가독성이 향상됩니다.

공백 페널티와 대체 행렬 적용하기

보다 현실적인 단백질 정렬을 수행하려면 다음 요소를 함께 고려해야 합니다.

• 아미노산 치환 가능성을 반영한 대체 행렬
• 새로운 공백을 생성할 때 적용하는 공백 열기 페널티
• 기존 공백을 연장할 때 적용하는 공백 연장 페널티

단백질 서열 정렬에서는 PAM 또는 BLOSUM 계열의 대체 행렬이 널리 사용됩니다.

다음 예제는 BLOSUM62 대체 행렬과 공백 페널티를 적용한 전역 정렬입니다.

from Bio import pairwise2
from Bio import SeqIO
from Bio.Align import substitution_matrices

blosum62 = substitution_matrices.load(
    "BLOSUM62"
)

seq1 = SeqIO.read(
    "alpha.faa",
    "fasta"
)

seq2 = SeqIO.read(
    "beta.faa",
    "fasta"
)

alignments = pairwise2.align.globalds(
    seq1.seq,
    seq2.seq,
    blosum62,
    -10,
    -0.5
)

globalds의 의미는 다음과 같습니다.

• global: 전역 정렬
• d: 대체 행렬을 사용하여 일치 및 불일치 점수를 계산
• s: 공백 열기와 공백 연장에 서로 다른 페널티를 적용

여기서는 공백 열기 페널티로 -10, 공백 연장 페널티로 -0.5를 사용합니다.

정렬 결과의 수는 다음과 같이 확인합니다.

print(len(alignments))

첫 번째 정렬 결과를 출력하면 다음과 같습니다.

print(
    pairwise2.format_alignment(
        *alignments[0]
    )
)

출력 예시는 다음과 같습니다.

MV-LSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTY...KYR
|| |.|..|..|.|.|||| ......|............|.......||.
MVHLTPEEKSAVTALWGKV-NVDEVGGEALGRLLVVYPWTQRFF...KYH
Score=292.5

이 결과는 단순히 동일한 문자 수를 세는 globalxx보다 생물학적으로 더 타당한 정렬을 제공합니다.

점(.)은 서로 다른 아미노산이지만 BLOSUM62 행렬에서 어느 정도 유사한 치환으로 평가된 위치를 나타냅니다.

동일한 서열, 대체 행렬, 공백 페널티를 사용하면 EMBOSS의 needle과 같은 점수를 얻을 수 있습니다.

지역 정렬 수행하기

지역 정렬은 함수 이름에서 global 대신 local을 사용합니다.

다음 예제는 두 개의 짧은 단백질 서열에서 가장 유사한 부분을 찾습니다.

from Bio import pairwise2
from Bio.Align import substitution_matrices

blosum62 = substitution_matrices.load(
    "BLOSUM62"
)

alignments = pairwise2.align.localds(
    "LSPADKTNVKAA",
    "PEEKSAV",
    blosum62,
    -10,
    -1
)

첫 번째 결과를 출력합니다.

print(
    pairwise2.format_alignment(
        *alignments[0]
    )
)

출력 예시는 다음과 같습니다.

3 PADKTNV
  |..|..|
1 PEEKSAV
  Score=16

지역 정렬에서는 전체 서열이 아니라 가장 높은 점수를 얻은 일부 구간만 출력됩니다.

출력 왼쪽의 숫자는 정렬된 영역이 각 원본 서열에서 시작되는 위치를 나타냅니다.

정렬되지 않은 영역까지 출력하기

지역 정렬 결과를 출력할 때 정렬되지 않은 앞뒤 영역까지 함께 확인하려면 full_sequences=True 옵션을 사용합니다.

print(
    pairwise2.format_alignment(
        *alignments[0],
        full_sequences=True
    )
)

출력 예시는 다음과 같습니다.

LSPADKTNVKAA
  |..|..|
--PEEKSAV---
Score=16

이 형식은 정렬된 부분이 원본 서열의 어느 위치에 존재하는지 전체적인 맥락에서 확인할 때 유용합니다.

지역 정렬 점수의 특성

Smith–Waterman 방식의 지역 정렬은 양수 점수를 갖는 구간만 결과로 제공합니다.

최고 정렬 점수가 0 이하라면 의미 있는 지역 정렬이 존재하지 않는 것으로 간주되어 결과가 반환되지 않을 수 있습니다.

또한 점수가 0인 문자 대응은 정렬 영역의 시작이나 끝을 확장하는 데 사용되지 않을 수 있습니다. 이는 지역 정렬이 양의 누적 점수를 갖는 핵심 영역을 중심으로 결과를 구성하기 때문입니다.

일치·불일치 점수를 직접 지정하기

대체 행렬을 사용하지 않고 일치와 불일치 점수를 직접 설정할 수도 있습니다.

다음 예제에서는 일치에 5점, 불일치에 -4점, 공백 열기에 -2점, 공백 연장에 -0.5점을 적용합니다.

alignments = pairwise2.align.localms(
    "AGAACT",
    "GAC",
    5,
    -4,
    -2,
    -0.5
)

정렬 결과를 출력합니다.

print(
    pairwise2.format_alignment(
        *alignments[0]
    )
)

출력 예시는 다음과 같습니다.

2 GAAC
  | ||
1 G-AC
  Score=13

localms에서 m은 사용자가 일치 점수와 불일치 점수를 직접 지정한다는 의미이고, s는 공백 열기와 연장에 서로 다른 값을 적용한다는 의미입니다.

점수만 계산하기

실제 정렬 결과가 필요하지 않고 최고 점수만 확인하려면 score_only=True 옵션을 사용할 수 있습니다.

score = pairwise2.align.globalds(
    seq1.seq,
    seq2.seq,
    blosum62,
    -10,
    -0.5,
    score_only=True
)

print(score)

정렬 경로와 정렬 문자열을 생성하지 않기 때문에 실행 속도가 빨라지고 메모리 사용량도 줄어듭니다.

많은 서열 쌍의 유사도를 빠르게 비교하거나, 정렬 점수를 기준으로 후보를 선별할 때 유용합니다.

하나의 최적 정렬만 반환하기

최적 점수를 갖는 정렬이 여러 개 존재하더라도 하나의 결과만 필요하다면 one_alignment_only=True 옵션을 사용할 수 있습니다.

alignments = pairwise2.align.globalds(
    seq1.seq,
    seq2.seq,
    blosum62,
    -10,
    -0.5,
    one_alignment_only=True
)

이 옵션은 불필요한 정렬 결과 생성을 줄여 처리 속도와 메모리 효율을 높일 수 있습니다.

pairwise2 사용 시 고려 사항

pairwise2는 간단한 사용법과 다양한 점수 설정 기능을 제공하지만, 최신 Biopython에서는 더 새로운 PairwiseAligner 클래스를 권장합니다.

특히 긴 서열을 정렬하거나 새로운 프로젝트를 작성할 때는 PairwiseAligner가 더 적합합니다.

pairwise2는 기존 코드와 예제를 이해하거나, 사용자 정의 점수 함수와 특수한 정렬 조건을 적용해야 할 때 유용하게 사용할 수 있습니다.