如何用Bash提取GFF文件指定染色体的Gene行基因ID与链信息

需求说明 #

我有一个GFF文件，希望通过Bash脚本实现以下功能：

• 提取feature type为gene的行（第3列）
• 提取这些行的基因ID（第9列的ID=字段）和链信息（第7列）
• 支持按指定染色体过滤

尝试的脚本 #

#!/bin/bash
gff_file="$1"
sequence_id="$2"
gene_ids=$(zcat "$gff_file" | awk -F'\t' -v seq_id="$sequence_id" '$1 == seq_id && $3 == "gene"             {match($9, /ID=([^;]+)/, arr); print arr[1]$7}')
if [ -z "$gene_ids" ]; then
    exit 0i
sorted_gene_ids=$(echo "$gene_ids" | sort -nk2 | awk '{print $1}')
echo "$sorted_gene_ids"

GFF文件示例 #

chr1    v1.0    gene    289 3692    .   -   .   ID=Oeu061231.1;tid=PAC:37727357;id=gOeu061231.1;Name=Oeu061231.1;gene_id=Oeu061231.1
chr1    v1.0    mRNA    289 3692    .   -   .   ID=Oeu061231.1;Parent=Oeu061231.1;Name=Oeu061231.1;gene_id=Oeu061231.1
chr1    v1.0    exon    289 349 .   -   .   ID=Oeu061231.1:exon:1;Parent=Oeu061231.1;Name=Oeu061231.1;gene_id=Oeu061231.1
chr1    v1.0    CDS 289 349 .   -   1   ID=Oeu061231.1:CDS;Parent=Oeu061231.1;Name=Oeu061231.1;gene_id=Oeu061231.1
chr1    v1.0    exon    473 787 .   -   .   ID=Oeu061231.1:exon:2;Parent=Oeu061231.1;Name=Oeu061231.1;gene_id=Oeu061231.1
chr1    v1.0    CDS 473 787 .   -   1   ID=Oeu061231.1:CDS;Parent=Oeu061231.1;Name=Oeu061231.1;gene_id=Oeu061231.1
...
chr2    v1.0    gene    21189213    21190423    .   +   .   ID=Oeu046640.1;tid=PAC:37723918;id=gOeu046640.1;Name=Oeu046640.1;gene_id=Oeu046640.1
chr2    v1.0    mRNA    21189213    21190423    .   +   .   ID=Oeu046640.1;Parent=Oeu046640.1;Name=Oeu046640.1;gene_id=Oeu046640.1
chr2    v1.0    exon    21189213    21189336    .   +   .   ID=Oeu046640.1:exon:1;Parent=Oeu046640.1;Name=Oeu046640.1;gene_id=Oeu046640.1
chr2    v1.0    CDS 21189213    21189336    .   +   0   ID=Oeu046640.1:CDS;Parent=Oeu046640.1;Name=Oeu046640.1;gene_id=Oeu046640.1
chr2    v1.0    exon    21189890    21189977    .   +   .   ID=Oeu046640.1:exon:2;Parent=Oeu046640.1;Name=Oeu046640.1;gene_id=Oeu046640.1
chr2    v1.0    CDS 21189890    21189977    .   +   2   ID=Oeu046640.1:CDS;Parent=Oeu046640.1;Name=Oeu046640.1;gene_id=Oeu046640.1
chr2    v1.0    exon    21190084    21190150    .   +   .   ID=Oeu046640.1:exon:3;Parent=Oeu046640.1;Name=Oeu046640.1;gene_id=Oeu046640.1
chr2    v1.0    CDS 21190084    21190150    .   +   1   ID=Oeu046640.1:CDS;Parent=Oeu046640.1;Name=Oeu046640.1;gene_id=Oeu046640.1
chr2    v1.0    exon    21190370    21190423    .   +   .   ID=Oeu046640.1:exon:4;Parent=Oeu046640.1;Name=Oeu046640.1;gene_id=Oeu046640.1
chr2    v1.0    CDS 21190370    21190423    .   +   0   ID=Oeu046640.1:CDS;Parent=Oeu046640.1;Name=Oeu046640.1;gene_id=Oeu046640.1

脚本调用方式 #

gene_id_extracter.sh example.gff.gz "chr2"

预期输出 #

Oeu046640.1+

脚本问题修正 #

原脚本存在两处问题：

1. exit 0后有乱码字符，导致脚本无法正常执行
2. 排序逻辑无效：print arr[1]$7将基因ID和链信息直接拼接为一个字符串，后续sort -nk2无法识别$2字段，排序操作失效

修正后的脚本（无需排序） #

#!/bin/bash
gff_file="$1"
sequence_id="$2"
zcat "$gff_file" | awk -F'\t' -v seq_id="$sequence_id" '$1 == seq_id && $3 == "gene" {
    match($9, /ID=([^;]+)/, arr)
    if (arr[1] != "") print arr[1]$7
}'

修正后的脚本（需按链排序） #

如果需要按链信息排序，可先将基因ID和链信息用分隔符分开，排序后再拼接：

#!/bin/bash
gff_file="$1"
sequence_id="$2"
zcat "$gff_file" | awk -F'\t' -v seq_id="$sequence_id" '$1 == seq_id && $3 == "gene" {
    match($9, /ID=([^;]+)/, arr)
    if (arr[1] != "") print arr[1], $7
}' | sort -k2 | awk '{print $1$2}'

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