python - 使用 python subprocess.call 将 fasta 序列的计数写入文件

Question

我有超过 14000 个 fasta 文件，我只想保留包含 5 个序列的文件。我知道我可以使用以下 bash 命令来获取单个 fasta 文件中的序列数:

grep -c "^>" filename.fasta

所以我的方法是将每个文件中的文件名和序列计数写入一个文本文件，然后我可以使用它来仅隔离我想要的序列。要对这么多文件运行 grep 命令，我使用的是 subprocess.call:

import subprocess
import os


with open("five_seqs.txt", "w") as f:
    for file in os.listdir("/Users/vivaksoni1/Downloads/DA_CDS/fasta_files"):
        f.write(file),
        subprocess.call(["grep", "-c", "^>", file], stdout = f)

我的部分问题是 grep 命令是“^>”，但子进程要求每个参数都有自己的引号。当我本质上将作为参数输入时如何使用“^>”:“”^>”。

此外，我是否必须在 f.write(file) 之后添加 f.write("\n") ？目前我的输出只是一个文本文件，每个条目都一个接一个地出现，而 subprocess 命令只是将每个文件名打印到终端并声明没有找到这样的文件:

grep: MZ23900789.fasta: 没有那个文件或目录

Answer 1

尝试以下代码，它应该适用于您的示例。它将写入文件名加上制表符分隔符和序列数(即 > 字符)。使用 Popen 和 communicate 可以更灵活地处理输出。在 Ubuntu 上测试。

import subprocess
import os

fasta_dir = "/Users/vivaksoni1/Downloads/DA_CDS/fasta_files/"

with open("five_seqs.txt", "w") as f:
    for file in os.listdir(fasta_dir):
        f.write(file + '\t')
        grep = subprocess.Popen(["grep", "-c", "^>", fasta_dir + file], stdout = subprocess.PIPE)
        out, err = grep.communicate()
        f.write(out + '\n')