Shell脚本为GTF文件添加Entrez_ID字段失败，求排查错误原因

问题背景 #

GTF是用于基因组注释的文件格式。我编写Shell脚本，读取包含Entrez_ID字段的参考文件，匹配GTF文件中的gene_id后，将对应Entrez_ID添加到GTF条目末尾，但实际输出中添加的Entrez_ID值是字段名而非参考文件中的数值。

GTF文件示例 #

cat Canis_lupus_familiaris.ROS_Cfam_1.0.108.gtf | head
#!genome-build ROS_Cfam_1.0
#!genome-version ROS_Cfam_1.0
#!genome-date 2020-09
#!genome-build-accession GCA_014441545.1
#!genebuild-last-updated 2020-10
X       ensembl gene    24550462        24552226        .       -       .       gene_id "ENSCAFG00845015183"; gene_version "1"; gene_source "ensembl"; gene_biotype "protein_coding";
X       ensembl transcript      24550462        24552226        .       -       .       gene_id "ENSCAFG00845015183"; gene_version "1"; transcript_id "ENSCAFT00845027108"; transcript_version "1"; gene_source "ensembl"; gene_biotype "protein_coding"; transcript_source "ensembl"; transcript_biotype "protein_coding"; tag "Ensembl_canonical";
X       ensembl exon    24552206        24552226        .       -       .       gene_id "ENSCAFG00845015183"; gene_version "1"; transcript_id "ENSCAFT00845027108"; transcript_version "1"; exon_number "1"; gene_source "ensembl"; gene_biotype "protein_coding"; transcript_source "ensembl"; transcript_biotype "protein_coding"; exon_id "ENSCAFE00845128634"; exon_version "1"; tag "Ensembl_canonical";
X       ensembl CDS     24552206        24552226        .       -       0       gene_id "ENSCAFG00845015183"; gene_version "1"; transcript_id "ENSCAFT00845027108"; transcript_version "1"; exon_number "1"; gene_source "ensembl"; gene_biotype "protein_coding"; transcript_source "ensembl"; transcript_biotype "protein_coding"; protein_id "ENSCAFP00845021332"; protein_version "1"; tag "Ensembl_canonical";
X       ensembl start_codon     24552224        24552226        .       -       0       gene_id "ENSCAFG00845015183"; gene_version "1"; transcript_id "ENSCAFT00845027108"; transcript_version "1"; exon_number "1"; gene_source "ensembl"; gene_biotype "protein_coding"; transcript_source "ensembl"; transcript_biotype "protein_coding"; tag "Ensembl_canonical";

参考文件（含Entrez_ID） #

gene_id Entrez_ID
ENSCAFG00845006432 399518
ENSCAFG00845002136 399530
ENSCAFG00845029798 399544
ENSCAFG00845011460 399545
ENSCAFG00845001610 399653
ENSCAFG00845013158 403157
ENSCAFG00845014982 403168
ENSCAFG00845021967 403170
ENSCAFG00845019241 403400

期望输出示例 #

#!genome-build ROS_Cfam_1.0
#!genome-version ROS_Cfam_1.0
#!genome-date 2020-09
#!genome-build-accession GCA_014441545.1
#!genebuild-last-updated 2020-10
X       ensembl gene    24550462        24552226        .       -       .       gene_id "ENSCAFG00845006432"; gene_version "1"; gene_source "ensembl"; gene_biotype "protein_coding"; Entrez_ID "399518";
X       ensembl transcript      24550462        24552226        .       -       .       gene_id "ENSCAFG00845006432"; gene_version "1"; transcript_id "ENSCAFT00845027108"; transcript_version "1"; gene_source "ensembl"; gene_biotype "protein_coding"; transcript_source "ensembl"; transcript_biotype "protein_coding"; tag "Ensembl_canonical"; Entrez_ID "399518";
X       ensembl exon    24552206        24552226        .       -       .       gene_id "ENSCAFG00845006432"; gene_version "1"; transcript_id "ENSCAFT00845027108"; transcript_version "1"; exon_number "1"; gene_source "ensembl"; gene_biotype "protein_coding"; transcript_source "ensembl"; transcript_biotype "protein_coding"; exon_id "ENSCAFE00845128634"; exon_version "1"; tag "Ensembl_canonical"; Entrez_ID "399518";
X       ensembl CDS     24552206        24552226        .       -       0       gene_id "ENSCAFG00845006432"; gene_version "1"; transcript_id "ENSCAFT00845027108"; transcript_version "1"; exon_number "1"; gene_source "ensembl"; gene_biotype "protein_coding"; transcript_source "ensembl"; transcript_biotype "protein_coding"; protein_id "ENSCAFP00845021332"; protein_version "1"; tag "Ensembl_canonical"; Entrez_ID "399518";
X       ensembl start_codon     24552224        24552226        .       -       0       gene_id "ENSCAFG00845006432"; gene_version "1"; transcript_id "ENSCAFT00845027108"; transcript_version "1"; exon_number "1"; gene_source "ensembl"; gene_biotype "protein_coding"; transcript_source "ensembl"; transcript_biotype "protein_coding"; tag "Ensembl_canonical"; Entrez_ID "399518";

编写的脚本 #

#!/bin/bash

input_file="$1"
gff_file="$2"
output_file="$3"

# Create a temporary file for storing the intermediate results
temp_file=$(mktemp)

# Function to append temporary file contents to output file
cleanup() {
    if [[ -f "$temp_file" ]]; then
        cat "$temp_file" >> "$output_file"
        rm "$temp_file"
    fi
}

# Trap termination signals and execute cleanup function
trap cleanup EXIT TERM

# Read the input file line by line
while IFS=$'\t' read -r gene_id Entrez_ID; do
    # Search for the corresponding line in the GFF file and append the annotation
    awk -v id="$gene_id" -v eid="$Entrez_ID" 'BEGIN {OFS="\t"} $0 ~ id { sub(/;$/, "; Entrez_ID \"" eid "\";", $0) } 1' "$gff_file" >> "$temp_file"
done < "$input_file"

# Combine the original GFF file with the annotated lines and store the result in the output file
cat "$temp_file" > "$output_file"

# Clean up the temporary file
cleanup

echo "Mapping and annotation complete. Output file: $output_file"

实际输出结果 #

cat output.gtf | head
#!genome-build ROS_Cfam_1.0
#!genome-version ROS_Cfam_1.0
#!genome-date 2020-09
#!genome-build-accession GCA_014441545.1
#!genebuild-last-updated 2020-10
X       ensembl gene    24550462        24552226        .       -       .       gene_id "ENSCAFG00845015183"; gene_version "1"; gene_source "ensembl"; gene_biotype "protein_coding"; Entrez_ID "Entrez_ID";
X       ensembl transcript      24550462        24552226        .       -       .       gene_id "ENSCAFG00845015183"; gene_version "1"; transcript_id "ENSCAFT00845027108"; transcript_version "1"; gene_source "ensembl"; gene_biotype "protein_coding"; transcript_source "ensembl"; transcript_biotype "protein_coding"; tag "Ensembl_canonical"; Entrez_ID "Entrez_ID";
X       ensembl exon    24552206        24552226        .       -       .       gene_id "ENSCAFG00845015183"; gene_version "1"; transcript_id "ENSCAFT00845027108"; transcript_version "1"; exon_number "1"; gene_source "ensembl"; gene_biotype "protein_coding"; transcript_source "ensembl"; transcript_biotype "protein_coding"; exon_id "ENSCAFE00845128634"; exon_version "1"; tag "Ensembl_canonical"; Entrez_ID "Entrez_ID";
X       ensembl CDS     24552206        24552226        .       -       0       gene_id "ENSCAFG00845015183"; gene_version "1"; transcript_id "ENSCAFT00845027108"; transcript_version "1"; exon_number "1"; gene_source "ensembl"; gene_biotype "protein_coding"; transcript_source "ensembl"; transcript_biotype "protein_coding"; protein_id "ENSCAFP00845021332"; protein_version "1"; tag "Ensembl_canonical"; Entrez_ID "Entrez_ID";
X       ensembl start_codon     24552224        24552226        .       -       0       gene_id "ENSCAFG00845015183"; gene_version "1"; transcript_id "ENSCAFT00845027108"; transcript_version "1"; exon_number "1"; gene_source "ensembl"; gene_biotype "protein_coding"; transcript_source "ensembl"; transcript_biotype "protein_coding"; tag "Ensembl_canonical"; Entrez_ID "Entrez_ID";

错误原因与修正方案 #

核心错误点 #

1. 未跳过参考文件表头：脚本读取参考文件时，第一行表头gene_id Entrez_ID被当作有效数据处理，导致Entrez_ID变量被赋值为字符串"Entrez_ID"，最终写入GTF文件。
2. 重复处理GTF文件：每读取一行参考文件就完整遍历一次GTF文件并追加到临时文件，最终输出会包含GTF文件的多份拷贝，完全不符合需求。
3. 匹配逻辑不严谨：用$0 ~ id会匹配任何包含该gene_id的字符串，可能误匹配到transcript_id、exon_id等其他字段。

修正后的脚本 #

#!/bin/bash

input_file="$1"
gff_file="$2"
output_file="$3"

# 使用awk一次性完成映射与替换，提升效率
awk '
BEGIN {
    OFS = "\t"
    # 读取参考文件构建gene_id到Entrez_ID的映射表
    while ((getline line < "'"$input_file"'") > 0) {
        if (NR == 1) continue;  # 跳过表头
        split(line, arr, /\t/);
        map[arr[1]] = arr[2];
    }
    close("'"$input_file"'");
}
# 处理GTF文件
{
    if ($0 ~ /^#/) {  # 保留注释行
        print;
        next;
    }
    # 精准提取gene_id字段
    if (match($0, /gene_id "([^"]+)"/, gene)) {
        gene_id = gene[1];
        if (gene_id in map) {
            # 在行末尾添加Entrez_ID字段
            sub(/;?$/, "; Entrez_ID \"" map[gene_id] "\";", $0);
        }
    }
    print;
}
' "$gff_file" > "$output_file"

echo "Mapping and annotation complete. Output file: $output_file"

修正说明 #

1. 跳过表头：读取参考文件时通过NR == 1跳过第一行表头，避免无效数据干扰。
2. 一次性处理：用awk的getline读取参考文件构建映射表，然后一次性遍历GTF文件完成替换，避免重复处理导致的输出冗余。
3. 精准匹配：用正则表达式gene_id "([^"]+)"提取gene_id，确保只匹配目标字段，避免误匹配。
4. 简化流程：去掉不必要的临时文件和清理逻辑，直接将结果写入输出文件，提升脚本效率。

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