Python/Unix高效删除超大FASTA文件指定行及大文件处理技巧

问题描述 #

我是生物信息学从业者，当前需处理大小为64.2GB的FASTA格式人类参考基因组文件GRCh38.fa，经wc统计该文件参数为：1052454251 1052456168 64199706147 GRCh38.fa。核心需求是删除文件中以>chrM开头的线粒体序列标题行，及其下属所有对应的DNA序列行，保留其余染色体的序列内容。

此前尝试两类方案均存在问题：

• 使用readlines()一次性读取全文件，因内存占用过高导致系统挂起
• 采用逐行读取、过滤写入新文件的方案，内存占用可控但处理64GB文件耗时过长

目前检索到的Python删除文件指定行方案几乎都需要全量重写文件，因此想确认是否存在可直接删除文件指定行的专用Python/Unix命令。此外日常工作多围绕超大型数据集分析展开，希望获取Python编写内存高效、时间高效代码处理大文件的相关技巧，解决代码因RAM不足被系统杀死的问题。

原有逐行过滤代码如下：

from memory_profiler import profile  

@profile                       #This is just to check how much memory is used.
def memcheck():
    g=open("chrMremoved.fa",'w')
    to_write=1
    with open("GRCh38.fa") as f:
        for i in f:
            if(i[0]==">"):
                if(i[0:5]==">chrM"):
                    to_write=0
                else:
                    to_write=1
            if(to_write==1):
                g.write(i)
    g.close()

if __name__ == "__main__":
    memcheck()

解决方案 #

注意：不存在不需要重写文件就能直接删除中间指定内容的通用方法。磁盘上的文件是以连续块形式存储的，删除中间段内容后，后续所有保留的内容块都需要前移覆盖，所有宣称"原地修改"的工具本质都是写临时文件后替换原文件，不需要在这个方向上浪费时间找方案。

最快过滤方案（Unix流工具，优先选） #

纯文本过滤类任务优先用系统自带的C实现流处理工具，性能是手写Python的2-5倍，内存占用稳定在KB级，完全不会出现内存溢出问题：

1. mawk方案（速度最快）：mawk是做了极致性能优化的awk实现，处理64GB文件在SSD上仅需2-3分钟，机械盘约8-10分钟，命令如下：

   mawk 'BEGIN{write=1} /^>/{write=!/^>chrM/} write' GRCh38.fa > chrM_removed.fa
   

2. sed方案（系统自带无需额外安装）：大多数Linux发行版默认预装sed，不需要额外装包，命令如下：

   sed '/^>chrM/,/^>/{/^>/!d; /^>chrM/d}' GRCh38.fa > chrM_removed.fa
   

两个命令的逻辑完全匹配需求：识别>chrM开头的线粒体记录起点，直到下一条序列的标题行出现为止，跳过区间内所有线粒体序列内容，其余内容全部写入新文件。

Python代码优化方案 #

原有Python代码速度慢主要是三个原因：一是加载了memory_profiler装饰器，该工具会给每行代码插入内存监控钩子，带来30%以上的性能损耗；二是未手动设置读写缓冲区，默认小缓冲区导致IO次数过多；三是条件判断冗余。优化后代码性能可达到mawk的60%-70%，比原代码快2-3倍：

def filter_chrm(input_path: str, output_path: str) -> None:
    # 设置64MB读写缓冲区，大幅减少磁盘IO切换开销
    buf_size = 64 * 1024 * 1024
    to_write = True
    with open(input_path, 'r', buffering=buf_size) as fin, \
         open(output_path, 'w', buffering=buf_size) as fout:
        for line in fin:
            # 遇到新序列标题行时，判断是否为chrM，更新写入标记
            if line.startswith('>'):
                to_write = not line.startswith('>chrM')
            if to_write:
                fout.write(line)

if __name__ == "__main__":
    filter_chrm("GRCh38.fa", "chrMremoved.fa")

超大型文件处理通用高效技巧 #

• 永远不要用read()、readlines()一次性加载全文件，Python的文件对象本身是迭代器，逐行/逐块遍历只会把当前处理的块加载到内存，内存占用完全可控
• 大文件读写手动设置16MB-128MB的缓冲区，默认的几KB缓冲区会导致频繁的内核态/用户态切换，IO效率极低
• 处理大文件的生产任务时，移除所有调试类装饰器、冗余打印语句，这类非业务逻辑的性能损耗往往远超预期
• 简单的文本过滤、替换任务优先用sed、mawk、grep等系统流工具，这类工具经过几十年的优化，纯文本处理性能远高于通用脚本语言
• 必须用Python实现复杂处理逻辑时，优先用内置字符串方法（比如startswith()、endswith()）做匹配，不要用正则处理简单的前缀/后缀判断，内置方法是C实现，速度比正则快一个数量级
• 处理无换行结构的二进制文件时，不要逐行读取，每次固定读取64MB左右的块，在块内完成逻辑处理，性能提升明显
• 循环内不要放不变的逻辑、重复的变量赋值，所有固定计算全部提到循环外部，减少循环内的运算量
• 单线程性能遇到瓶颈时，可以按记录边界切分文件（比如FASTA按>开头的序列切分），用多进程并行处理后合并结果，不要直接按固定字节数切分避免破坏格式结构

内容的提问来源于stack exchange，提问作者Joshua Jebaraj