通过 Python 进行密码子对齐？答案

【问题标题】：Codon alignment via Python?通过 Python 进行密码子对齐？
【发布时间】：2013-12-30 16:52:10
【问题描述】：

我有成对的编码 DNA 序列，我希望通过 Python 执行成对的密码子比对，我已经“完成了一半”这个过程。

到目前为止..

我使用 Biopython 包从 genbank 中检索直系同源 DNA 序列对。
我将直系同源对翻译成肽序列，然后使用EMBOSS Needle 程序比对它们。

我希望..

将肽序列中的缺口转移到原始 DNA 序列中。

问题

我将不胜感激有关程序/代码（从 Python 调用）的建议，这些程序/代码可以将空位从对齐的肽序列对转移到相应核苷酸序列对的密码子上。或者可以从头开始进行成对密码子比对的程序/代码。

【问题讨论】：

好吧，我一直在做一些搜索，发现 PAL2NAL 可以满足我的需要，但需要通过网络服务器 - 但我需要通过 Python 来完成。我很想制作自己的脚本，但是我想它可能已经存在，或者可能在 Biopython 中有一些明显的功能（例如，沿着“addGapsFromPeptide()”之类的行）
有一个谷歌夏天的代码项目在 Biopython 中工作 - lists.open-bio.org/pipermail/biopython-dev/2013-July/…。也许这会有所帮助？
@Stedy 尽管有这个名字，但我不认为 CodonAlign（BioPython fork 在 GSoC2013 期间由 Zheng Ruan 开发）实际上是一个“真正的”密码子对齐器。看起来像是对密码子对齐序列进行分析的简洁工具，它还可以执行“对齐蛋白质并映射回核苷酸”技巧。有关更多信息，请参阅下面我的其他评论。

标签： python bioinformatics biopython dna-sequence sequence-alignment

【解决方案1】：

您需要做的就是将核苷酸序列分成三联体。每个氨基酸是一个三联体，每个空位是三个空位。所以在伪代码中：

for x in range(0, len(aminoacid)):
    if x != "-":
       print nucleotide[3x:3x+3]
    else:
       print "---"

【讨论】：

谢谢你，我设法用你的答案让它发挥作用:)
我也在使用这种方法，但它并不是真正的 trye 密码子对齐算法。这只是映射回密码子和it won't work 'right' in some cases 的氨基酸比对。我发现PRANK 可以工作，但是通过 Python 脚本与之交谈很尴尬。一个主题可能更适合BioStars。
抱歉上面的错字：trye => true

【解决方案2】：

您可以通过添加缺少的字符将肽映射到核苷酸：

codons = str.maketrans({'M' : 'ATG',
                        'R' : 'CGT',
                        ...,
                        '-' : '---'}) # Your missing character

peptide = 'M-R'
result = peptide.translate(codons)

然后翻译完整的序列。

【讨论】：

我认为这里的问题是氨基酸可以由多个密码子编码，即 CGT 不是唯一可以编码 R 的东西，但你给了我一些灵感来尝试我自己的脚本所以谢谢你

【解决方案3】：

我知道你三年前问过这个问题，但这篇文章是我在谷歌搜索“密码子对齐 python”中找到的第一件事。因此，我想为所有可能偶然发现仍在寻找库来执行此操作的人做出回应。

您可以为此使用 PyCogent 库。

他们在他们的网站上解释得很好：http://pycogent.org/examples/align_codons_to_protein.html

【讨论】：

【解决方案4】：

最后自己做了一个Python函数，觉得还是分享一下吧。

它需要一个对齐的肽序列和相应的未对齐的核苷酸序列，并给出一个对齐的核苷酸序列：

功能

def gapsFromPeptide( peptide_seq, nucleotide_seq ):
    """ Transfers gaps from aligned peptide seq into codon partitioned nucleotide seq (codon alignment) 
          - peptide_seq is an aligned peptide sequence with gaps that need to be transferred to nucleotide seq
          - nucleotide_seq is an un-aligned dna sequence whose codons translate to peptide seq"""
    def chunks(l, n):
        """ Yield successive n-sized chunks from l."""
        for i in xrange(0, len(l), n):
            yield l[i:i+n]
    codons = [codon for codon in chunks(nucleotide_seq,3)]  #splits nucleotides into codons (triplets) 
    gappedCodons = []
    codonCount = 0
    for aa in peptide_seq:  #adds '---' gaps to nucleotide seq corresponding to peptide
        if aa!='-':
            gappedCodons.append(codons[codonCount])
            codonCount += 1
        else:
            gappedCodons.append('---')
    return(''.join(gappedCodons))

用法

>>> unaligned_dna_seq = 'ATGATGATG'
>>> aligned_peptide_seq = 'M-MM'
>>> aligned_dna_seq = gapsFromPeptide(aligned_peptide_seq, unaligned_dna_seq)
>>> print(aligned_dna_seq)

    ATG---ATGATG

【讨论】：

@Stylize 这是我根据你的想法制作的完整脚本 :)