多序列比对(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的多序列比对的其他用法,感兴趣的童鞋可以去官网进一步了解。
http://biopython.org/
撰稿:李心如
编辑:王丽燕
转自:锐翌基因
- 本文固定链接: https://maimengkong.com/kyjc/1562.html
- 转载请注明: : 萌小白 2023年6月11日 于 卖萌控的博客 发表
- 百度已收录