By liyupeng, 28 August, 2025
Forums

gff_rename.sh(论坛)

第一版

  1. 设计原理
    1. 根据https://www.kdocs.cn/l/ckSyLDAjzqyi流程化封装
    2. 功能:对输入的注释文件,先按照染色体进行分割;然后,提取基因的id、坐标位置信息,进行排序;然后,根据排序后的基因id文档,重新提取注释文件(完整内容的匹配,比如提取id123,提取的结果就是id123,不能把id123456也提取出来,这个坑踩了两次了);最后,利用convert3.sh进行重命名。
    3. 这里算是嵌套了两个脚本,同时封装了一部分流程;能够在rename版本上使用,修改后,对evm版本和pasaupdate后的版本也具有适用性
    4. 暂时没有发现其他问题,如有需要,后续更新
    5. 使用前把两个脚本添加到环境变量中
  2. 使用命令行示例
    1. i,指定输入的gff文件
    2. o,指定输出的gff文件名
    3. ID的命名,比如TgF11G06970,默认已经在代码中写好Tg,所以F是后面,用于区分雌树或者雄树,因为没有第二个物种,所以这里就只设置了这些,也够用了
bash gff_rename.sh -i Tgra.chr.convert2.gff -t F -o Tgra.chr.convert2_order.gff
  1. 代码内容
    1. gff_rename.sh
#!/bin/bash
set -euo pipefail

######################################
# 参数解析
######################################
usage() {
    echo "用法: $0 -i <input.gff> -t <tag> -o <output.gff>"
    echo "  -i 输入注释文件 (必需)"
    echo "  -t 标签 (如 F, 默认: F)"
    echo "  -o 最终输出文件 (默认: renamed_ordered.gff)"
    exit 1
}

INPUT=""
TAG="F"
OUTPUT="renamed_ordered.gff"

while getopts "i:t:o:h" opt; do
    case "$opt" in
        i) INPUT=$OPTARG ;;
        t) TAG=$OPTARG ;;
        o) OUTPUT=$OPTARG ;;
        h) usage ;;
        *) usage ;;
    esac
done

if [[ -z "$INPUT" ]]; then
    usage
fi

# 创建统一的工作目录
WORKDIR="gff_rename_dir"
mkdir -p "$WORKDIR"

echo "[1] 分割 GFF 文件到 $WORKDIR ..."
awk '{if ($1 ~ /^Chr[0-9]+$/) print > "'"$WORKDIR"'/"$1".gff"}' "$INPUT"

for i in $(seq -w 01 11); do
    echo "[2] 提取 Chr$i 基因坐标..."
    awk '$3=="gene"{
        if ($9 ~ /^ID=evm\.TU\./) {
            sub(/^ID=evm\.TU\./,"",$9)
            split($9,a,";")
            id=a[1]
        } else {
            split($9,a,";")
            sub(/^ID=/,"",a[1])
            id=a[1]
        }
        print id "\t" $4 "\t" $5 "\t" $7
    }' $WORKDIR/Chr$i.gff > $WORKDIR/chr${i}.txt

    echo "[3] 排序 Chr$i 基因..."
    { awk '$4=="+"' $WORKDIR/chr${i}.txt | sort -k2,2n -k3,3n; \
      awk '$4=="-"' $WORKDIR/chr${i}.txt | sort -k2,2n -k3,3n; } > $WORKDIR/chr${i}_sorted.txt

    echo "[4] 提取 Chr$i 排序后注释..."
    while read id start end strand; do
        grep -P "\b${id}\b" $WORKDIR/Chr$i.gff
    done < $WORKDIR/chr${i}_sorted.txt > $WORKDIR/Chr${i}_sorted.gff

    echo "[5] 调用 convert3.sh 修改 ID..."
    bash convert3.sh -i $WORKDIR/Chr${i}_sorted.gff \
        -c ${i}G -t $TAG -o $WORKDIR/${i}rename.order.gff3
done

echo "[6] 合并结果到 $WORKDIR/$OUTPUT ..."
cat $(for i in $(seq -w 01 11); do echo $WORKDIR/${i}rename.order.gff3; done) > $WORKDIR/$OUTPUT

echo "[7] 在基因之间增加空行 ..."
awk 'NR>1 && $3=="gene"{print ""} {print}' $WORKDIR/$OUTPUT > $WORKDIR/${OUTPUT}.tmp && mv $WORKDIR/${OUTPUT}.tmp $WORKDIR/$OUTPUT

echo "[8] 检查 gene 数量..."
GENE_COUNT=$(grep -cP '^\S+\s+\S+\s+gene\s' $WORKDIR/$OUTPUT)
echo "基因数量: $GENE_COUNT"

echo "✅ 全部完成!最终文件在: $WORKDIR/$OUTPUT"
  1. convert3.sh
#!/bin/bash

# 默认值
input_file=""
chromosome=""
suffix=""
output_file=""

# 解析命令行参数
while getopts "i:c:t:o:" opt; do
    case ${opt} in
        i) input_file=${OPTARG} ;;
        c) chromosome=${OPTARG} ;;
        t) suffix=${OPTARG} ;;
        o) output_file=${OPTARG} ;;
        \?)
            echo "Usage: $0 [-i input_file] [-c chromosome] [-t suffix] [-o output_file]"
            exit 1
            ;;
    esac
done

# 检查输入文件是否指定
if [ -z "$input_file" ]; then
    echo "Error: Input file not specified. Use -i option."
    exit 1
fi

# 如果没有手动指定 chromosome,通过输入文件名解析染色体编号
if [ -z "$chromosome" ]; then
    base_filename=$(basename "$input_file")
    if [[ $base_filename =~ ^Chr([0-9]+)\. ]]; then
        chrom_num=$(printf "%02d" "${BASH_REMATCH[1]}")
        chromosome="${chrom_num}G"
    else
        echo "Error: Could not extract chromosome number from file name."
        exit 1
    fi
fi

# 如果没有指定输出文件名,则默认使用输入文件名加上 "_convert.gff" 后缀
if [ -z "$output_file" ]; then
    output_file="${input_file%.gff}_convert.gff"
fi

# 清除输入文件中的注释行并处理 GFF 文件
awk '!/^#/' "$input_file" | awk -v chromosome="$chromosome" -v suffix="$suffix" '
BEGIN { gene_id = 1 }
$3 == "gene" {
    id = sprintf("Tg%s%s%04d0", suffix, chromosome, gene_id);
    sub(/ID=[^;]+/, "ID=" id);
    sub(/Name=[^;]+/, "Name=" id);
    print;
    gene_map[gene_id] = id;  # 记录基因 ID
    gene_id++;
}
$3 == "mRNA" {
    parent_id = gene_map[gene_id-1];
    id = sprintf("%s-mRNA%d", parent_id, mRNA_count[gene_id-1]+1);
    sub(/ID=[^;]+/, "ID=" id);
    sub(/Parent=[^;]+/, "Parent=" parent_id);
    sub(/Name=[^;]+/, "Name=" id);
    print;
    mRNA_map[parent_id] = id;  # 记录 mRNA ID
    mRNA_count[gene_id-1]++;
}
$3 != "gene" && $3 != "mRNA" {
    parent_id = mRNA_map[gene_map[gene_id-1]];
    sub(/ID=[^;]+/, "ID=" parent_id);
    sub(/Parent=[^;]+/, "Parent=" parent_id);
    print;
}
' | awk -F'\t' '
BEGIN { OFS="\t"; }
{
    if ($3 == "exon" || $3 == "CDS") {
        split($9, id_arr, ";");
        split(id_arr[1], id_parts_arr, "=");
        id = id_parts_arr[2];

        if (!seen[id]) {
            exon_count[id] = 0;
            cds_count[id] = 0;
            seen[id] = 1;
        }

        if ($3 == "exon") {
            exon_count[id]++;
            new_id = id ":exon" exon_count[id];
        } else if ($3 == "CDS") {
            cds_count[id]++;
            new_id = id ":CDS" cds_count[id];
        }

        sub(/ID=[^;]+/, "ID=" new_id, $9);
        print $1, $2, $3, $4, $5, $6, $7, $8, $9;
    } else {
        print $0;
    }
}
' > "$output_file"

echo "Processing complete. Output saved to $output_file"