在没有输入文件的情况下在 Biopython 中创建对齐答案

【问题标题】：create alignment in Biopython without input file在没有输入文件的情况下在 Biopython 中创建对齐
【发布时间】：2023-03-17 07:57:01
【问题描述】：

我在字典中有一个蛋白质序列比对（id_prot 作为键，比对序列作为值；可能是另一种格式），我想用这个比对来用 Biopython 构建 NJ 树

但是，根据文档，加载用于系统发育分析的序列的唯一方法是从输入文件中加载。例如：

aln = AlignIO.read('Tests/TreeConstruction/msa.phy', 'phylip')

有谁知道如何在不读取输入文件的情况下在 aln 变量中加载序列？

【问题讨论】：

后退一两步 - 首先你的对齐方式是如何被读入字典的？
@Chris：我使用自定义 python 脚本创建了它。
您的自定义脚本可以输出标准对齐格式，那么您的问题就解决了
如下所述，我有成千上万的对齐方式......我想将所有内容都保存在 python 中，没有输出文件。否则我也可以使用标准的 NJ 程序而不是 Biopython。

标签： biopython

【解决方案1】：

从您的字典中创建 Phylip 格式的所需输入，并使用StringIO 加载它。请注意，序列/蛋白质 ID 最长可达 10 个字符。 IDs less than 10 characters must have spaces appended to them to reach the 10 character fixed width.

from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor
from Bio import AlignIO
from io import StringIO

dct = {'Alpha':   'AACGTGGCCACAT',
       'Beta':    'AAGGTCGCCACAC',
       'Gamma':   'CAGTTCGCCACAA',
       'Delta':   'GAGATTTCCGCCT',
       'Epsilon': 'GAGATCTCCGCCC'}


count = len(dct)
length = max(map(len, dct.values()))

msa = f" {count} {length}\n"
msa += '\n'.join(f"{prot_id:<10} {sequence}" for prot_id, sequence in dct.items())
print(msa)
print()

aln = AlignIO.read(StringIO(msa), 'phylip')
print(aln)
print()

calculator = DistanceCalculator('identity')
dm = calculator.get_distance(aln)
print(dm)
print()

constructor = DistanceTreeConstructor(calculator, 'nj')
tree = constructor.build_tree(aln)
print(tree)

输出：

 5 13
Alpha      AACGTGGCCACAT
Beta       AAGGTCGCCACAC
Gamma      CAGTTCGCCACAA
Delta      GAGATTTCCGCCT
Epsilon    GAGATCTCCGCCC

SingleLetterAlphabet() alignment with 5 rows and 13 columns
AACGTGGCCACAT Alpha
AAGGTCGCCACAC Beta
CAGTTCGCCACAA Gamma
GAGATTTCCGCCT Delta
GAGATCTCCGCCC Epsilon

Alpha   0
Beta    0.23076923076923073 0
Gamma   0.3846153846153846  0.23076923076923073 0
Delta   0.5384615384615384  0.5384615384615384  0.5384615384615384  0
Epsilon 0.6153846153846154  0.3846153846153846  0.46153846153846156 0.15384615384615385 0
    Alpha   Beta    Gamma   Delta   Epsilon

Tree(rooted=False)
    Clade(branch_length=0, name='Inner3')
        Clade(branch_length=0.18269230769230765, name='Alpha')
        Clade(branch_length=0.04807692307692307, name='Beta')
        Clade(branch_length=0.04807692307692307, name='Inner2')
            Clade(branch_length=0.27884615384615385, name='Inner1')
                Clade(branch_length=0.051282051282051266, name='Epsilon')
                Clade(branch_length=0.10256410256410259, name='Delta')
            Clade(branch_length=0.14423076923076922, name='Gamma')

【讨论】：

非常感谢您的回复。但是，由于我想做几十万次，我不想使用输入文件（以避免 I/O 限制）。