使用 biopython 将 dna 对齐转换为 numpy 数组答案

【问题标题】：Transform dna alignment into numpy array using biopython使用 biopython 将 dna 对齐转换为 numpy 数组
【发布时间】：2016-09-29 20:32:35
【问题描述】：

我有几个已对齐的 DNA 序列，我想只保留在特定位置可变的碱基。

如果我们首先将对齐方式转换为数组，这也许可以做到。我尝试使用 Biopython 教程中的代码，但它给出了错误。

import numpy as np
from Bio import AlignIO
alignment = AlignIO.parse("ma-all-mito.fa", "fasta")
align_array = np.array([list(rec) for rec in alignment], np.character)
print("Array shape %i by %i" % align_array.shape)

我得到的错误：

Traceback (most recent call last):

File "C:/select-snps.py", line 8, in <module>
    print("Array shape %i by %i" % align_array.shape)
TypeError: not all arguments converted during string formatting

【问题讨论】：

你试过简单的print(align_array.shape)吗？
输出 (1, 99, 16926)
所以它可能奏效了

标签： python numpy biopython

【解决方案1】：

AlignIO 似乎不是您想要完成这项工作的工具。您的文件可能包含许多序列，而不是许多多序列比对，因此您可能希望使用SeqIO，而不是AlignIO (source)。这就是数组的形状为 (1, 99, 16926) 的原因，因为您有 99 个长度为 16926 的序列的 1 个比对。

如果您只想要一个序列数组（看起来您是从提供给np.array 的np.character dtype 执行的），请执行以下操作：

import numpy as np
from Bio import SeqIO
records = SeqIO.parse("ma-all-mito.fa", "fasta")
align_array = np.array([record.seq for record in records], np.character)
print("Array shape %i by %i" % align_array.shape)
# expect to be (99, 16926)

请注意，从技术上讲，records 的每个元素也是一个 BioPython SeqRecord，除了元数据之外还包括序列。 list(record) 是获取序列的快捷方式，另一种方式是record.seq。两者都可以，但我选择使用属性方式，因为它更明确。

【讨论】：

【解决方案2】：

我是在回答您的问题，而不是修复您的代码。如果你只想保留某些位置，你想使用AlignIO：

FASTA 样本al.fas:

>seq1
CATCGATCAGCATCGACATGCGGCA-ACG
>seq2
CATCGATCAG---CGACATGCGGCATACG
>seq3
CATC-ATCAGCATCGACATGCGGCATACG
>seq4
CATCGATCAGCATCGACAAACGGCATACG

现在假设您只想保留某些位置。 MultipleSeqAlignment 允许您像 numpy 数组一样查询对齐：

from Bio import AlignIO


al = AlignIO.read("al.fas", "fasta")

# Print the 11th column
print(al[:, 10])

# Print the 12-15 columns
print(al[:, 11:14])

如果您想知道对齐的形状，请使用len 和get_alignment_length：

>>> print(len(al), al.get_alignment_length())
4 29

当您使用AlignIO.parse() 加载对齐时，它假定要解析的文件可能包含多个对齐（PHYLIP 会这样做）。因此，解析器返回每个对齐的迭代器，而不是代码所暗示的记录。但是您的 FASTA 文件每个文件只包含一个对齐，parse() 只产生一个 MultipleSeqAlignment。所以对你的代码的修复是：

alignment = AlignIO.read("ma-all-mito.fa", "fasta")
align_array = np.array(alignment, np.character)
print("Array shape %i by %i" % align_array.shape)

【讨论】：