3.5 Dealing with alignments

특정 sequence 들을 align 할 수 있는 것은 종종 아주 유용하다. 나는 이것을 sequence 사이의 관계에 대한 "quick and dirty" 아이디어를 얻기 위해서 매우 자주 한다.

따라서, alignment를 해서 당신에게 다루기 쉬운 object를 돌려주는 python script를 빠르게 쓸 수 있는 것은 매우 멋지다. Biopython에서 alignment 와 관련된 code 란 python-level에서 alignment program 에 대한 access를 허용하는 것을 의미하며, 당신은 script 안에서 빠르게 alignment를 실행할 수 있다.

3.5.1 Clustalw

Clustalx (http://www-igbmc.u-strasbg.fr/BioInfo/ClustalX/Top.html) 는 multiple alignment를 하는데 매우 멋진 프로그램이다.

Biopython 은 Clustalx 에 의해 생성된 clustal format에서 alignment에 대한 access를 허용한다. (이것들은 보통 *.aln 확장자를 가진다.) 그것은 역시 clustalw 에 대한 access를 허용한다. 이것은 clustalx 의 command line version 이다.

Clustalw 와 interacting 하는 첫 step 은 program 에 pass 하고 싶은 command line을 set up 하는 것이다. Clustalw 은 많은 command line option을 가진다. 그리고 많은 parameter를 설정한다면, 당신은 많은 ol' command line을 typing 하다가 끝날 수 있다. 이 command line class 는 모든 옵션을 설정될 수 있는 class 의 attribute 로 만듬으로써 command line 은 모델링한다. 어떤 parameter를 설정하기 위한 몇 개의 편리한 기능이 역시 존재하여, parameter 에 대한 몇 개의 error checking이 이루어질 수 있다.

Clustalw multiple alignment를 하는 command line object를 만들기 위해서, 우리는 아래와 같이 한다.

import os
from Align.Clustalw import MultipleAlignCL

cline = MultipleAlignCL(os.path.join(os.curdir, 'opuntia.fasta'))
cline.set_output('test.aln')

(역자주 : 위 명령 그대로 하면 error 납니다. from Align.Clustalw import MultipleAlignCL 을 from Bio.Clustalw import MultipleAlignCL 으로 바꿔 줘야 합니다. )

처음에 우리는 MultipleAlignCL object를 import 한다. 이것은 clustalw 로부터 multiple alignment를 하는 것을 모델링한다. 우리는 그 다음에 우리가 alignment를 위해서 사용할 fasta file 의 single argument를 가지고 command line을 initialize 한다.

초기화 기능은 역시 옵션으로 clustalw 실행자의 위치를 정하는 두 번째 argument를 취한다. 디폴트로 command line 은 당신이 그것을 PATH 어딘가에 가지고 있다고 가정하고, 단지 'clustalw'를 불러낸다.

두 번째 argument 는 출력물을 test.aln 으로 가도록 설정한다. MultipleAlignCL object 는 역시 출력물의 format, gap cost 와 같은 것들을 설정하는 많은 parameter를 가지고 있다.

우리는 MultipleAlignCL class 의 str member attribute을 불러냄으로써 우리가 생성한 command line을 볼 수 있다. 이것은 str(cline)을 call 하거나 간단하게 print cline을 가지고 command line을 print 함으로써 이루어진다. 이 경우에, 이것은 다음의 output을 줄 것이다.

clustalw ./opuntia.fasta -OUTFILE=test.aln

(역자주 : 이것은 명령어가 아닙니다. 출력물입니다. 이것은 유닉스 쉘에서의 명령 그대로입니다. 하지만, 윈도우용 biopython 에서는 clustalw 와의 연동이 제대로 되지 않는 것 같습니다. Biopython 에서는 에러 나는 것이 print cline 했을때 나오는 것과 똑같은 것을 도스창에서 실행하면 잘 됩니다.)

• clustalx.exe를 fasta, py file과 같은 폴더에 옮겨주고 실행하면 윈도우에서도 이 예제를 실행할 수 있습니다.

이제, 우리는 간단한 command line을 set up 했다. 이제 우리는 그 commandline을 실행하여, 우리가 처리할 수 있는 결과를 얻기를 원한다. 이것은 Clustalw 의 do_alignment function을 다음과 같이 사용함으로써 이루어질 수 있다.

from Bio import Clustalw
alignment = Clustalw.do_alignment(cline)

당신이 이것을 실행시켰을 때 일어나는 것은 Biopython 이 당신의 command line과 clustalw를 주어진 parameter를 가지고 실행시키는 것이다.

(역자주 : What happens when you run this if that Biopython executes your command line and runs clustalw with the given parameters. 이 문장 중간의 if 가 is 의 오기가 아닐까 의심됨)

Biopython 이 parse 할 수 있는 format 이라면, (현재는 단지 clustal format 뿐이다.) 그것은 그 다음에 output을 grab 하고, 그 결과를 parse 해서 그것들을 적절한 형태의 alignment object 로 돌려준다. 그래서 이 경우에, 우리는 clustal format 의 default 로 결과를 얻었기 때문에, 반환된 alignment object 는 ClustalAlignment type 일 것이다.

우리가 한번 이 alignment를 얻으면, 우리는 alignment 에 관련된 모든 sequence 에 대한 seq_record object를 얻는 것과 같은 흥미로운 것들을 할 수 있다.

all_records = alignment.get_all_seqs()

print 'description:', all_records[0].description
print 'sequence:', all_records[0].seq

이것은 alignment에서 첫 번째 seuqence에 대한 description 과 sequence object를 print out 할 것이다.

description: gi|6273285|gb|AF191659.1|AF191
sequence: Seq('TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAAAAAAATGAAT 
...', IUPACAmbiguousDNA())

당신은 역시 alignment 의 maximum length를 계산할 수 있다.

length = alignment.get_alignment_length()

최종적으로 alignment object를 original format 으로 쓰기 위해서 우리는 단지 str function을 access 하기만 하면 된다. Alingment를 print 하면 다음과 같을 것이다.

CLUSTAL X (1.81) multiple sequence alignment

gi|6273285|gb|AF191659.1|AF191      TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA
gi|6273284|gb|AF191658.1|AF191      TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAA
...

이것은 모든 original 정보를 손상하지 않고 당신의 alignment를 file 로 쓰는 것을 쉽게 만든다.

만약 당신이 alignment를 가지고 보다 흥미로운 것들을 하기를 원한다면, 가장 멋진 것은 alignment를 section 3.5.2 에 기술된 SummaryInfo object 와 같은 alignment information generating object 로 넘기는 것이다.

3.5.2 summary information 계산

당신이 일단 alignment를 가지면, 그에 관한 정보를 찾기를 바랄 것이 틀림없다. Alignment object 그 자체에 있는 alignment 에 대한 정보를 생성할 수 모든 funtion을 갖기 보다는, 우리는 기능을 alignment 에 작동하는 독립적인 class 로 분리하고자 시도하였다.

한 object 에 대한 summary 정보를 계산하기 위해 준비하는 것은 쉽다. alignment 라 불리는 alignment object를 가지고 있다고 하자. 우리가 해야할 것은 summary information을 계산할 object를 얻어야 한다.

from Bio.Align import AlignInfo
summary_align = AlignInfo.SummaryInfo(alignment)

Summary_align object 는 매우 유용하고, 다음과 같은 신선한 것들을 당신을 위해서 해 줄 것이다.

Quick consensus sequence 의 계산 -- see section 3.5.3 Alignment 에 대한 position specific score matix 얻기 -- see section 3.5.4 Alignment 에 대한 information content 의 계산 -- see section 3.5.5 Alingment 에서의 substitution 에 대한 정보의 생성 -- 이것을 사용에서 substitution matrix를 생성하는 것을 section 3.6 이 상술한다.

3.5.3 Quick consensus sequence 의 계산

Section 3.5.2 에 기술된 SummaryInfo object 는, alignment 의 quick consensus를 계산하는 기능을 제공한다. 우리가 summary_align 이라 불리는 SummaryInfo object를 가졌다고 가정하면, 우리는 consensus를 다음과 같이 함으로써 계산할 수 있다.

consensus = summary_align.dumb_consensus()

이름이 제시하는 것과 같이 이것은 정말로 간단한 consensus 계산기로서, 단지 consensus 의 각 point에서 모든 residue 를 더할 뿐이다.

만약 가장 공통된 value 가 어떤 threshold value 보다 높다면 (default 는 .3 이다.) consensus 에 common residue를 더할 것이다. 만약 그것이 threshold 에 이르지 않는다면, 그것은 consensus 에 ambiguity character를 더할 것이다. 반환된 consensus object 는 Seq object 이고, 그 alphabet 은 consensus를 구성한 sequence 들의 alphabet 으로부터 추론된 것이다. 그러므로 consensus 에 대해 print 하면 다음과 같다.

consensus Seq('TATACATNAAAGNAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAAAAAAATGAAT 
...', IUPACAmbiguousDNA())

당신은 optional parameter를 줌으로써 어떻게 dumb_consensus 가 작동하는 지를 조절할 수 있다.

threshold 한 특정 residue 가 얼마나 공통되게 그것이 add 되기 전에 position 에 있어야 하는지를 정하는 threshold 이다. Default 는 .7 이다.

ambiguous character 이것은 사용할 ambiguity character 이다. Default 는 'N'이다.

consensus alphabet 이것은 consensus sequence의 사용을 위한 alphabet 이다. 만약 어떤 alphabet 이 정해지지 않았다면, 우리는 alignment 에 있는 sequence 의 alphabet 에 근거한 alphabet를 추측하려 할 것이다.

(역자주 : If an alphabet is not specified than we will try to guess the alphabet based on the alphabets of the sequences in the alignment. 이 문장에서 'than'이 'then' 으로 바뀌어야 되지 않을까 싶습니다.)

3.5.4 Position Specific Score Matrices

Position specific score matrices (PSSMs) 는 consensus 와는 다른 방법으로 alignment information을 요약하고, 다른 작업들을 위해서 유용할 수 있다. 기본적으로 PSSM 이 count matrix 이다. Alignment 각각의 column 에 대해, 각각의 alphabet 문자가 세어져서 더해진다. 왼쪽 축을 따른 어떤 대표적인 sequence 에 대해서 total 이 표시된다. 이 sequence 는 consensus sequence 일 수도 있지만, alignment 의 어떤 sequence 일수도 있다. 예를 들어, 다음의 alignment 일 때,

GTATC
AT--C
CTGTC

이 alignement 에 대한 PSSM 은 다음과 같다.

      G A T C
    G 1 1 0 1
    T 0 0 3 0
    A 1 1 0 0
    T 0 0 2 0
    C 0 0 0 3

우리가 c_align 이라는 alignment object를 가지고 있다고 가정하자. Side를 따라 consensus sequence 의 PSSM을 얻기 위하여, 우리는 처음에 summary object를 얻어서 consensus sequence를 계산한다.

summary_align = AlignInfo.SummaryInfo(c_align)
consensus = summary_align.dumb_consensus()

이제, 우리는 이것을 계산할 때 어떤 N 의 애매한 residue 도 무시한 PSSM을 만들기를 원한다.

my_pssm = summary_align.pos_specific_score_matrix(consensus, chars_to_ignore = ['N'])

이에 관해서 두가지를 주목해야 한다.

알파벳에서 엄격함을 유지하기 위해서 당신은 alignment object 의 알파벳에 있는 PSSM 의 top을 따라 character만 포함시킬 수 있다. Gap 은 PSSM의 top axis를 따라 포함되지 않는다.

Axis 의 왼쪽 side를 따라 표시되도록 넘겨진 sequence 는 consensus 일 필요가 없다. 예를 들어, 당신이 이 axis를 따라 alignment 에서의 두 번째 sequence를 표시하기를 원한다면, 이렇게 할 필요가 있다.

second_seq = alignment.get_seq_by_num(1)
my_pssm = summary_align.pos_specific_score_matrix(second_seq, chars_to_ignore = ['N'])

위의 command 는 PSSM object를 반환한다. 위에서 보여진 바와 같이 PSSM을 출력하기 위해서 우리는 단지 my_pssm을 출력하기만 하면 된다. 이것은 다음과 같은 결과를 보여준다.

    A   C   G   T
T  0.0 0.0 0.0 7.0
A  7.0 0.0 0.0 0.0
T  0.0 0.0 0.0 7.0
A  7.0 0.0 0.0 0.0
C  0.0 7.0 0.0 0.0
A  7.0 0.0 0.0 0.0
T  0.0 0.0 0.0 7.0
T  1.0 0.0 0.0 6.0
...

당신은 your_pssm[sequence_number][residue_count_name] 와 같이 기술함으로써 PSSM의 어떤 element 도 access 할 수 있다. 예를 들어, 위의 PSSM 의 second element 에 있는 'A' residue 에 대한 count를 얻기 위해서, 다음과 같이 한다.

>>> print my_pssm[1]['A']
7.0

PSSM class 의 구조는 다행스럽게도 element를 access 하는 것과 matrix를 예쁘게 프린트하기에 쉽다.

3.5.5 Information Content

진화적 conservation 의 가능한 유용한 measure 는 sequence 의 information content 이다.

분자생물학자에게 정보 이론에 대한 유용한 소개는 http://www.lecb.ncifcrf.gov/~toms/paper/primer/ 에서 볼 수 있다. 우리의 목적을 위해서, 우리는 consensus sequence 의 information content나, consensus sequence 의 부분을 볼 것이다. 우리는 아래의 공식을 사용하여 multiple sequence alignment 의 특정 column 에서 information content를 계산한다.

 공식 생략 

ICj -- alignment 의 jth jth column 에서의 information content Na -- alphabet 에서의 문자의 수 Pij -- column 에서의 특정한 문자의 빈도 (만약 aligment column에서 G 가 6번 중에 3번 나오면 이것은 0.5이다.) Qi -- 문자의 expected frequency. 이것은 optional argument 이다. 이것의 사용은 사용자의 판단에 맡겨진다. Default 로, 그것은자동적으로 protein alphabet을 위해서는 0.05로, nucleic acid alphabet를 위해서는 0.25로 설정된다. 이것은 prior distribution 에 대한 어떤 가정도 없는 information content를 얻기 위한 것이다. Prior를 가정할 때나, non-standard alphabet을 사용할 때, 사용자는 Qi 의 수치를 제공해야 한다.

이제, 우리는 어떻게 Biopython에서 information content 가 계산되는지에 대해 알았다. Alignment 의 특정 region에서 어떻게 그것을 얻는지에 대해서 보자.

처음에, 우리는 alignment summary object를 얻기 위해서 우리의 alignment를 사용할 필요가 있다. 이것을 우리는 summary_align 으로 가정할 것이다. (이것을 어떻게 얻는지에 대해서는 section 3.5.2를 보라) 한번 우리가 이 object를 얻은 다음에는, 한 region 에 대해서 information content를 계산하는 것은 다음과 같이 쉽다.

info_content = summary_align.information_content(5, 30, chars_to_ignore = ['N'])

와! 저것은 위의 공식보다 훨씬 쉽다. info_content 라는 변수는 이제 특정 region 에 대한 (5개에서 30개의 alignment 로부터 나온) information content인 float value를 가지고 있다. 이 수치는 우리의 alphabet 에 포함되지 않았기 때문에, inforamtion content를 계산할 때, 우리는 특히 ambiguity residue 인 'N'을 무시한다. (따라서 우리는 그것을 보는데 관심을 가질 필요가 없다.)

위에서 언급된 것과 같이, expected frequencies를 제공함으로써, 우리는 relative information content를 역시 계산할 수 있다.

expect_freq = {
    'A' : .3,
    'G' : .2,
    'T' : .3,
    'C' : .2}

기대값은 raw dictionary 로 pass 되어서는 안되고, 대신에, SubsMat.FreqTable object 로 pass 되어야 한다. (FreqTables 에 대한 더 많은 정보는 section 4.4.2를 보라.) FreqTables object 는 Biopython Seq class 가 작동하는 방법과 유사하게, dictionary 와 Alphabet을 연관짓는 기준을 제공한다.

Frequency dictionary 로부터 FreqTable object를 생성하기 위해서, 우리는 다음과 같이 할 필요가 있다.

from Bio.Alphabet import IUPAC
from Bio.SubsMat import FreqTable

e_freq_table = FreqTable.FreqTable(expect_freq, FreqTable.FREQ, IUPAC.unambigous_dna)

이제 우리는 그것을 얻었다. 우리의 alignment 의 region 에 대한 relative information content를 계산하는 것은 다음과 같이 쉽다.

info_content = summary_align.information_content(5, 30, e_freq_table = e_freq_table, chars_to_ignore = ['N'])

이제, info_content 는 expected frequencies에 관련하여 region 에 대한 relative information content를 가진다.

반환된 값은 위의 공식에서 base 2 인 logarithm을 사용해서 계산되어진다. 당신은 이것을 log_base parameter를 당신이 원하는 base 로 넘김으로써 수정할 수 있다.

info_content = summary_align.information_content(5, 30, log_base = 10, chars_to_ignore = ['N'])

이제, 당신은 information content를 계산할 준비가 되었다. 당신이 이것은 어떤 real life problem 에 적용하기를 원한다면, 아마도 그것이 어떻게 사용되는지에 대한 아이디어를 얻기 위해 information content 에 대한 문헌을 파는 것이 가장 좋을 것이다. 당신의 파는 것이 이 function 의 coding 에서의 어떤 실수를 밝혀내지 않기를 바란다.

3.5.6 Alignment format 사이의 Translating

항상 당신이 해야만 하는 것으로 끝나는 하나의 것은 서로 다른 format 사이의 변환이다. Biopython 은 이것을 alignment object를 위한 FormatConverter class를 사용해서 한다. 처음에, 우리가 단지 clustal format 의 alignment를 ClustalAlignment object 로 parsing 했다고 하자.

import os
from Bio import Clustalw

alignment = Clustalw.parse_file(os.path.join(os.curdir, 'test.aln'))

이제, 이 alignment를 FASTA format 으로 변환하지, 처음에 우리는 converter object를 생성한다.

from Bio.Align.FormatConvert import FormatConverter
converter = FormatConverter(alignment)

우리는 converter 에 우리가 convert 하고자 하는 alignment를 pass 한다. 이제, 이것을 FASTA alignment format 으로 얻기 위해서, 우리는 단지 다음과 같이 한다.

fasta_align = converter.to_fasta()

새로이 생성된 fasta_align object를 fasta_align을 프린트하여 보면 다음과 같다.

>gi|6273285|gb|AF191659.1|AF191
TATACATTAAAGAAGGGGGATGCGGATAAATGGAAAGGCGAAAGAAAGAATATATA----
------ATATATTTCAAATTTCCTTATATACCCAAATATAAAAATATCTAATAAATTAGA
...

변환 과정은 물론 alignment format 특이적인 정보를 잃는다. 그러나, 대부분의 기본적인 alignment 에 대한 정보는 유지될 것이다.

더 많은 format 이 더해짐에 따라, converter 는 모든 이러한 서로 다른 format을 읽고 쓰도록 강화될 것이다.