By liyupeng, 28 August, 2025
Forums
Researches

注释文件id重命名v4

情况介绍

之前的注释，比如NCC和BA，虽然改了名字，但是这些排序上相邻的基因在基因组上的坐标不是相邻的，这样可能有些不太美观，主要是没法看duplication的基因，观察这些规律等等。因此，v4将重新对注释文件中的基因按照坐标进行排序，然后重新命名。至于已有的NCC和BA版本的注释文件，考虑到已经流传开来，这里将制作新的对应表。注释的内容实际上没有改变，只是ID名字再次更新。

ID名稍微做一些调整，具体就是这样，就是增加一个-

NCC版注释：ID=Tgra01G00010

NCC排序版注释：ID=Tgra01G-00010

操作流程（NCC）

  1. 分割NCC的注释文件成11个染色体的注释，一个contig的注释
awk '{if ($1 ~ /^Chr[0-9]+$/) print > $1".gff"}'  Tgra.chr.convert2.gff
awk '{if ($1 !~ /^Chr[0-9]+$/) print}'  Tgra.chr.convert2.gff > ptg.gff
  1. 提取基因ID和起始、终止位置和正反链4列内容
    1. 这里只放一个染色体
awk '$3=="gene"{split($9,a,";"); id=a[1]; sub(/^ID=/,"",id); print id "\t" $4 "\t" $5 "\t" $7}' Chr1.gff > chr1.txt
  1. 按照坐标对基因进行排序
    1. 把内容拆分成正义链和反义链，正义链放在上半部分，反义链放在下半部分
    2. 按照起始和终止位置，对基因进行排序，递增（起始和终止）
{ awk '$4=="+"' gene_coords.txt | sort -k2,2n -k3,3n; awk '$4=="-"' chr1.txt | sort -k2,2n -k3,3n; } > chr1_sorted.txt
  1. 根据先前排序号的基因id，提取注释文件
    1. 按照拍过顺序的基因id提取得到的注释文件也是排序好的
while read id start end strand; do
    grep "$id" Chr1.gff
done < chr1_sorted.txt > Chr1_sorted.gff
  1. 修改注释文件的id
    1. 脚本convert2.sh，https://www.kdocs.cn/l/cdJRauELhrfX用到过的，算是ncc版特制的
    2. 示例：Tgra01G-00020，前缀Tgra已经写好，这里需要输入的只有01G-（染色体编号和后缀-）部分，至于基因ID编号则是按照文本从上到下自动命名
bash convert2.sh -i Chr1_sorted.gff -c 01G-
  1. 单独处理ptg.gff
    1. 提取id，这一步也需要改变
      1. contig部分的序列太多，这里的命名方式以contig的编号为主，代替原先的基因id
awk '$3=="gene"{print $1 "\t" $4 "\t" $5 "\t" $7}' ptg.gff > ptg.txt
  1. 排序，这里按照序列片段的编号进行排序，然后再按照坐标进行排序
    1. 先按照正义链和反义链分为两部分，正义链的在上，反义链的在下
    2. 序列编号有两个数字部分，以_分为三部分，第一部分的数字夹杂在ptg和l中间，作为第一顺序依据；第三部分的数字则是作为第二顺序排序依据。
{
  awk '
  {
    split($1, parts, "_");
    # 提取 ptg 和 l 之间数字
    if (match(parts[1], /^ptg([0-9]+)l$/, arr)) {
      prefix_num = arr[1] + 0;  # 转数字，去前导0
    } else {
      prefix_num = 0;
    }
    fragment_num = parts[3] + 0;  # fragment 后数字
    if ($4 == "+") {
      print prefix_num "\t" fragment_num "\t" $0
    }
  }' ptg.txt | sort -k1,1n -k2,2n | cut -f3-
  
  awk '
  {
    split($1, parts, "_");
    if (match(parts[1], /^ptg([0-9]+)l$/, arr)) {
      prefix_num = arr[1] + 0;
    } else {
      prefix_num = 0;
    }
    fragment_num = parts[3] + 0;
    if ($4 == "-") {
      print prefix_num "\t" fragment_num "\t" $0
    }
  }' ptg.txt | sort -k1,1n -k2,2n | cut -f3-
} > ptg_sorted.txt
  1. 序列id去重
    1. 先前忘记考虑一条序列上可能存在多个基因，以至于序列id重复出现，这导致在后面的提取注释内容的时候，得到的基因数量也随之翻倍
    2. 这里去重后保留下来的所有序列，一共331条，ptg.gff的基因数量为620
awk '!seen[$1]++' ptg_sorted.txt > ptg_sorted_unique.txt
  1. 根据排序后的序列编号提取注释内容
    1. 上面的没法用，让G老师重新写了一个
    2. 根据ptg_sorted.txt第一列，序列编号进行查询；序列编号需要完全匹配，而不是部分匹配。
while read id start end strand; do
  grep -E "^${id}[[:space:]]" ptg.gff
done < ptg_sorted_unique.txt > ptg_sorted.gff
  1. id重命名
bash convert2.sh -i ptg_sorted.gff -c _ptg-
  1. 合并所有的注释文件
cat Chr1_sorted_convert.gff \
    Chr2_sorted_convert.gff \
    Chr3_sorted_convert.gff \
    Chr4_sorted_convert.gff \
    Chr5_sorted_convert.gff \
    Chr6_sorted_convert.gff \
    Chr7_sorted_convert.gff \
    Chr8_sorted_convert.gff \
    Chr9_sorted_convert.gff \
    Chr10_sorted_convert.gff \
    Chr11_sorted_convert.gff \
    ptg_sorted_convert.gff \
    > Tgra.chr.convert2_order.gff
  1. 检查一下数量，看看有没有漏掉的
    1. 32056，数量ok
grep -cP '^\S+\s+\S+\s+gene\s' Tgra.chr.convert2_order.gff
  1. 制作NCC版本和NCC排序版的基因id对应表
python compare_gff.py -i Tgra.chr.convert2.gff  -l Tgra.chr.convert2_order.gff -o ncc_order_comparison.tsv
  1. 每个gene之间插入一个空行，这样视觉上更好看一下，每个基因也分的更加清楚，对注释内容本身没有影响
awk 'NR>1 && $3=="gene"{print ""} {print}' Tgra.chr.convert2_order.gff > female_ncc_order.gff3

附件

ncc_order注释路径

/data2/liyupeng/alice/output/female_anno/female_ncc_order_gff

ba_order注释路径

/data2/liyupeng/alice/output/female_anno/female_mark1_rename_order_gff

结果图