【问题标题】:Inability of (Bio.SeqUtils.molecular_weight function) to calculate molecular weight of unambiguous nucleotide sequence(Bio.SeqUtils.molecular_weight 函数)无法计算明确核苷酸序列的分子量
【发布时间】:2022-01-10 04:00:34
【问题描述】:

我正在尝试在 python 中创建一个函数,该函数从 fasta 文件中读取明确和不明确的核苷酸序列并返回序列 ID 和分子量。

我尝试了以下代码:

import Bio
from Bio import SeqUtils, SeqIO

def function(filename):
    nucleotides ={'A','T', 'C', 'G'}
    with open(filename) as file:
        for record in SeqIO.parse(file, "fasta"):
            for nucl in record:
                if nucl in nucleotides:
                    continue
                else:
                    print(str(record.id)+": is ambiguous")
                    break
            else:
                mol_weight= Bio.SeqUtils.molecular_weight(record)
                print(str(record.id)+": is unambiguous & molecular weight = "+str(mol_weight))

function("random.fasta")

如果我在模棱两可的序列上使用这段代码,绝对没有问题,我会得到我想要的结果。但是,如果我包含明确的序列,我会得到 "ValueError: 'I' is not a valid unambiguous letter for DNA" 由 'Bio.SeqUtils.molecular_weight(record)' 函数引起,在我看来'没有意义,因为如果我正确理解,字母 'I' 通常会导致第一个 else 语句中的中断。

我还使用了一个自定义函数来手动计算分子量(基于固定值),在这种情况下没有错误,并且该函数对于模棱两可和明确的序列都很好。然而有人指出,我的手动计算不如内置函数准确。

我对 Python 有点陌生,但我希望有人知道为什么 ValueError 仍然会出现以及如何解决它。

【问题讨论】:

  • 我修复了函数体的缩进 - 在表单中粘贴代码时很容易出错。我还包括了使代码运行所需的导入 - 所需要的只是一个有效的 fasta 文件。

标签: python sequence biopython fasta


【解决方案1】:

我认为错误发生在表达式mol_weight= Bio.SeqUtils.molecular_weight(record)

记录被隐式转换为它的字符串表示ID: SEQUENCE_1 ...。所以 ValueError 中的'I' 根本不对应一个核苷酸,而是一个相当随机的字符串值。

试试这一行:mol_weight = Bio.SeqUtils.molecular_weight(record.seq)

对我来说,它返回了想要的结果:

SEQUENCE_1: is unambiguous & molecular weight = 331.2218

对于一个简单的序列(只是字母 A),它会在引发您描述的错误之前。

>SEQUENCE_1
A

更新:

完整的解决方案如下所示:

from Bio import SeqIO, SeqUtils


def function(filename):
    with open(filename) as file:
        for record in SeqIO.parse(file, "fasta"):
            try:
                mol_weight = SeqUtils.molecular_weight(record.seq)
                print(f"{record.id}: is unambiguous & molecular weight = {mol_weight}")
            except ValueError as exception:
                print(f"{record.id}: {exception}")


if __name__ == '__main__':
    function("random.fasta")

我从@pippo1980 的回答中加入了执行处理,因为它确实更简单、更快,并且可以说更pythonic。

【讨论】:

  • 非常感谢!毕竟,这对我来说是一个非常愚蠢的错误。
  • 并非如此。由于鸭子类型系统,使用 Python 很容易犯这样的错误。我必须启动调试器才能发现这个。如果它解决了您的问题,也请将答案标记为“已接受的答案”。
【解决方案2】:

更短:

from Bio import SeqIO, SeqUtils



def function(filename):
    with open(filename) as file:
        for record in SeqIO.parse(file, "fasta"):
            try:
                mol_weight= SeqUtils.molecular_weight(record.seq, 'DNA')
                print(str(record.id)+": is unambiguous & molecular weight = "+str(mol_weight))
            
            except Exception as exception:
                print(record.id+':',str(type(exception).__name__), str(exception))
            

function("random.fasta")
            

function("random.fasta")

输入random.fasta文件:

>first
AAAAAAAAAAAAATTTTTTTTTTTTTTGGGGGGGGGGGGGGCCCCCCCCCCCCCCGTA
>second
AAAAAAAAAAAAAAAAAAAAATTTTTTTTTTTTTTTTICCCCCCCCCCCCCCCCCCCCCCCCCCC
>third
AAAAAAAAAQQQQQQQQQQQQQQQQQQTTTTTTTTTT
>fourth
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAATTTTTTTTTTTTTTTTTGGGGGGGGGGGGGGGGGCGCGCGCGCGC

输出:

first: is unambiguous & molecular weight = 17952.438
second: ValueError 'I' is not a valid unambiguous letter for DNA
third: ValueError 'Q' is not a valid unambiguous letter for DNA
fourth: is unambiguous & molecular weight = 25755.56080000001

【讨论】: