TechOverflow Logo
TechOverflow

一个简单的 FASTA 统计工具

English Deutsch 简体中文

问题

使用现有软件计算 FASTA 文件的简单统计出奇地困难。我最近需要计算 FASTA 格式单个序列的核苷酸计数和相对 GC 频率，但除非你安装依赖繁重的原生软件如 FASTX或你自己使用 BioPython 或类似工具开发，否则似乎没有针对这组简单问题的简单、无依赖解决方案。

解决方案：fasta-stats.py

我开发了一个简单的 Python 脚本，仅使用标准 Python 库，能够为一个或多个 FASTA 文件（明文或 Gzipped）生成这些统计，每个文件有一个或多个序列。计算字符（字符通常是氨基酸代码或核苷酸）的绝对和相对频率，并为输入文件的每个序列打印。如果文件包含核苷酸序列（即只包含 ATGC 字符），也会自动打印 A+T 和 G+C 计数。

默认情况下，字符以不区分大小写的方式计数

用法示例：

fasta_stats_usage.sh
#多个文件的简单统计，自动检测 .gz 输入
fasta-stats.py sequence1.fa sequence.2.fa.gz
#区分大小写计数字符
fasta-stats.py --case-sensitive sequence1.fa sequence.2.fa.gz
#仅计数 ATGC 字符，忽略所有其他字符
fasta-stats.py --only ATGC sequence1.fa sequence.2.fa.gz

脚本源码：

fasta_stats.py
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
fasta-stats.py：计算 FASTA 文件中核苷酸/氨基酸数量的实用脚本。
更新日志：
    1.1：Python3 就绪
"""
from __future__ import with_statement
import sys
import argparse
import gzip
#Counter 用于每字符统计
from collections import Counter
__author__    = "Uli Koehler & Anton Smirnov"
__copyright__ = "Copyright 2013 Uli Koehler & Anton Smirnov"
__license__   = "Apache v2.0"
__version__   = "1.1"

def printSequenceStats(fileName, sequenceName, charOccurrenceMap, totalCharCount):
    """
    将序列详情打印到 stdout。
    从 parseFile() 内部调用。

    关键字参数：
        sequenceName：要打印的序列名称
        charOccurrenceMap：包含 char --> 出现次数映射的类字典对象
        totalCharCount：序列的总体字符计数
    """
    print ("Sequence '{}' from FASTA file '{}' contains {} sequence characters:".format(
        sequenceName, fileName, totalCharCount))
    for char in sorted(charOccurrenceMap.keys()):
        charCount = charOccurrenceMap[char]
        relativeFrequency = charCount * 100.0 / totalCharCount
        print ("\t{} : {} = {}%".format(char, charCount, relativeFrequency))
    #对于核苷酸序列（仅 ATGC），还打印 A+T vs G+C 计数
    if sorted(charOccurrenceMap.keys()) == ["A","C","G","T"]:
        #打印 A+T 计数
        atCount = charOccurrenceMap["A"] + charOccurrenceMap["T"]
        atRelFrequency = atCount * 100.0 / totalCharCount
        print ("\tA+T : {} = {}%".format(atCount, atRelFrequency))
        #打印 G+C 计数
        gcCount = charOccurrenceMap["G"] + charOccurrenceMap["C"]
        atRelFrequency = gcCount * 100.0 / totalCharCount
        print ("\tG+C : {} = {}%%".format(gcCount, atRelFrequency))

def parseFile(filename, caseSensitive=False, charWhitelist=None):
    """
    解析 FASTA 文件并调用 printRe
    """
    #设置为标题行，从开头移除 ">"
    sequenceName = None
    #键：字符，值：出现次数
    charOccurrenceMap = Counter()
    #当前序列中的字符数，不包括 \n
    charCount = 0
    #跟踪连续注释，因为它们被追加
    previousLineWasComment = False
    #打开并迭代文件，自动检测 gzip
    openFunc = gzip.open if filename.endswith(".gz") else open
    with openFunc(filename, "r") as infile:
        for line in infile:
            line = line.strip()
            #与原始规范超级兼容
            if line.startswith(">") or line.startswith(";"):
                #处理前一个序列（如果有）
                if sequenceName is not None:
                    printSequenceStats(filename, sequenceName, charOccurrenceMap, charCount)
                    charOccurrenceMap = Counter()
                    charCount = 0
                #将整个注释行作为（新）序列 ID（去除 ">"）
                #连接连续序列行
                if previousLineWasComment: #追加 -- 在中间添加一个空格以规范化空格计数
                    sequenceName += " " + line[1:].strip()
                else:
                    sequenceName = line[1:].strip()
                previousLineWasComment = True
            else: #行属于序列
                previousLineWasComment = False
                #行之前已被剥离，所以我们可以直接计数
                #递增每字符统计（字符出现）
                for char in line:
                    #跳过白名单中不存在的任何字符，如果启用了白名单（--only）
                    if charWhitelist is not None and not char in charWhitelist:
                        continue
                    #我们只能在白名单过滤后计数
                    charCount += 1
                    #在不区分大小写模式（默认）下仅计数大写字符
                    char = char if caseSensitive else char.upper()
                    charOccurrenceMap[char] += 1
        #最后一行已读取，打印最后一条记录（如果有）
        if sequenceName is not None:
            printSequenceStats(filename, sequenceName, charOccurrenceMap, charCount)

if __name__ == "__main__":
    #允许指定单个或多个文件
    parser = argparse.ArgumentParser(description='Compute simple statistics for FASTA files.')
    parser.add_argument('infiles', nargs='+', help='要生成统计的 FASTA 文件（.fa, .fa.gz）')
    parser.add_argument('--case-sensitive', action='store_true', help='以区分大小写的方式计数字符。默认禁用。')
    parser.add_argument('-o','--only', help='如果提供此选项（例如设置为 \'ATGC\'），不在集合中的字符将在所有统计中被忽略')
    args = parser.parse_args()
    #处理所有 FASTA 文件
    for infile in args.infiles:
        parseFile(infile, caseSensitive=args.case_sensitive, charWhitelist=args.only)
Check out similar posts by category: Bioinformatics, Python