바이오파이썬 서열정렬 예제

바이오파이썬 튜토리얼 6장 내용입니다. 

 

 

from Bio.Seq import Seq

from Bio.SeqRecord import SeqRecord
from Bio.Align import MultipleSeqAlignment
alignment = MultipleSeqAlignment(
[
SeqRecord(Seq("ACTCCTA"), id='seq1'),
SeqRecord(Seq("AAT-CTA"), id='seq2'),
SeqRecord(Seq("CCTACT-"), id='seq3'),
SeqRecord(Seq("TCTCCTC"), id='seq4'),
]
)

print(alignment)

 

Alignment with 4 rows and 7 columns
ACTCCTA seq1
AAT-CTA seq2
CCTACT- seq3
TCTCCTC seq4

substitutions = alignment.substitutions
print(substitutions)

 

    A    C    T
A 2.0  4.5  1.0
C 4.5 10.0  0.5
T 1.0  0.5 12.0

서열 정렬의 쌍은 불규칙적이기 때문에, 개수는 대각선의 위와 아래에서 동일하게 나눠집니다. 예를 들어, A와 C의 서열은 [‘A’, ‘C’]에서 4.5, [‘C’, ‘A’]에서 4.5가 저장되어 있습니다. 이러한 계산 방식은 이 결과들을 이용해 20.4.2에 나타난 것처럼 대체 행렬을 계산할 때 더욱 편리하게 해줍니다. alignment.substitutions는 정렬에 나타나는 문자만을 포함하고 있습니다. 정렬에 나타나지 않는 문자의 성분을 추가하려면 select 메서드를 사용할 수 있습니다.

 

m = substitutions.select("ATCG")
print(m)

 

    A    T    C   G
A 2.0  1.0  4.5 0.0
T 1.0 12.0  0.5 0.0
C 4.5  0.5 10.0 0.0
G 0.0  0.0  0.0 0.0

Clustal W는 다중 서열 정렬의 인기있는 command line tool입니다 (Clustal X라는 그래픽 인터페이스도 존재합니다). 바이오파이썬의 Bio.Allign.Applications 모듈은 이 정렬 도구(그리고 다른 것들의)의 wrapper를 가지고 있습니다. 파이썬에서 ClustalW를 사용하기 전에, 스스로 ClustalW를 command line에서 직접 사용해 다른 옵션에 익숙하게 만드는 것을 먼저 시도해 보기 바랍니다. 그러면 바이오파이썬 wrapper가 실제 명령줄 API에 충실하다는 것을 발견할 수 있을겁니다

 

from Bio.Align.Applications import ClustalwCommandline
help(ClustalwCommandline)

 

  http://www.clustal.org/

 

opuntia.fasta

0.01MB

일곱 가시부채 선인장(부채선인장 속)의 DNA 염기서열인 opuntina.fasta입니다. 

 

from Bio.Align.Applications import ClustalwCommandline
cline = ClustalwCommandline("clustalw2", infile="opuntia.fasta")
print(cline)

 

 

http://www.clustal.org/download/current/

Index of /download/current

위 링크에서 운영체제 버전에 맞는 파일을 다운도르 받으면 됩니다. 

 

import os
from Bio.Align.Applications import ClustalwCommandline
clustalw_exe = r"C:\Program Files (x86)\clustalw2\clustalw2.exe"
clustalw_cline = ClustalwCommandline(clustalw_exe, infile="opuntia.fasta")
assert os.path.isfile(clustalw_exe), "Clustal W executable missing"
stdout, stderr = clustalw_cline()

 

에러가 나타나면 

C:\Program Files (x86)\clustalw2\clustalw2.exe

 

설치 위치가 잘못된 경우입니다. 본인 컴퓨터 어느 위치에 설치가 되어 있는지 찾아봐야 합니다. 

에러가 없으면 자동으로 opuntia.aln 파일이 저장됩니다. 

 

from Bio import AlignIO
align = AlignIO.read("opuntia.aln", "clustal")
print(align)

 

 

Alignment with 7 rows and 906 columns
TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273285|gb|AF191659.1|AF191
TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273284|gb|AF191658.1|AF191
TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273287|gb|AF191661.1|AF191
TATACATAAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273286|gb|AF191660.1|AF191
TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273290|gb|AF191664.1|AF191
TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273289|gb|AF191663.1|AF191
TATACATTAAAGGAGGGGGATGCGGATAAATGGAAAGGCGAAAG...AGA gi|6273291|gb|AF191665.1|AF191

 

계통수를 자동으로 나타냅니다. 뉴윅버전으로 나타내는데 간단하게 실행해 봅시다. 

 

from Bio import Phylo
tree = Phylo.read("opuntia.dnd", "newick")
Phylo.draw_ascii(tree)

 

                             _______________ gi|6273291|gb|AF191665.1|AF191665
  __________________________|
 |                          |   ______ gi|6273290|gb|AF191664.1|AF191664
 |                          |__|
 |                             |_____ gi|6273289|gb|AF191663.1|AF191663
 |
_|_________________ gi|6273287|gb|AF191661.1|AF191661
 |
 |__________ gi|6273286|gb|AF191660.1|AF191660
 |
 |    __ gi|6273285|gb|AF191659.1|AF191659
 |___|
     | gi|6273284|gb|AF191658.1|AF191658

 

MUSCLE은 Clustal W보다 더 최근에 만들어진 다중 서열 정렬 도구이고, 바이오파이썬은 이를 위한 Wrapper를 Bio.Align.Applications 모듈에서 지원합니다. 

 

 

from Bio.Align.Applications import MuscleCommandline
help(MuscleCommandline)

 

http://www.drive5.com/muscle/

 

from Bio.Align.Applications import MuscleCommandline
cline = MuscleCommandline(input="opuntia.fasta", out="opuntia.txt")
print(cline)

muscle -in opuntia.fasta -out opuntia.txt

 

https://github.com/rcedgar/muscle/releases/tag/v5.1

Release Version 5.1 · rcedgar/muscle

 

다운로드하여서 위에 clustal과 같은 방법으로 실행해 보면 됩니다. 

 

 

EMBOSS

 

https://www.ebi.ac.uk/Tools/emboss/

EMBOSS programs < EMBL-EBI

 

Smith-Waterman 알고리즘 로컬 서열과 Needleman Wunsch global alignment를 위해 water와 needle로 접근해서 서열분석을 할 수 있습니다.

 

 

beta.faa

0.00MB

alpha.faa

0.00MB

 

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

 

 

from Bio.Emboss.Applications import NeedleCommandline
help(NeedleCommandline)

 

 

 

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

 

 

정렬

 

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)

 

len(alignments)

 

80

print(alignments[0])

 

Alignment(seqA='MV-LSPADKTNV---K-A--A-WGKVGAHAG---EY-GA-EALE-RMFLSF----PTTK-TY--FPHFDL-SH-G--S---AQVK-G------HGKKV--A--DA-LTNAVAHVD-DMPNALS----A-LSD-LHAH--KLR-VDPV-NFKL-LSHCL---LVT--LAAHLPA----EFTPA-VH-ASLDKFLAS---VSTV------LTS--KYR-', seqB='MVHL-----T--PEEKSAVTALWGKV-----NVDE-VG-GEAL-GR--L--LVVYP---WT-QRF--F--ES-FGDLSTPDA-V-MGNPKVKAHGKKVLGAFSD-GL----AH--LD--N-L-KGTFATLS-EL--HCDKL-HVDP-ENF--RL---LGNVLV-CVL-AH---HFGKEFTP-PV-QA------A-YQKV--VAGVANAL--AHKY-H', score=72.0, start=0, end=217)

 

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

 

 

 

https://tsyoon.tistory.com/37?category=988004 

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)
len(alignments)

 

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

 

MV-LSPADKTNVKAAWGKVGAHAGEYGAEALERMFLSFPTTKTYFPHF-DLS-----HGSAQVKGHGKKVADALTNAVAHVDDMPNALSALSDLHAHKLRVDPVNFKLLSHCLLVTLAAHLPAEFTPAVHASLDKFLASVSTVLTSKYR
|| |.|..|..|.|.||||  ...|.|.|||.|.....|.|...|..| |||     .|...||.|||||..|.....||.|........||.||..||.|||.||.||...|...||.|...||||.|.|...|..|.|...|..||.
MVHLTPEEKSAVTALWGKV--NVDEVGGEALGRLLVVYPWTQRFFESFGDLSTPDAVMGNPKVKAHGKKVLGAFSDGLAHLDNLKGTFATLSELHCDKLHVDPENFRLLGNVLVCVLAHHFGKEFTPPVQAAYQKVVAGVANALAHKYH
  Score=292.5

 

서열 비교의 간단한 예제입니다. 

 

 

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 sequence 명령어를 넣으면 전체 서열 비교에 대한 내용이 나타납니다. 

 

 

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], full_sequences=True))

 

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

4번째 3번째부터 시작되는 예제입니다.

 

from Bio import pairwise2
from Bio.Align import substitution_matrices
blosum62 = substitution_matrices.load("BLOSUM62")
alignments = pairwise2.align.localds("LSSPADKTNVKKAA", "DDPEEKSAVNN",blosum62, -10, -1)
print(pairwise2.format_alignment(*alignments[0]))

 

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

 

서열 정렬에 대한 내용은 다음 링크에 잘 나와 있습니다. 

 

https://tsyoon.tistory.com/32?category=988004 

서열분석

 

 

biopython_chapter6.ipynb

0.05MB