如何快速从GTDB代表基因组匹配目标蛋白对应CDS序列？

问题背景 #

需要从GTDB代表基因组的fna.gz文件中，匹配自定义Hit list（蛋白名称格式为NCBI基因组登录号+自然数，如AE017199.1_1、AE017199.1_2）对应的CDS序列。由于这类ID无法直接通过GTDB/NCBI接口获取对应序列，已编写bash脚本实现匹配，但单线程处理所有细菌文件耗时约6个月，急需优化效率。

原脚本核心逻辑片段：

for gz_file in "$zip_directory"/*.fna.gz; do
    current_file=$((current_file + 1))
    echo "Processing file $current_file of $total_files"
    gzip -dc "$gz_file" | while IFS= read -r line; do
        if [[ $line == ">"* ]]; then
            for word in "${words[@]}"; do
                if [[ $line == *"$word "* ]]; then
                    echo "$line - $(basename "$gz_file" .fna.gz)" >> "$temp_dir/$(basename "$gz_file" .fna.gz).temp"
                    while IFS= read -r next_line && [[ ! $next_line == ">"* ]]; do
                        echo "$next_line" >> "$temp_dir/$(basename "$gz_file" .fna.gz).temp"
                    done
                    break
                fi
            done
        fi
    done
done

原脚本的性能瓶颈 #

• 单线程串行处理：逐个文件处理，完全未利用多核CPU资源
• 低效匹配逻辑：每个序列头都要遍历整个Hit list数组，时间复杂度为O(N*M)（N为序列数，M为Hit list条目数）
• 冗余文件IO：每个文件生成单独临时文件，频繁磁盘写入拖慢速度
• bash逐行处理缺陷：bash本身不适合大规模文本处理，逐行读取解析效率极低

优化方案 #

1. 并行化处理（快速见效） #

利用GNU Parallel或xargs实现多文件并行处理，直接将处理速度提升至接近CPU核心数的倍数：

GNU Parallel示例

# 先将Hit list保存为hit_list.txt（每行一个ID）
# 定义单个文件处理函数
process_file() {
    gz_file=$1
    base_name=$(basename "$gz_file" .fna.gz)
    gzip -dc "$gz_file" | seqkit grep -f hit_list.txt -r > "results/${base_name}_matches.fasta"
}
export -f process_file

# -j指定并行数（建议设为CPU核心数的70-80%，避免资源耗尽）
find "$zip_directory" -name "*.fna.gz" | parallel -j 16 process_file {}

xargs示例

find "$zip_directory" -name "*.fna.gz" | xargs -P 16 -I {} bash -c '
    base_name=$(basename "{}" .fna.gz)
    gzip -dc "{}" | seqkit grep -f hit_list.txt -r > "results/${base_name}_matches.fasta"
'

2. 用专业序列工具替换bash逐行处理 #

使用seqkit（高性能FASTA/Q处理工具）替代bash的逐行解析，底层基于Go实现，性能远超bash脚本：

• -f hit_list.txt：指定匹配的ID列表文件
• -r：允许序列头中包含目标ID（适配你的ID是序列头一部分的场景）

可通过conda快速安装：conda install -c bioconda seqkit

3. 用awk优化匹配逻辑（无额外工具依赖） #

如果无法安装第三方工具，用awk实现哈希匹配（避免遍历数组），性能远优于bash循环：

# 将Hit list转为awk关联数组，实现O(1)查找
awk -v hit_list="hit_list.txt" '
BEGIN {
    while ((getline line < hit_list) > 0) {
        hits[line] = 1
    }
    close(hit_list)
}
/^>/ {
    # 从序列头提取目标ID（假设ID是空格分隔的第一个字段，可根据实际格式调整）
    split(substr($0, 2), parts, " ")
    target_id = parts[1]
    in_target = (target_id in hits)
    if (in_target) {
        print $0 " - " ARGV[1]
    }
    next
}
in_target { print }
' <(gzip -dc "$gz_file") >> results/all_matches.fasta

可结合GNU Parallel让该脚本并行处理多个文件。

4. 减少IO开销 #

• 避免生成大量临时文件，直接将所有匹配结果写入统一输出文件（或按分类合并）
• 使用SSD存储处理文件和结果，大幅提升磁盘读写速度

预期效果 #

通过并行化+专业工具的组合，处理时间可从6个月压缩至数天甚至数小时（取决于CPU核心数和存储性能）。

内容的提问来源于stack exchange，提问作者Marburg-Researcher