Data Science/Bioinformatics with Biopython

8. Chapter7. Multiple Sequence Alignment

HJChung 2020. 3. 20. 19:13

 

 

 

 

 

다중 서열 정렬: Multiple Sequence Alignment

Multiple Sequence Alignment객체는 각 행들이 SeqRecord 객체로 된 행렬과 같은 객체입니다. 여기서는

  1. sequence alignment결과 파일을 읽고/ 파싱하고/서열 중간을 슬라이싱하는 방법에 대해 알아봅니다.
  2. 직접 서열들을 가지고 정렬을 진행해 봅니다.
  3. 보존 서열을 알 수 있는 WebLogo를 그리고 종 간 유사도를 알아볼 수 있는 계통수를 그려봅니다.

1. Multiple Sequence Alignment 준비 과정

Multiple Sequence Alignment를 하기 위해서는

  1. 비교할 여러 서열이 필요합니다. : 여기서는 헤모글로빈의 HBA단백질에 대해 인간을 비롯한 다른 동물들의 HBA서열을 단백질 데이터베이스(www.uniprot.org)에서 받아옵니다.
  2. Multiple Sequence Alignment과정을 진행합니다. : 여기서는 Alignment를 수행하는 tool로, MUSCLE을 사용합니다.
  3. 위의 과정을 시각화해 봅니다. : 여기서는 WebLogo와 계통수를 그려봅니다.
 

2.sequence alignment결과 파일을 읽고/ 파싱하고/서열 중간을 슬라이싱하는 방법에 대해 알아봅니다.

이를 위해서는 바이오파이썬의 AlignIO모듈을 사용합니다. 그리고 파일을 읽을 때는 여태까지의 경우와 마찬가지로 AlignIO.read('파일이나 파일을 연 객체', '파일 종류')를 넣어줍니다.

In [1]:
#1. read MSA 
from Bio import AlignIO

alignment = AlignIO.read("HBA.aln", "clustal")
print(alignment)
 
SingleLetterAlphabet() alignment with 7 rows and 142 columns
MVLSAADKNNVKGIFTKIAGHAEEYGAETLERMFTTYPPTKTYF...KYR sp|P01994|HBA_CHICK
MVLSPTDKSNVKATWAKIGNHGAEYGAEALERMFMNFPSTKTYF...KYR sp|P18971|HBA_BALAC
MVLSPADKTNIKTAWEKIGSHGGEYGAEAVERMFLGFPTTKTYF...KYR sp|P01948|HBA_RABIT
MVLSGEDKSNIKAAWGKIGGHGAEYGAEALERMFASFPTTKTYF...KYR sp|P01942|HBA_MOUSE
MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYF...KYR sp|P69907|HBA_PANTR
MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYF...KYR sp|P69905|HBA_HUMAN
MVLSAADKTNVKAAWSKVGGHAGEYGAEALERMFLGFPTTKTYF...KYR sp|P01958|HBA_HORSE
In [2]:
#2. 읽어온 Alignment 다루기 - slice sequence나누기)
#위의 MSA 객체에서 서열 부분과 ID부분을 나누고자 한다면
for record in alignment:
    print("%s-%s" %(record.seq, record.id))
 
MVLSAADKNNVKGIFTKIAGHAEEYGAETLERMFTTYPPTKTYFPHFDLSHGSAQIKGHGKKVVAALIEAANHIDDIAGTLSKLSDLHAHKLRVDPVNFKLLGQCFLVVVAIHHPAALTPEVHASLDKFLCAVGTVLTAKYR-sp|P01994|HBA_CHICK
MVLSPTDKSNVKATWAKIGNHGAEYGAEALERMFMNFPSTKTYFPHFDLGHDSAQVKGHGKKVADALTKAVGHMDNLLDALSDLSDLHAHKLRVDPANFKLLSHCLLVTLALHLPAEFTPSVHASLDKFLASVSTVLTSKYR-sp|P18971|HBA_BALAC
MVLSPADKTNIKTAWEKIGSHGGEYGAEAVERMFLGFPTTKTYFPHFDFTHGSEQIKAHGKKVSEALTKAVGHLDDLPGALSTLSDLHAHKLRVDPVNFKLLSHCLLVTLANHHPSEFTPAVHASLDKFLANVSTVLTSKYR-sp|P01948|HBA_RABIT
MVLSGEDKSNIKAAWGKIGGHGAEYGAEALERMFASFPTTKTYFPHFDVSHGSAQVKGHGKKVADALASAAGHLDDLPGALSALSDLHAHKLRVDPVNFKLLSHCLLVTLASHHPADFTPAVHASLDKFLASVSTVLTSKYR-sp|P01942|HBA_MOUSE
MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLSHGSAQVKGHGKKVADALTNAVAHVDDMPNALSALSDLHAHKLRVDPVNFKLLSHCLLVTLAAHLPAEFTPAVHASLDKFLASVSTVLTSKYR-sp|P69907|HBA_PANTR
MVLSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHFDLSHGSAQVKGHGKKVADALTNAVAHVDDMPNALSALSDLHAHKLRVDPVNFKLLSHCLLVTLAAHLPAEFTPAVHASLDKFLASVSTVLTSKYR-sp|P69905|HBA_HUMAN
MVLSAADKTNVKAAWSKVGGHAGEYGAEALERMFLGFPTTKTYFPHFDLSHGSAQVKAHGKKVGDALTLAVGHLDDLPGALSNLSDLHAHKLRVDPVNFKLLSHCLLSTLAVHLPNDFTPAVHASLDKFLSSVSTVLTSKYR-sp|P01958|HBA_HORSE
 

위의 결과로부터 알 수 있듯이 sequence alignment파일을 읽으면 여러개의 SeqRecord 객체가 각 행에 들어있습니다. SeqRecord 객체에는 seq속성이 있고 이는 서열을 담고 있는데 문자열과 같이 슬라이싱이 가능합니다. 예를들어, 각 SeqRecord객체의 서열에서 앞 10개의 서열을 가져온다고 하면, [start:end]를 사용하면 됩니다.

In [3]:
#3. 읽어온 Alignment 다루기 - slice sequence나누기)
#SeqRecord객체의 서열에서 앞 10개의 서열을 가져온다고 하면, 
for record in alignment:
    print("%s-%s" %(record.seq[:10], record.id))
 
MVLSAADKNN-sp|P01994|HBA_CHICK
MVLSPTDKSN-sp|P18971|HBA_BALAC
MVLSPADKTN-sp|P01948|HBA_RABIT
MVLSGEDKSN-sp|P01942|HBA_MOUSE
MVLSPADKTN-sp|P69907|HBA_PANTR
MVLSPADKTN-sp|P69905|HBA_HUMAN
MVLSAADKTN-sp|P01958|HBA_HORSE
 

3. 직접 서열들을 가지고 정렬을 진행해 봅니다.

FASTA 파일에서부터 MSA 객체를 만들기 위해 MSA를 수행하는 툴 중 유명한 MUSCLE 툴을 사용합니다.

  1. 우선 MUSCLE을 설치합니다. conda install -c bioconda muscle
  2. 바이오파이썬을 통해 MUSCLE을 실행합니다. MUSCLE은 콘솔 환경에서 단독 실행할 수도 있지만 바이오파이썬을 통해 실행 할 수도 있습니다.
In [4]:
#4. Mutiple Sequence Alignmnet 객체의 서열 앞에서 10개 서열 가져오기 (Muscle을 사용하여)
from Bio.Align.Applications import MuscleCommandline

muscle_exe = "/Users/jeonghyeonjeong/for github/Bioinformatics/muscle3.8.31_i86darwin64"
cmd_line = MuscleCommandline(muscle_exe, input="HBA.all.fasta", out="HBA.aln", clw=" ")

print(cmd_line)

stdour, stderr = cmd_line()
#import subprocess
#muscle_result = subprocess.check_output([muscle_exe, "-in", "HBA.all.fasta", "-out", "HBA.aln"])
 
"/Users/jeonghyeonjeong/for github/Bioinformatics/muscle3.8.31_i86darwin64" -in HBA.all.fasta -out HBA.aln -clw
 

4. WebLogo로 보존 서열 알아보기

서열을 정렬하였으니 이제 시각화를 해 볼 차례입니다.

In [5]:
#5. WebLogo
from Bio import AlignIO #파일 읽는데 필요
from Bio.Alphabet import IUPAC #4.3.2에서 배운 Alphabet 모듈로, 서열 정보 하는데 필요
from Bio.motifs import Motif
from Bio import motifs
from Bio.Seq import Seq
 

Motif 란

단백질 서열 내의 특징이 보존된 짧은 펩타이드 서열이다. 이 motif는 단백질 내에서 구조적(ex: Zinc finger motif), 기능적으로 중요하다. 보통 짧은 것이 특징이나 특이적으로 15 ~ 20 정도의 긴 서열도 존재한다. 기능이 보존된 패턴이기에 Domain과 비교가 된다. motif 는 수학적 모델로 생성한 단순한 패턴인데 반해 Domain 은 단백질 구조(1차, 2차, 3차, 4차 구조)와 단백질 패밀리(family, superfamily 등) 에 연관되어 있고, 서열데이터 측면으로 보면 아미노산 길이의 차이가 있다.

출처: http://www.incodom.kr/Motif#h_1b91e03212546bb5d67ce0b82fa83af1

 

Bio.Motifs

Bio.motifs 패키지는 서열의 모티프를 분석할 때 사용된다. AlignAce 혹은 MEME 같은 프로그램을 통해 모티프 서열 분석을 진행할 수 있는데 Bio.motifs 패키지에서는 이들 프로그램의 결과를 쉽게 읽어 들일 수 있는 다양한 기능을 제공한다.

먼저 해당 패키지를 사용하기 위해서는 아래와 같이 단순히 호출만 해주면 된다. Biopython에 포함된 motif 패키지를 호출만 하면 되는 것이다.

from Bio import motifs 먼저 간단한 예를 만들어 보자. 먼저 motif 역시 서열 정보이기 때문에 from Bio.Seq import Seq를 이용하여야 한다. 그 후 각 서열 생성을 위해 Seq("서열") 기능을 통해 서열을 만들고 이를 하나의 인스턴스. 즉, 서열의 리스트로 만들어준다.

출처: http://www.incodom.kr/Biopython/Motifs

In [6]:
alignment = AlignIO.read("HBA.aln", "clustal")
instance = []

for record in alignment:
    s = Seq(str(record.seq), IUPAC.protein) #서열 생성을 위해 Seq("서열") 기능을 통해 서열을 만들고 이를 하나의 인스턴스. 즉, 서열의 리스트로 만들어준다.
    instance.append(s)
    
m = motifs.create(instance) #만들어진 인스턴스를 motifs.create() 메소드를 이용하여 아래와 같이 실행하면 모티프 오브젝트가 생성된다.

Motif.weblogo(m, 'HBA_WebLogo.png') #만들어진 모티프 오브젝트(m)의 Web Logo를 'HBA_WebLogo.png'이 파일로 저장해라

#Jupyter notebook을 사용한다면 다음 두 줄로 노트북 내에서 그림 파일을 확인 할 수 있습니다. 
from IPython.display import Image
Image("HBA_WebLogo.png")
    
Out[6]:
 

5. 계통수 그려보기

앞에서 Multiple Sequence Alignment를 시각화 하는 방법 중에서 다른 하나는 계통수를 그려보는 것으로, 생물 종 간 유전적 유사도를 표현하는 것입니다.

Bio.Phylo에 대한 더 많은 정보는 http://www.incodom.kr/Biopython/Phylo 를 참고할 수 있습니다.

이는 newick 포멧을 가지고 그려집니다. (앞서서 그린 WebLogo는 MUSCLE로 진행한 cluster 포맷 파일이기때문에 이를 newick 파일 포맷으로 변형해야 합니다. 이는 https://www.ebi.ac.uk/Tools/phylogeny/simple_phylogeny/ 을 통해 할 수 있습니다. )

In [7]:
#6. phylogenetic tree
from Bio import Phylo

tree = Phylo.read("HBA.newick", "newick")
print(tree)

Phylo.draw(tree)
 
Tree(rooted=False, weight=1.0)
    Clade()
        Clade(branch_length=0.00537)
            Clade(branch_length=0.29387, name='sp|P01994|HBA_CHICK')
            Clade(branch_length=0.01203)
                Clade(branch_length=0.01282)
                    Clade(branch_length=0.1181, name='sp|P60529|HBA_CANLF')
                    Clade(branch_length=0.10357, name='sp|P01948|HBA_RABIT')
                Clade(branch_length=0.00617)
                    Clade(branch_length=0.05673)
                        Clade(branch_length=0.0, name='sp|P69907|HBA_PANTR')
                        Clade(branch_length=0.0, name='sp|P69905|HBA_HUMAN')
                    Clade(branch_length=0.07405, name='sp|P01958|HBA_HORSE')
        Clade(branch_length=0.13313, name='sp|P18971|HBA_BALAC')
        Clade(branch_length=0.06808, name='sp|P01942|HBA_MOUSE')
 
<Figure size 640x480 with 1 Axes>