生物塞隆排列的无上限指数

Question

我第一次用生物圈。如果这是个基本问题，请原谅我。在

我想输入序列，然后对齐它们，并且能够引用原始序列（uncapped）和对齐序列（gaped）的索引位置。在

我的真实例子是烯醇化酶（UniprotP37869和P0A6P9）。底物结合赖氨酸在大肠杆菌中为392，在枯草芽孢杆菌中为389。如果一个人做了两个肌肉排列，这个赖氨酸在排列中的指数是394。我想能够很容易地在gappend和ungapped之间转换。在

例1:大肠杆菌残基的对齐指数是多少？（按排列顺序回答394）。在

例2我在394处发现一个保守残基。在原始（未封顶）序列中的哪一个？（在大肠杆菌中回答392，在枯草芽孢杆菌中回答389）。在

>>>cline = MuscleCommandline(input="enolase.txt", out="enolase.aln")
>>>cline()

>>> alignment = AlignIO.read("enolase.aln","fasta")
>>> print(alignment[:,390:])
SingleLetterAlphabet() alignment with 2 rows and 45 columns
AGQIKTGAPSRTDRVAKYNQLLRIEDQLAETAQYHGINSFYNLNK sp|P37869|ENO_BACSU
AGQIKTGSMSRSDRVAKYNQLIRIEEALGEKAPYNGRKEIKGQA- sp|P0A6P9|ENO_ECOLI
>>> print(alignment[:,394])
KK

Answer 1

有趣的问题！据我所知，没有什么东西内置在BioPython中。下面是我将如何解决它。在

让我们从示例2开始。如果您将您的两个文件enolase.txt和{}分别与原始的未封顶序列和对齐序列一起以FASTA格式保存，那么我们可以循环压缩的记录，计算对齐序列中的间隙数，并计算未封顶序列中剩余部分的索引：

from Bio import SeqIO, AlignIO

def find_in_original(index, original_path, alignment_path):
    alignment = AlignIO.read(alignment_path, 'fasta')
    original = SeqIO.parse(original_path, 'fasta')

    for original_record, alignment_record in zip(original, alignment):
        alignment_seq = str(alignment_record.seq)
        original_seq = str(original_record.seq)

        gaps = alignment_seq[:index].count('-')
        original_index = len(alignment_seq[:index]) - gaps
        assert alignment_seq[index] ==  original_seq[original_index]
        yield ("The conserved residue {} at location {} in the alignment can be"
               " found at location {} in {}.".format(alignment_seq[index],
               index, original_index, original_record.id.split('|')[-1]))

结果是：

^{pr2}$

对于反向操作，我们查看对齐中的所有可能的索引，如果减去间隙，看看哪个索引等于未封顶的序列：

def find_in_alignment(index, organism, original_path, alignment_path):
    alignment = AlignIO.read(alignment_path, 'fasta')
    original = SeqIO.parse(original_path, 'fasta')

    for original_record, alignment_record in zip(original, alignment):
        if organism in original_record.id:
            alignment_seq = str(alignment_record.seq)
            original_seq = str(original_record.seq)

            residue = original_seq[index]
            for i in range(index, len(alignment_seq)):
                if alignment_seq[i] == residue and \
                   alignment_seq[:i].replace('-', '') == original_seq[:index]:
                    return ("The residue {} at location {} in {} is at location"
                            " {} in the alignment.".format(residue, index,
                            organism, i))

结果是：

^{4}$