如何统计FASTA文件中各23bp序列的非重叠出现次数并添加至TSV(BED)文件新列

针对你处理大规模BED和FASTA文件的需求，我来梳理一套高效的解决方案：

你的核心需求分为两点：

• 统计BED文件第7列的每个23bp序列在hg19全基因组FASTA中的非重叠出现次数（需支持跨FASTA行的序列匹配）
• 将统计结果追加到BED文件对应行的新列，输出指定的TSV格式

你已经能用grep -o <seq> <fasta> | grep -c ""统计单个序列的出现次数，但这套方案在批量处理和大文件场景下存在明显瓶颈：

• grep/awk原生无法高效处理跨多行的序列匹配，且如果遍历50万条序列逐一扫描6000万行的FASTA，重复IO会导致效率极低
• Perl方案无法适配超大规模文件，内存占用和处理速度都会出现瓶颈

结合你提到的Ted Lyngmo的协作方向，这里整理核心实现要点，保证处理效率适配你的文件规模：

1. FASTA预处理：
   • 读取FASTA时，跳过所有以>开头的头部行，将每条染色体的序列拼接成连续的字符串（消除换行符），彻底解决跨分行匹配的问题
   • 现代服务器内存足够直接加载hg19全基因组（约3GB），如果内存有限也可以按染色体分块存储，按需加载
2. 序列匹配与计数：
   • 对每个23bp目标序列，在每条染色体的连续序列中执行非重叠滑动匹配：匹配到目标序列后，直接跳过23bp继续查找，直到序列末尾
   • 利用C++标准库的string::find循环即可实现高效匹配，也可以用KMP算法进一步优化（对短序列来说，find的效率已经足够）
3. 批量处理BED文件：
   • 逐行读取BED文件，提取第7列的序列，先收集所有序列并去重，避免重复统计相同序列（如果BED中有大量重复序列，这一步能大幅节省时间）
   • 用哈希表（std::unordered_map）存储每个序列的总计数，最后遍历原始BED行，将对应计数追加到行尾输出

#include <iostream>
#include <fstream>
#include <string>
#include <unordered_map>
#include <vector>
#include <sstream>

// 读取FASTA文件，返回每条染色体的连续序列
std::vector<std::string> read_fasta(const std::string& fasta_path) {
    std::vector<std::string> chroms;
    std::ifstream fasta(fasta_path);
    std::string line, current_chrom;

    while (std::getline(fasta, line)) {
        if (line.empty()) continue;
        // 跳过染色体头部行
        if (line.front() == '>') {
            if (!current_chrom.empty()) {
                chroms.push_back(current_chrom);
                current_chrom.clear();
            }
            continue;
        }
        // 拼接连续序列
        current_chrom += line;
    }
    // 加入最后一条染色体的序列
    if (!current_chrom.empty()) chroms.push_back(current_chrom);
    return chroms;
}

// 统计单个序列在所有染色体中的非重叠出现次数
int count_non_overlapping(const std::string& target_seq, const std::vector<std::string>& chroms) {
    int total_count = 0;
    const size_t seq_len = target_seq.size();

    for (const auto& chrom_seq : chroms) {
        size_t current_pos = 0;
        // 循环查找非重叠匹配
        while ((current_pos = chrom_seq.find(target_seq, current_pos)) != std::string::npos) {
            total_count++;
            current_pos += seq_len;
        }
    }
    return total_count;
}

int main(int argc, char* argv[]) {
    if (argc != 3) {
        std::cerr << "Usage: " << argv[0] << " <input_bed> <hg19_fasta>\n";
        return 1;
    }

    // 1. 读取并预处理FASTA文件
    auto chrom_sequences = read_fasta(argv[2]);

    // 2. 读取BED文件，收集所有序列并去重统计
    std::ifstream bed_file(argv[1]);
    std::unordered_map<std::string, int> seq_count_map;
    std::vector<std::vector<std::string>> bed_records; // 保存原始BED的每行字段

    std::string line;
    while (std::getline(bed_file, line)) {
        std::vector<std::string> fields;
        std::stringstream ss(line);
        std::string field;
        // 分割TSV字段
        while (std::getline(ss, field, '\t')) {
            fields.push_back(field);
        }
        bed_records.push_back(fields);
        const std::string& target_seq = fields[6];
        // 仅统计未处理过的序列
        if (seq_count_map.find(target_seq) == seq_count_map.end()) {
            seq_count_map[target_seq] = count_non_overlapping(target_seq, chrom_sequences);
        }
    }

    // 3. 输出带计数的BED文件
    for (const auto& record : bed_records) {
        for (size_t i = 0; i < record.size(); ++i) {
            if (i > 0) std::cout << '\t';
            std::cout << record[i];
        }
        std::cout << '\t' << seq_count_map[record[6]] << '\n';
    }

    return 0;
}

• 编译：g++ -O3 -o seq_counter seq_counter.cpp（-O3开启最高优化，大幅提升处理速度）
• 运行：./seq_counter your_input.bed hg19.fasta > output_with_counts.bed

这套方案的核心优势：

• 仅扫描FASTA文件一次，后续所有统计都基于内存中的连续序列，彻底避免重复IO开销
• 非重叠匹配逻辑准确，完美支持跨FASTA行的序列匹配
• 对BED中的重复序列做去重处理，减少不必要的计算
• C++原生性能完全能应对50万行BED和3GB级别的hg19基因组

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