Biopython给出了ValueError：即使序列的长度相同，序列的长度也必须相同

Question

我试图通过从我的数据生成一个.phy文件来创建一个系统发育树

我有一个数据帧

ndf=

ESV trunc

1 esv1 TACGTAGGTG...

2 esv2 TACGGAGGGT...

3 esv3 TACGGGGGG...

7 esv7 TACGTAGGGT...

我检查了列“trunc”元素的长度：

length_checker = np.vectorize(len)

arr_len = length_checker(ndf['trunc'])

结果arr_len为所有元素提供了相同的长度（=253）

我将此数据帧保存为.phy文件，如下所示：

23 253

1. esv1 TACGTAGGTGGCGAGCGTTATCCGGAATTATTGGGCGTAAAGCGCGCGCAGGCGGTTTCTTAAGTCTGATGTGAAAGCCCACGGCTCAACCGTGGAGGGTCATTGGAAACTGGGAAACTTGAGTGCAGAAGAGGAAAGCGGAATTCCACGTGTAGCGGTGAAATGCGTAGAGATGTGGAGGAACACCAGTGGCGAAGGCGGCTTTCTGGTCTGTAACTGACGCTGAGGCGCGAAAGCGTGGGGAGCAAACAGG

2. esv2 TACGGAGGGTGCAAGCGTTATCCGGATTCACTGGGTTTAAAGGGTGCGTAGGTGGGTTGGTAAGTCAGTGGTGAAATCTCCGAGCTTAACTTGGAAACTGCCATTGATACTATTAATCTTGAATATTGTGGAGGTTAGCGGAATATGTCATGTAGCGGTGAAATGCTTAGAGATGACATAGAACACCAATTGCGAAGGCAGCTGGCTACACATATATTGACACTGAGGCACGAAAGCGTGGGGATCAAACAGG

3. esv3 TACGGGGGGGGCAAGCGTTGTTCGGAATTACTGGGCGTAAAGGGCGCGTAGGCGGCCAGACCAAGTCGAGTGTGAAATTGCAGGGCTTAACTTTGCAGGGTCGCTCGATACTGGTCGGCTAGAGTGTGGAAGAGGGTACTGGAATTCCCGGTGTAGCGGTGAAATGCGTAGATATCGGGAGGAACACCAGCGGCGAAGGCGGGTACCTGGGCCAACACTGACGCTGAGGCGCGAAAGCTAGGGGAGCAAACAG

这与this tutorial中使用的文件类似

但是，当我运行命令时 aln = AlignIO.read('msa.phy', 'phylip')

我得到“ValueError：序列必须都是相同长度的”

我不知道为什么会这样，也不知道如何修复它。非常感谢您的帮助

谢谢

Answer 1

一般来说，phylip是不同程序之间系统发育最复杂的格式。有严格的phylip格式和放松的phylip格式等。。。t不容易知道所使用的分隔符、空格字符和/或回车符

我认为您似乎在分类单元的名称（即序列标签）和序列名称之间留了一个空格，即

2. esv2

Phylip格式正在监视标签和序列数据之间的空间。在本例中，序列长度为3bp。使用“.”通常也不是一个好主意。整数似乎不表示行号

另一个问题是，您可以/应该尝试将序列与标签保持在同一行，并删除回车符，即

esv2 TACGGAGGGTGCAAGCGTTATCCGGATTCACTGGGTTTAAAGGGTGCGTAGGTGGGTTGGTAAGTCAGTGGTGAAATCTCCGAGCTTAACTTGGAAACTGCCATTGATACTATTAATCTTGAATATTGTGGAGGTTAGCGGAATATGTCATGTAGCGGTGAAATGCTTAGAGATGACATAGAACACCAATTGCGAAGGCAGCTGGCTACACATATATTGACACTGAGGCACGAAAGCGTGGGGATCAAACAGG

有时回车符不起作用（这可能是phylip格式），传统格式使用空格字符“”。我总是保持统一数量的空格来保持对齐。。。不确定是否需要这样做

注意如果分类单元名称超过10个字符，则需要宽松的phylip格式，这种格式在任何情况下都是一个好主意

最终的解决方案是，所有其他失败的都是转换为fasta，作为fasta导入，然后转换为phylip。如果这一切都失败了。。。回来后还有更多的问题需要解决

---

Fasta格式删除“23 254”标题，然后每个序列如下所示

>esv2
TACGGAGGGTGCAAGCGTTATCCGGATTCACTGGGTTTAAAGGGTGCGTAGGTGGGTTGGTAAGTCAGTGGTGAAATCTCCGAGCTTAACTTGGAAACTGCCATTGATACTATTAATCTTGAATATTGTGGAGGTTAGCGGAATATGTCATGTAGCGGTGAAATGCTTAGAGATGACATAGAACACCAATTGCGAAGGCAGCTGGCTACACATATATTGACACTGAGGCACGAAAGCGTGGGGATCAAACAGG

在“>；esv2”和序列之间始终存在回车符。此外，始终存在“>；”作为标签（分类单元名称）的前缀，无任何spae。您可以简单地通过Python中的reg ex或“re”进行转换。使用perl单行程序，它将是s/^([az]+[0-9]+)/>$1/g类型的代码。我敢肯定，他们将是一个在线网站，将做到这一点

然后，只需在import命令中将“phylip”替换为“fasta”。一旦导入，您要求BioPython转换为您想要的任何格式，它应该不会有任何问题

Answer 2

首先，请阅读How to make good reproducible pandas examples的答案。今后请提供minimal reproducibl example

其次，Michael G绝对正确，phylip是一种语法非常独特的格式

下面的代码将允许您从数据框架生成一个系统发育树

首先进行一些导入，让我们重新创建数据帧

import pandas as pd
from Bio import Phylo
from Bio.Phylo.TreeConstruction import DistanceCalculator
from Bio.Phylo.TreeConstruction import DistanceTreeConstructor
from Bio import AlignIO

data = {'ESV' : ['esv1', 'esv2', 'esv3'],
        'trunc': ['TACGTAGGTGGCGAGCGTTATCCGGAATTATTGGGCGTAAAGCGCGCGCAGGCGGTTTCTTAAGTCTGATGTGAAAGCCCACGGCTCAACCGTGGAGGGTCATTGGAAACTGGGAAACTTGAGTGCAGAAGAGGAAAGCGGAATTCCACGTGTAGCGGTGAAATGCGTAGAGATGTGGAGGAACACCAGTGGCGAAGGCGGCTTTCTGGTCTGTAACTGACGCTGAGGCGCGAAAGCGTGGGGAGCAAACAGG',
                 'TACGGAGGGTGCAAGCGTTATCCGGATTCACTGGGTTTAAAGGGTGCGTAGGTGGGTTGGTAAGTCAGTGGTGAAATCTCCGAGCTTAACTTGGAAACTGCCATTGATACTATTAATCTTGAATATTGTGGAGGTTAGCGGAATATGTCATGTAGCGGTGAAATGCTTAGAGATGACATAGAACACCAATTGCGAAGGCAGCTGGCTACACATATATTGACACTGAGGCACGAAAGCGTGGGGATCAAACAGG',
                 'TACGGGGGGGGCAAGCGTTGTTCGGAATTACTGGGCGTAAAGGGCGCGTAGGCGGCCAGACCAAGTCGAGTGTGAAATTGCAGGGCTTAACTTTGCAGGGTCGCTCGATACTGGTCGGCTAGAGTGTGGAAGAGGGTACTGGAATTCCCGGTGTAGCGGTGAAATGCGTAGATATCGGGAGGAACACCAGCGGCGAAGGCGGGTACCTGGGCCAACACTGACGCTGAGGCGCGAAAGCTAGGGGAGCAAACAG']

}

ndf = pd.DataFrame.from_dict(data)
print(ndf)

输出：

    ESV                                              trunc
0  esv1  TACGTAGGTGGCGAGCGTTATCCGGAATTATTGGGCGTAAAGCGCG...
1  esv2  TACGGAGGGTGCAAGCGTTATCCGGATTCACTGGGTTTAAAGGGTG...
2  esv3  TACGGGGGGGGCAAGCGTTGTTCGGAATTACTGGGCGTAAAGGGCG...

接下来，以正确的格式编写phylip文件

with open("test.phy", 'w') as f:
    f.write("{:10} {}\n".format(ndf.shape[0], ndf.trunc.str.len()[0]))
    for row in ndf.iterrows():
        f.write("{:10} {}\n".format(*row[1].to_list()))

输出test.phy：

         3 253
esv1       TACGTAGGTGGCGAGCGTTATCCGGAATTATTGGGCGTAAAGCGCGCGCAGGCGGTTTCTTAAGTCTGATGTGAAAGCCCACGGCTCAACCGTGGAGGGTCATTGGAAACTGGGAAACTTGAGTGCAGAAGAGGAAAGCGGAATTCCACGTGTAGCGGTGAAATGCGTAGAGATGTGGAGGAACACCAGTGGCGAAGGCGGCTTTCTGGTCTGTAACTGACGCTGAGGCGCGAAAGCGTGGGGAGCAAACAGG
esv2       TACGGAGGGTGCAAGCGTTATCCGGATTCACTGGGTTTAAAGGGTGCGTAGGTGGGTTGGTAAGTCAGTGGTGAAATCTCCGAGCTTAACTTGGAAACTGCCATTGATACTATTAATCTTGAATATTGTGGAGGTTAGCGGAATATGTCATGTAGCGGTGAAATGCTTAGAGATGACATAGAACACCAATTGCGAAGGCAGCTGGCTACACATATATTGACACTGAGGCACGAAAGCGTGGGGATCAAACAGG
esv3       TACGGGGGGGGCAAGCGTTGTTCGGAATTACTGGGCGTAAAGGGCGCGTAGGCGGCCAGACCAAGTCGAGTGTGAAATTGCAGGGCTTAACTTTGCAGGGTCGCTCGATACTGGTCGGCTAGAGTGTGGAAGAGGGTACTGGAATTCCCGGTGTAGCGGTGAAATGCGTAGATATCGGGAGGAACACCAGCGGCGAAGGCGGGTACCTGGGCCAACACTGACGCTGAGGCGCGAAAGCTAGGGGAGCAAACAG

现在我们可以开始创建我们的系统发育树

# Read the sequences and align
aln = AlignIO.read('test.phy', 'phylip')
print(aln)

输出：

SingleLetterAlphabet() alignment with 3 rows and 253 columns
TACGTAGGTGGCGAGCGTTATCCGGAATTATTGGGCGTAAAGCG...AGG esv1
TACGGAGGGTGCAAGCGTTATCCGGATTCACTGGGTTTAAAGGG...AGG esv2
TACGGGGGGGGCAAGCGTTGTTCGGAATTACTGGGCGTAAAGGG...CAG esv3

计算距离矩阵：

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

输出：

esv1    0
esv2    0.3003952569169961  0
esv3    0.6086956521739131  0.6245059288537549  0

使用UPGMA算法构建系统发育树，并以ascii格式绘制树

constructor = DistanceTreeConstructor()
tree = constructor.upgma(dm)
Phylo.draw_ascii(tree)

输出：

  ________________________________________________________________________ esv3
_|
 |                                     ___________________________________ esv2
 |____________________________________|
                                      |___________________________________ esv1

或者画一幅漂亮的树图：

Phylo.draw(tree)

输出：