将一个多FASTA文件拆分成具有相同数量的访问号的文件。

Question

我有一个包含成千上万个存取号码的文件：

长这样...

>NC_033829.1 Kallithea virus isolate DrosEU46_Kharkiv_2014, complete genome
AGTCAGCAACGTCGATGTGGCGTACAATTTCTTGATTACATTTTTGTTCCTAACAAAATGTTGATATACT

>NC_020414.2 Escherichia phage UAB_Phi78, complete genome
TAGGCGTGTGTCAGGTCTCTCGGCCTCGGCCTCGCCGGGATGTCCCCATAGGGTGCCTGTGGGCGCTAGG

如果想将此文件拆分为每个接入号一个文件，那么我可以使用以下代码。

awk -F '|' '/^>/ {F=sprintf("%s.fasta",$2); print > F;next;} {print >> F;}' < yourfile.fa

我有一个包含成千上万个接入号（也称为>NC_*）的文件，想要将其分割，使每个文件包含大约5000个接入号。由于我是awk/bash/python的新手，因此我很难找到一个简洁的解决方案。

欢迎提出任何想法或评论。

- LDT

Answer 1

从您的问题中并不清楚"accession number"是每个输入块唯一的（不要假设阅读您问题的人对您的领域有任何了解-对我们来说，它只是文本行）。如果您的问题仅表达您想要每个输出文件5000个新行分隔块，而不是5000个接入号码，那么问题将更清晰明了。

看到您发布的答案后，现在清楚这是您应该使用的内容。

awk -v RS= -v ORS='\n\n' '
    (NR%5000) == 1 { close(out); out="myseq"(++n_seq)".fa" }
    { print > out }
' my_sequences.fa

- Ed Morton - SO stop bullying

---

谢谢@Ed！这非常有帮助，而且讲解得很清楚。我感谢你的时间。 - LDT

Answer 2

最好使用Biopython的Bio.SeqIO来处理读取和写入FASTA文件。然后，您只需要以所需的方式对记录（SeqRecord对象）进行分组即可。我更喜欢将分组函数产生迭代器：

from itertools import chain, islice

from Bio import SeqIO


def grouper(n, iterable):
    it = iter(iterable)
    while True:
        chunk_it = islice(it, n)
        try:
            first = next(chunk_it)
        except StopIteration:
            return
        yield chain((first,), chunk_it)


for idx, group in enumerate(grouper(5000, SeqIO.parse('input.fa', 'fasta')), 1):
    SeqIO.write(group, f'out-{idx}.fa', 'fasta')

- Steve

---

你为什么需要在这里使用 grouper 函数？你能直接使用 itertools.islice() 吗？ - Chris_Rands

@Chris_Rands 好吧，itertools.islice()只是返回所选元素的迭代器。需要的是一种遍历所有元素的方式，因此需要这个函数。 - Steve

Answer 3

假设：章节之间由空行分隔。

算法：

• 将文件分成章节
• 从章节中提取访问号
• 将章节输出到以访问号命名的文件中。

Awk术语：一个“记录”将是我们的章节——由两个换行符分隔的文件部分。一个“字段”通常由空格分隔——通过用空格或>字符分隔，第二个字段将是访问号。

只需将记录分隔符设置为两个换行符，字段分隔符设置为>或空格，然后将该行输出到以第二个字段命名的文件中：

awk -v RS='' -v FS='[> ]' '{f=($2 ".txt"); print >> f; close(f)}'

@edit将>更改为>>，将RS='\n\n'更改为RS=''

@edit还添加了关闭

- KamilCuk

---

我认为我还没有理解如何在同一文件中维护多个接入号。比方说，如果我有四个NC，我该如何分割我的文件，以便将两个NC放到一个文件中，另外两个NC放到另一个文件中？我需要改变哪个参数？对于我的愚蠢问题，我感到很抱歉。 - LDT

@LDT，你是在说你接受的答案实际上并不能满足你的需求吗？它会按照你在问题中所描述的去执行，所以如果你想要的不是你所问的，请提出一个新的问题。 - Ed Morton

亲爱的 @EdMorton，没错。我努力添加了解决方案。非常感谢你的时间，我学到了很多，并为你的答案点赞。 - LDT

@LDT 不用谢，但我之前还没有发布答案。现在我已经发布了。 - Ed Morton

Answer 4

非常感谢您的回答，我学到了很多。我真正想做的是将多重fasta文件拆分成具有相同访问号数量的文件。在经过长时间的奋斗和同事的帮助后，这是我的答案。

awk 'BEGIN {n_seq=0;} /^>/ {if(n_seq%5000==0){file=sprintf("myseq%d.fa",n_seq);} print >> file; n_seq++; next;} { print >> file; }' < my_sequences.fa

在这里，您可以创建新的fasta文件，每个文件都有5000个访问号，也称为头文件。

谢谢大家。

- LDT

---

您不需要BEGIN部分，也不应该在它们之间有一个next的相同打印语句，因为那是冗余代码。如果您的输入文件很大，在某些awk中会出现“太多打开的文件”的错误，而在gawk中，由于您没有在进行操作时关闭输出文件，因此速度会变慢。同时，从shell重定向输入并不是必要的，awk完全可以打开一个输入文件，然后使用FILENAME。 - Ed Morton