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