1

我正在尝试编写代码来为文件夹中的 100 多个单独的 fasta 对齐文件中的每一个获得一致序列。一开始我只是想获得一个序列的共识（然后我将使用一个 for 循环来处理所有序列），但是我在共识的字母表上遇到了麻烦。我的测试fasta对齐是：

>seq1
ACGTACGATCGTTACTCCTA
>seq2
ACGTACGA---TTACTCGTA

我希望达成的共识是：

ACGTACGANNNTTACTCSTA

我希望任何包含间隙的列都用“N”表示，而任何没有 100% 相同核苷酸的列都用歧义代码表示。

我的代码不起作用是：

from Bio import AlignIO
from Bio.Align import AlignInfo
from Bio.Alphabet import IUPAC, Gapped
alphabet = Gapped(IUPAC.ambiguous_dna)

alignment = AlignIO.read(open("fasta_align_for_consensus.fa"), "fasta")
summary_align = AlignInfo.SummaryInfo(alignment)
consensus = summary_align.gap_consensus(threshold = 1.0, ambiguous = 'N', consensus_alpha  \
= alphabet, require_multiple = 2)

对象“歧义”只接受一个字符串并将一个“N”放置在一致性中存在对齐多态性的任何位置，我似乎无法解决这个问题。任何有关如何纠正此问题的建议将不胜感激。谢谢！

4

1 回答 1

1

当前简单的共识方法并不能满足您的要求。听起来你在问 IUPAC 歧义代码（也许有一些阈值？）和间隙的特殊处理。您必须自己编写一些代码，也许基于现有的方法。

于 2013-05-09T17:14:00.893 回答