By liyupeng, 31 March, 2026
Forums
Researches

香榧雌树基因组注释第四版（TGA_V4_VA）

基本介绍

  1. 这个版本的注释命名为VA（Vesta Argo），名字起源于《在地下城寻求邂逅是否搞错了什么恶》贝尔·克朗你的招式【圣火英烈斩】，或者翻译成【英雄圣火】，原理为双重蓄力，整个适合这个将BA和NCC结合起来的新的注释版本VA
  2. VA版本的注释很简单，基本知识将BA和NCC融合起来，主体是BA的部分https://www.kdocs.cn/l/caE41kwA97SB。至于融合的流程，这设计到EVM的基本原理，市面上的注释流程普遍采取一套综合性的方法，及结合从头注释（模型预测）、同源注释（蛋白序列比对）和RNA-seq比对这三种来源的结果文件，最后用EVM根据不同来源的注释文件以及可信度设置权重（有点类似分红），整合多种渠道的注释文件，最后整合成一个结果文件
    1. 这么说肯定很难理解，打个比喻好了，地三鲜是用茄子、土豆和青椒三种食材组成的，把地三鲜当做我们常用的注释文件（结果文件）的话，那么茄子、土豆和青椒在这里就是三种注释方法的过程中间文件，直接吃（拿来用）也可以，但没有做成地三鲜后好吃，差不多是这么一回事

前期准备

  1. 修改ncc版本注释文件第二列，将EVM来源更改为ncc来源，为了以防万一（因为之类是打算用evm进行整合的，避免可能存在的冲突）
awk 'BEGIN{OFS="\t"} {$2="ncc"; print}' Tgra.chr.convert2b.gff > Tgra.chr.convert2b.ncc.gff
  1. 拆分NCC版本的注释文件
    1. 因为我是分染色体进行跑流程的，这样比较快
    2. 还有，这里contig的注释内容被丢掉了
      1. contig的序列只有161mb，比较小
      2. contig序列的注释内容不怎么关注
      3. 因为序列少，存在的基因上，BA中没法用braker进行注释，所以也没有这部分的内容，当然在VA当中是不可能存在的，今后倒是有考虑在MA里面放进去，不过那是几个月后的事情了，暂且不在这里过多赘述
awk '{if ($1 ~ /^Chr[0-9]+$/) print > $1".gff3"}' Tgra.chr.convert2b.ncc.gff

EVM整合

这里使用脚本evm.sh来完成所有整合步骤：https://www.kdocs.cn/l/cay1X2HyUNvI

参数解读：

  1. b，braker的注释文件，gtf格式
  2. d，RNA注释步骤最后一步（提取开放阅读框）得到的gff3
  3. g，同源注释得到的gff3
  4. p，RNA注释步骤倒数第二步得到的注释文件，转录本比对到参考基因组上的注释信息
  5. f，参考基因组文件
  6. n，指定染色体编号，用于修改weights文件用
  7. t，线程信息，这里指定是10
  8. 补充：这里还需要将weights_base.txt放置在容器目录中（因为脚本中这么写的，好像是/home目录下，可以更改，考虑在后续版本增加参数，指定这个文件，现在是属于默认状态，没有的话直接报错），这是权重文件，已经提前写好了一部分，剩余的部分脚本在运行的时候会自动编入，根据不同的染色体（这个是有RNA注释的gff3输入文件决定）。
  9. 输出文件，比如这里是femalee的1号染色体，所以输出文件时：01evm.out.gff3
nohup bash evm.sh \
  -b 01braker.gtf \
  -d 01.assemblies.fasta.transdecoder.genome.gff3 \
  -g 01gth_output.gff \
  -p 01.pasa_assemblies.gff3 \
  -f Chr1.fna \
  -n 01 \
  -t 10 &
  1. weights_base.txt
    1. 因为ncc是筛选过的高质量的基因，所以在这里赋予它和RNA-seq比对的注释一样的最高权重10
ABINITIO_PREDICTION   braker           1
PROTEIN               genomeThreader   3
OTHER_PREDICTION      transdecoder     5
TRANSCRIPT            assembler-01    10
OTHER_PREDICTION      ncc             10
  1. 注释文件预处理（脚本路径：/home/soft/EVidenceModeler/EvmUtils/misc/）
    1. 格式转换：将其他soft生成的gff/gtf转换成evm支持输入格式的gff3
braker_GTF_to_EVM_GFF3.pl 01braker.gtf > denovo.gff3
genomeThreader_to_evm_gff3.pl 01gth_output.gff > evm_pro.gff3
cat 01.pasa_assemblies.gff3 > evm_rna.gff3
  1. gff3文件的处理与合并
    1. 修改修改braker来源的注释文件，将其中三个来源的全部换成braker来源
      1. 这样可以减少权重分配，保证每个软件预测得到的基因都被保留下来，最后整合过的基因数量不会过低
    2. transdecoder是基于转录本的序列特征来预测开放阅读框，并从这些预测的ORFs中推断基因编码区域，并非蛋白/RNA-seq的比对，所以归类到gene prediction中
awk 'BEGIN{OFS="\t"} {if(NF > 0) $2="braker"; print}' denovo.gff3 > denovo2.gff3
cat denovo2.gff3 01.assemblies.fasta.transdecoder.genome.gff3 Chr1.gff3 > evm_denovo.gff3
  1. 分割原始数据，用于后续并行运行
    1. 分区大小（segmentsize）一般在10mb-50mb之间
    2. 重叠区域（overlapsize）一般是分区大小的10%
      1. 参数：分区20mb，重叠2mb
partition_EVM_inputs.pl \
  --genome Chr1.fna \
  --gene_predictions evm_denovo.gff3 \
  --protein_alignments evm_pro.gff3 \
  --transcript_alignments evm_rna.gff3 \
  --segmentSize 12000000 --overlapSize 1200000 \
  --partition_listing partitions_list.out \
  --partition_dir partitions_output
  1. 创建并行命令并执行
write_EVM_commands.pl --genome Chr1.fna \
    --weights /home/server/work/va/01/weights.txt \
    --gene_predictions evm_denovo.gff3 \
    --protein_alignments evm_pro.gff3 \
    --transcript_alignments evm_rna.gff3 \
    --output_file_name evm.out \
    --partitions partitions_list.out > commands.list
nohup parallel --jobs 30 < commands.list &
  1. 合并并行结果
recombine_EVM_partial_outputs.pl \
  --partitions partitions_list.out \
  --output_file_name evm.out
  1. 格式转换
convert_EVM_outputs_to_GFF3.pl --partitions partitions_list.out --output_file_name evm.out --genome Chr1.fna
  1. 用于检测生成文件的基因数量的命令（根据第3列gene出现次数进行判定）
    1. 基因数量：45230
grep -cP '^\S+\s+\S+\s+gene\s'

 

PASA优化

后期加工（重命名、筛选、提取）

  1. 更改染色体编号的格式，Chr1-Chr9改成Chr01-Chr09这种格式，其他内容不变
awk 'BEGIN{OFS="\t"} $1 ~ /^Chr[1-9]$/ {$1 = "Chr0" substr($1, 4)} {print}' female_mark2.gff3 > female_mark2.renamed.gff3
  1. 清除注释文件中的#注释行
grep -v '^#' female_mark2.renamed.gff3 > female_mark2.no_comment.gff3
  1. 根据染色体编号，将注释文件拆分成11个
awk '$1 ~ /^Chr[0-9]+$/ {print > $1".gff3"}' female_mark2.no_comment.gff3
  1. 修改基因ID
    1. 这里再参数c的内容中加入v4，代表第4版，这样可以用于和BA版作区分，不然两个版本的基因ID命名格式一样，容易弄混淆
bash convert3.sh -i Chr01.gff3 -c 01G -t Fv4 -o 01evm.update_rename.gff3
bash convert3.sh -i Chr02.gff3 -c 02G -t Fv4 -o 02evm.update_rename.gff3
bash convert3.sh -i Chr03.gff3 -c 03G -t Fv4 -o 03evm.update_rename.gff3
bash convert3.sh -i Chr04.gff3 -c 04G -t Fv4 -o 04evm.update_rename.gff3
bash convert3.sh -i Chr05.gff3 -c 05G -t Fv4 -o 05evm.update_rename.gff3
bash convert3.sh -i Chr06.gff3 -c 06G -t Fv4 -o 06evm.update_rename.gff3
bash convert3.sh -i Chr07.gff3 -c 07G -t Fv4 -o 07evm.update_rename.gff3
bash convert3.sh -i Chr08.gff3 -c 08G -t Fv4 -o 08evm.update_rename.gff3
bash convert3.sh -i Chr09.gff3 -c 09G -t Fv4 -o 09evm.update_rename.gff3
bash convert3.sh -i Chr10.gff3 -c 10G -t Fv4 -o 10evm.update_rename.gff3
bash convert3.sh -i Chr11.gff3 -c 11G -t Fv4 -o 11evm.update_rename.gff3
  1. 合并11个修改后的注释文件
cat 01evm.update_rename.gff3 02evm.update_rename.gff3 03evm.update_rename.gff3 04evm.update_rename.gff3 05evm.update_rename.gff3 06evm.update_rename.gff3 07evm.update_rename.gff3 08evm.update_rename.gff3 09evm.update_rename.gff3 10evm.update_rename.gff3 11evm.update_rename.gff3 > merged.gff3
  1. 插入空行（因为看起来更好）
    1. 把每个基因的注释内容分开来，这样看起来更舒服一些
awk '{if ($3=="gene") print ""; print}' merged.gff3 > female_mark2_rename.gff3
  1. 提取cds/蛋白序列
gffread female_mark2_rename.gff3 -g female_chr_rename.fna -x cds.fa
gffread female_mark2_rename.gff3 -g female_chr_rename.fna -y protein.fa
  1. 用脚本cds.check.py分离结构完整的和不完整的cds序列
    1. complete.cds.fa：64667
    2. incomplete.cds.fa：86
python cds.check.py -i cds.fa
  1. 提取完整cds序列的id
grep '^>' complete.cds.fa | sed 's/^>//' > complete.cds.id.txt
  1. 过滤蛋白序列
    1. 长度<50的过滤掉
    2. 存在未知符号“.”的过滤掉
    3. qualified_protein：64724
    4. unqualified_protein：29
python protein_select.py -i protein.fa
  1. 提取合格的蛋白序列的id
grep '^>' qualified_protein.fa | sed 's/^>//' > qualified_protein.id.txt
  1. 取合格的cds的id和蛋白的id的交集，这部分id将用作提取最终的cds/蛋白序列
    1. 交集id数量：64639
comm -12 <(sort qualified_protein.id.txt) <(sort complete.cds.id.txt) > id_intersection.txt
  1. 根据交集id文件提取cds和蛋白序列
seqkit grep -f id_intersection.txt complete.cds.fa > female_mark2_cds.fa
seqkit grep -f id_intersection.txt protein.fa > female_mark2_protein.fa
  1. 重新运行一遍两个脚本，检查无误

  1. 筛选掉的基因，留个文档
    1. 拿个所有蛋白序列的id文档（64753）
grep '^>' protein.fa | sed 's/^>//' > protein.id.txt
  1. 把所有蛋白序列的id文档和最终蛋白/cds序列id文档比对，取差集
    1. 114
comm -23 <(sort protein.id.txt) <(sort id_intersection.txt) > dropped_mrna_ids.txt
  1. 文件存储路径，VA版注释、提取蛋白、cds and so on
/data2/liyupeng/alice/output/female_anno/va_anno/rename_gff3_extract

补丁（2026/3）：基因排序，然后重命名

  1. 分割VA的注释文件成11个染色体的注释
awk '{if ($1 ~ /^Chr[0-9]+$/) print > $1".gff"}'  female_mark2_rename.gff3
  1. 提取基因ID和起始、终止位置和正反链4列内容
    1. 这里只放一个染色体
awk '$3=="gene"{split($9,a,";"); id=a[1]; sub(/^ID=/,"",id); print id "\t" $4 "\t" $5 "\t" $7}' Chr01.gff > chr01.txt
  1. 按照坐标对基因进行排序
    1. 按照起始和终止位置，对基因进行排序，递增（起始和终止）
sort -k2,2n -k3,3n chr01.txt > chr01_sorted.txt
  1. 根据先前排序号的基因id，提取注释文件
    1. 按照排过顺序的基因id提取得到的注释文件也是排序好的
while read id start end strand; do
    grep "$id" Chr01.gff
done < chr1_sorted.txt > Chr01_sorted.gff
  1. 每个gene之间插入一个空行，这样视觉上更好看一下，每个基因也分的更加清楚，对注释内容本身没有影响
awk 'NR>1 && $3=="gene"{print ""} {print}' female_ba.gff3 > female_ba_last.gff3
  1. 修改基因ID
    1. cds不计数
    2. 前缀改为TGA
    3. 插入空行命令行也加进去
    4. 添加5端和3端的注释内容，就是区别分开，用于基于jbrower的可视化框架使用
bash convert4.sh -i Chr01_sorted.gff -c 01G -t F -o Chr01_rename.gff3
bash convert4.sh -i Chr02_sorted.gff -c 02G -t F -o Chr02_rename.gff3
bash convert4.sh -i Chr03_sorted.gff -c 03G -t F -o Chr03_rename.gff3
bash convert4.sh -i Chr04_sorted.gff -c 04G -t F -o Chr04_rename.gff3
bash convert4.sh -i Chr05_sorted.gff -c 05G -t F -o Chr05_rename.gff3
bash convert4.sh -i Chr06_sorted.gff -c 06G -t F -o Chr06_rename.gff3
bash convert4.sh -i Chr07_sorted.gff -c 07G -t F -o Chr07_rename.gff3
bash convert4.sh -i Chr08_sorted.gff -c 08G -t F -o Chr08_rename.gff3
bash convert4.sh -i Chr09_sorted.gff -c 09G -t F -o Chr09_rename.gff3
bash convert4.sh -i Chr10_sorted.gff -c 10G -t F -o Chr10_rename.gff3
bash convert4.sh -i Chr11_sorted.gff -c 11G -t F -o Chr11_rename.gff3
  1. 合并11个修改后的注释文件
cat Chr01_rename.gff3 Chr02_rename.gff3 Chr03_rename.gff3 Chr04_rename.gff3 Chr05_rename.gff3 Chr06_rename.gff3 Chr07_rename.gff3 Chr08_rename.gff3 Chr09_rename.gff3 Chr10_rename.gff3 Chr11_rename.gff3 > female_mark4.gff3
  1. 提取cds/蛋白序列
gffread female_mark4.gff3 -g female_chr_rename.fna -x female_mark4.cds.fa
gffread female_mark4.gff3 -g female_chr_rename.fna -y female_mark4.protein.fa

VA版本的评估

  1. 首先，基因数量是45084，介于NC和BA之间。真核生物的基因数量在物种间的变动不大，好像是4万+左右，这个数量也正好合格
  2. cds检查，不完整的cds（缺少起始或终止密码子的）只有86个，这个数值也在于其之内
  3. 接下来就是通过IGV可视化进行的评估
  4. 部分基因在NC（绿色）和BA（红色）中都没有，合并后的VA（红色）中也不存在，这里只能通过IGV添加进去，金色的部分就是正在编辑的注释文件（ZA）

  1. 一部分基因，NC注释到了一部分，但没有被整合进去。可能是因为这部分在ncc因为cds结构不完整被去掉了，也有可能是因为可信度较低被算法过滤了，这部分也只能通过手动添加

  1. 一部分nc和ba的基因结构重叠，最后是以BA的结构为主的；这部分基因，BA是正确的，NC虽然有注释到这个基因位点的部分内容，但过长或者过短

  1. 比如NC把一个基因位点分割成了两个不同的基因，这在长内含子基因中挺常见的，而BA则是正确的注释到了。helixer也注释到了基因位点的部分内容，不过在对预测这类长内含子的效果不是很好，能够预测到基因，但不完整

  1. VA版本的注释也存在部分这样的问题，通过计算，只获取了BA和NC的重叠部分（可能是因为在EVM整合时给的权重过高了，最高级），以至于整合后的VA基因就是NC的注释内容，但实际可视化后可以看到这个基因不完整（不正确），但BA是正确注释的，这就是问题所在，这里也需要手动调整

  1. 上面的是NC注释基因过短的，下面这个是NC注释基因过长的，从覆盖度、剪切位点各个方面来看，NC中的这个基因多了一个外显子，以至于产生一条长的内含子。

  1. nc和ba都不行的，合并后的va也注释不全的。这里，nc和ba都只注释到了一部分的基因结构，融合后的VA也是一样，得手动补充左边半部分的基因结构内容。

  1. VA注释中，有的基因位点上，存在两个基因的重叠的问题。这里因为存在一个正确的基因（蓝色），所以把金色的删除

  1. 总的来说，VA虽然补充了NC的数量优势和BA的质量优势，但部分基因结构还是需要调整。不过，使用的话应该没问题，基因顺序排列过了，重名命名过了，也能够通过jbrower的可视化框架结构。

附件

protein_selcet.py

#!/usr/bin/env python3

import argparse
from Bio import SeqIO

def parse_args():
    parser = argparse.ArgumentParser(
        description="筛选蛋白序列长度大于等于设定值的序列（默认为50），"
                    "默认输出文件为 qualified_protein.fa 和 unqualified_protein.fa",
        formatter_class=argparse.RawTextHelpFormatter,
        epilog="""
示例用法:
  python filter_protein.py -i all_protein.fa -l 50
  python filter_protein.py -i all_protein.fa -l 100 -o filtered.fa -u short.fa
""")
    parser.add_argument('-i', '--input', required=True, help='输入的蛋白序列文件（fasta格式）')
    parser.add_argument('-l', '--length', type=int, default=50, help='蛋白序列最小长度阈值（默认: 50）')
    parser.add_argument('-o', '--output', default='qualified_protein.fa', help='输出符合长度要求的蛋白文件（默认: qualified_protein.fa）')
    parser.add_argument('-u', '--unqualified', default='unqualified_protein.fa', help='输出未达到长度要求的蛋白文件（默认: unqualified_protein.fa）')
    return parser.parse_args()

def filter_protein_sequences(input_file, min_length, output_file, unqualified_file):
    qualified = []
    unqualified = []

    for record in SeqIO.parse(input_file, "fasta"):
        if len(record.seq) >= min_length:
            qualified.append(record)
        else:
            unqualified.append(record)

    SeqIO.write(qualified, output_file, "fasta")
    SeqIO.write(unqualified, unqualified_file, "fasta")

    print(f"总序列数: {len(qualified) + len(unqualified)}")
    print(f"符合长度要求（≥{min_length}）的序列数: {len(qualified)} -> {output_file}")
    print(f"不符合长度要求的序列数: {len(unqualified)} -> {unqualified_file}")

if __name__ == "__main__":
    args = parse_args()
    filter_protein_sequences(args.input, args.length, args.output, args.unqualified)

注释文件路径

/data2/liyupeng/alice/output/female_anno/mark2_anno/order_rename_gff3_extract

参考

  1. https://www.kdocs.cn/l/crlSVN3OLH24