上海锐翌生物科技有限公司

服务电话:021-51001612

邮箱:support@realbio.cn

技术课堂

Biopython教程之“多序列比对”
发布时间:2017-12-25 11:31   点击率:
多序列比对(Multiple Sequence Alignment, MSA),对多个序列进行对位排列。这通常需要保证序列间的等同位点处在同一列上,并通过引进小横线(-)以保证最终的序列具有相同的长度。

在生物信息分析中,我们有时需要进行多序列比对,Biopython可以帮我们实现,特别是使用linux系统的同学,biopython值得拥有。


 
两种读取方法

Biopython提供了两种方法读取多序列比对数据,即Biopython提供的AlignIO.read和AlignIO.parse模块。这两种方法跟SeqIO处理一个和多个数据的设计方式是一样的。AlignIO.read() 只能读取一个多序列比对,而AlignIO.parse()可以依次读取多个序列比对数据。
 

使用AlignIO.parse()将会返回一个多序列比对的迭代器(iterator)。迭代器往往在循环中使用,即使用for或while都可以逐行读取,例如:
For line in iterator:(iterator为迭代器的变量名)

在实际数据分析过程中,会时常处理包含有多个多序列比对的文件。例如PHYLIP中的seqboot,Bill Pearson的FASTA程序。如果你所遇到的文件仅仅包括一个多序列比对,你应该使用AlignIO.read(),这将返回一个MultipleSeqAlignment对象。

使用这两个函数都必须注明两个参数:
第一个参数为包含有多序列比对数据的句柄(handle)。在实际操作中,这往往是一个打开的文件()或者文件名。

第二个参数为多序列比对文件格式(小写)。与Bio.SeqIO模块一样,Biopython不会对将读取的文件格式进行猜测。Bio.AlignIO支持多种比对格式,包括所有Bio.AlignIO模块支持的多序列比对数据格式可以在 http://biopython.org/wiki/AlignIO中找到。

Bio.AlignIO 模块还接受一个可选参数seq_count 。它可以处理不确定的多序列比对格式,或者包含有多个序列的排列。另一个可选参数 alphabet 允许用户指定序列比对文件的字符(alphabet),它可以用来说明序列比对的类型(DNA,RNA或蛋白质)。因为大多数序列比对格式并不区别序列的类型,因此指定这一参数可能会对后期的分析产生帮助。

 

 
多个序列比对读取

下面,我们可以开始进行序列比对的数据读取。假设我们有一个PHYLIP格式的很小的序列比对:
5    6
tests     AACAAC
greep    AACCCC
Gamlp   ACCAAC
Delta    CCACCA
Epon    CCAAAC

如果你想用PHYLIP工具包来bootstrap一个系统发生树,其中的一个步骤是用 bootseq 程序来产生许多序列比对。将给出类似于以下格式的序列比对:
5     7
Test      AAACCAA
Bet       AAACCCA
Semma    ACCCCAC
Delta     CCCAACC
Epsilon    CCCAAAC
4     5
Klpha     AAACAA
Yeta      AAACCC
Dalta     CCCACC
Silon     CCCAAA
5     5
Alpha     AAAAA
Beta      AACCC
Tema     ACAAC
Delta     CCCCA
Epsilon   CCAAC
5     6
Reha      AAAACC
Beta      ACCCCC
Gamma    AAAACC
KLOlta    CCCCAA
Epsilon    CAAACC

如果你想用AlignIO来读取这个文件,你可以使用:
Import Bio
align = Bio.AlignIO.parse("resampled.phy", "phylip")
for line in align:
    print line
第一行是导入Biopython模块;
第二行是读取resampled.phy文件,使用phylip格式;
第三行是循环读取迭代器内容;
第四行就是打印读取的内容。

这将给出以下的输出(这时只显示缩略的一部分):
SingleLetterAlphabet() alignment with 5 rows and 7 columns
AAACCAA Test
AAACCCA Bet
ACCCCAC Semma
CCCAACC Delta
CCCAAAC Epsilon

SingleLetterAlphabet() alignment with 4 rows and 5 columns
AACAA Klpha
AACCC Beta
CCACC Dalta
CCAAA Silon

SingleLetterAlphabet() alignment with 5 rows and 5 columns
AAAAC Alpha
AACCC Yeta
ACAAC Gamma
CCCCA Delta
CCAAC Epsilon

SingleLetterAlphabet() alignment with 5 rows and 6 columns
AAAACC Qaha
ACCCCC Beta
AAAACC Telma
CCCCAA Delta
CAAACC Silon

与Bio.SeqIO.parse一样,Bio.AlignIO.parse()将返回一个迭代器(iterator)。如果你希望把所有的序列比对都读取到内存中,我们可以把它储存在数列里。
from Bio import AlignIO
align = list(AlignIO.parse("resampled.phy", "phylip"))
last_align = align[-1]
first_align = align[0]
第一行为从Biopython导入AlignIO模块;
第二行是将比对存储在变量为align的数列里;
第三行是读取数列的最后一个元素;
第四行是读取数列的第一个元素;
当然,你也可以任意获取你想要的任何元素,只要修改元素下标就好了。


 
序列比对的输出

前面我们讨论了Bio.AlignIO.read()和Bio.AlignIO.parse()来读取各种格式的序列比对,现在让我们来使用Bio.AlignIO.write()输出序列比对文件。

使用这个函数都必须注明三个参数:
1. 第一个参数为一个 MultipleSeqAlignment 对象(或者是一个 Alignment 对象),可以通过循环迭代已经生成的迭代器获得;
2. 第二个参数为多序列比对文件格式(小写);
3. 第三个参数为期望使用的文件名。
我们来看个例子:
from Bio.Align import MultipleSeqAlignment
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import generic_dna
from Bio.Seq import Seq
one = MultipleSeqAlignment([
             SeqRecord(Seq("ACTGCTAGCTAG", generic_dna), id="A"),
             SeqRecord(Seq("ACT-CTAGCTAG", generic_dna), id="B"),
             SeqRecord(Seq("ACTGCTAGDTAG", generic_dna), id="C"),
         ])
two = MultipleSeqAlignment([
             SeqRecord(Seq("GTCAGC-AG", generic_dna), id="D"),
             SeqRecord(Seq("GACAGCTAG", generic_dna), id="E"),
             SeqRecord(Seq("GTCAGCTAG", generic_dna), id="F"),
         ])

three = MultipleSeqAlignment([
             SeqRecord(Seq("ACTAGTACAGCTG", generic_dna), id="G"),
             SeqRecord(Seq("ACTAGTACAGCT-", generic_dna), id="H"),
             SeqRecord(Seq("-CTACTACAGGTG", generic_dna), id="I"),
         ])
align = [one,two,three]
这里我们不再解释每一行的用法,感兴趣的同学可以对python作一定了解,就可以明白。

假设我们有一个包含三个 MultipleSeqAlignment 对象的列表align,现在我们将它写出为PHYLIP格式:
from Bio import AlignIO
AlignIO.write(align, "test.phy", "phylip")

打开test.phy文件,你会看到以下内容:
3   12
A      ACTGCTAGCT AG
B      ACT-CTAGCT AG
C      ACTGCTAGDT AG
3   9
D      GTCAGC-AGs
E      GACAGCTAG
F      GTCAGCTAG
3   13
G      ACTAGTACAG CTG
H      ACTAGTACAG CT-
I       -CTACTACAG GTG

这样,我们通过Biopython简单地进行了多序列比对的读取和写入操作。Biopython的多序列比对的其他用法,感兴趣的童鞋可以去官网进一步了解。

 

锐翌原创文章,未经授权严禁转载。