使用 Biopython 基于 IDS 过滤 FASTA 文件答案

【问题标题】：Filtering FASTA file based on IDS with Biopython使用 Biopython 基于 IDS 过滤 FASTA 文件
【发布时间】：2016-11-22 16:40:29
【问题描述】：

我对 Python 编程非常陌生。我有包含某些植物物种蛋白质序列的 fasta 文件。

我想根据每个序列包含的氨基酸数量过滤它们。标准是那些>20个氨基酸的序列。

我可以通过biopython cookbook 上的资源获得超过 20 个氨基酸序列。但是，当我尝试将它们写入文件时，它给了我这个Error。我无法解决此错误。此外，我还想在输出文件中包含每个序列的 ID。请帮我！

代码：

import Bio
from Bio import SeqIO
for s_record in SeqIO.parse('arabidopsis_thaliana_proteome.ath.tfa','fasta'):
    name = s_record.id
    seq = s_record.seq
    seqLen = len(s_record)
    if seqLen >20:
        desired_proteins=seq
        output_file=SeqIO.write(desired_proteins, "filtered.fasta","fasta")
output_file

输入文件：Arabidopsis Thaliana

>AT5G16970
MTATNKQVILKDYVSGFPTESDFDFTTTTVELRVPEGTNSVLVKNLYLSCDPYMRIRMGKPDPSTAALAQAYTPGQPIQGYGVSRIIESGHPDYKKGDLLWGIVAWEEYSVITPMTHAHFKIQHTDVPLSYYTGLLGMPGMTAYAGFYEVCSPKEGETVYVSAASGAVGQLVGQLAKMMGCYVVGSAGSKEKVDLLKTKFGFDDAFNYKEESDLTAALKRCFPNGIDIYFENVGGKMLDAVLVNMNMHGRIAVCGMISQYNLENQEGVHNLSNIIYKRIRIQGFVVSDFYDKYSKFLEFVLPHIREGKITYVEDVADGLEKAPEALVGLFHGKNVGKQVVVVARE*

>AT4G32100
MATNACKFLCLVLLFAFVTQGYGDDSYSLESLSVIQSKTGNMVENKPEWEVKVLNSSPCYFTHTTLSCVRFKSVTPIDSKVLSKSGDTCLLGNGDSIHDISFKYVWDTSFDLKVVDGYIACS*

提前谢谢你:)

【问题讨论】：

你得到的错误信息是什么？
AttributeError: 'str' 对象没有属性 'id'
请输出准确的错误信息，以及复制错误的示例输入。
嗨@Vince，我已经编辑了问题，带有错误快照和示例输入。
你已经设置了desired_proteins = seq = s_record.seq，这是一个类似字符串的Seq对象。 write 函数需要一个 SeqRecord 对象，即您的示例中的 s_record 应该可以工作。

标签： python-2.7 bioinformatics biopython

【解决方案1】：

根据此处的 BioPython 教程：

http://biopython.org/wiki/SeqIO

SeqIO.parse() 的第一个参数应该是文件句柄，而不是文件名：

from Bio import SeqIO
with open("example.fasta", "rU") as handle:
    for record in SeqIO.parse(handle, "fasta"):
        print(record.id)

这应该可行：

import Bio
from Bio import SeqIO
fh=open('arabidopsis_thaliana_proteome.ath.tfa')
for s_record in SeqIO.parse(fh,'fasta'):
    name = s_record.id
    seq = s_record.seq
    seqLen = len(s_record)
    if seqLen >20:
        desired_proteins=seq
        output_file=SeqIO.write(desired_proteins, "filtered.fasta","fasta")
output_file
fh.close()

【讨论】：

SeqIO.parse 已经接受文件名或句柄很长时间了，但很好的一点 - 网页可以更明确地说明这一点。
这个解决方案和原始解决方案都会在每次调用 write 时保持覆盖过滤。fasta，所以最终结果是它只包含最后一个长蛋白质。在这里，最简单的方法是在 for 循环之前打开一个输出句柄，然后使用它。
更新：SeqIO 页面已更新以显示带有文件名的示例

【解决方案2】：

使用bioawk（对生物信息学有用的 awk 的修改版本）而不是 biopython 的解决方案：

bioawk -c fastx 'length($seq) > 20 {print ">"$name"\n"$seq}' arabidopsis_thaliana_proteome.ath.tfa > filtered.fasta

-c fastx 告诉 bioawk 将文件解析为 fastx/fastq 格式。这定义了一个 name 和一个 seq 变量，可以使用普通的 'condition {action}' awk 语法来使用。

我不确定 bioawk 是否为大多数 linux 发行版打包，但从源代码安装并不太复杂： https://silico-sciences.com/2015/12/13/install-bioawk-on-ubuntu/

【讨论】：

【解决方案3】：

在http://biopython.org/wiki/SeqIO http://biopython.org/wiki/SeqIO 上有一个示例（过滤长度小于 300 的序列），“输入/输出示例 - 按序列长度过滤”

这涉及在内存中建立记录列表的选项，而不是内存高效的生成器表达式。

第三种选择是使用输出句柄多次调用 SeqIO.write：

from Bio import SeqIO

with open("short_seqs.fasta", "w") as out_handle:
    for record in SeqIO.parse("cor6_6.gb", "genbank"):
        if len(record.seq) < 300:
            SeqIO.write(record, out_handle, "fasta")

但是，这种方法仅适用于简单的文件格式，如 FASTA、GenBank、EMBL，您可以在这些格式中继续将记录附加到现有文件中。这不适用于 XML 输出格式之类的东西 - 因此链接页面上显示的方法更好（使用记录的生成器/迭代器对 SeqIO.write 进行一次调用）。

【讨论】：