全自动基因结构注释脚本gene_anno.sh

By liyupeng, 14 April, 2026

Forums

示例命令

s，指定软屏蔽后的基因组序列
g，指定硬屏蔽后的基因组序列
p，指定用于注释的蛋白序列文件
d，指定转录组测序数据文件位置
1. 所有转录组测序数据应该以（_1/_2）.fa.gz结尾
t，指定使用的线程数量
所有输入文件必须使用绝对路径
重复使用该脚本需要删除文件：/opt/Augustus/config//species/torreya02，这里的在braker预测过程中生成的已有augustus物种参数模型需要删除掉

bash test.sh \
  -s /home/server/work/test/Chr2.fna.soft.masked.fasta \
  -g /home/server/work/test/Chr2.fna.hard.masked.fasta \
  -p /home/server/work/test/allproteinv2.fasta \
  -d /home/server/work/test/trans/ \
  -t 30

README

All software for gene structure annotation was installed using Miniconda or cloned from GitHub and installed in /home/soft/.

Usage

Run the script repeat_anno.sh as follows to start the gene annotation pipeline:

bash repeat_anno.sh -i test.fasta -t 40
Options:
-s: Specify the soft-masked genome sequence file (FASTA format).
-g: Specify the hard-masked genome sequence file (FASTA format).
-p: Specify the protein sequence file for annotation.
-d: Specify the directory where transcriptome sequencing data files are stored. All transcriptome data files should end with _1.fa.gz and _2.fa.gz for paired-end reads.
-t: Specify the number of threads to use for parallel processing.
Note:
All input files must be specified using absolute file paths.
All transcriptome sequencing data should end with _1.fa.gz and _2.fa.gz for paired-end reads.
Output

The following outputs are generated by the pipeline:

Trinity: Transcriptome sequencing data aligned with the genome (SAM/BAM files), along with assembled transcript sequences.
PASA: Gene annotation based on transcript alignment with the genome, and its associated database.
Braker: Gene prediction results from Braker3.
GTH: Gene predictions from GenomeThreader.
EVM: Annotations generated by EVidenceModeler by integrating various prediction results (from PASA, Braker, GTH, etc.).
PASA_Update: Updated PASA annotations after two rounds of PASA pipeline, including annotations for alternative splicing and untranslated regions.
Results: Final annotated files generated from PASA_Update. These files are renamed and moved to the results directory.

输出文件

trinity：转录组测序数据和基因组序列的比对文件（sam/bam），转录本
pasa：基于转录本序列和基因组序列比对预测的基因注释内容，数据库
braker：braker3的预测结果
gth：genomethreader的预测结果
evm：evm整合各种来源预测结果，生成的一般注释文件
pasa_update：经过pasa两轮更新后，添加了可变剪切和非翻译区域的注释内容
results：自动化注释流程的最终生成文件，从pasa_update里面复制并改名了的

源代码

#!/bin/bash

# 初始化变量
genome1=""
genome2=""
threads=""
trans_dir=""
protein=""

# 使用 getopts 来处理 -i, -t 和 -d 参数，将参数后面的值赋给变量
while getopts "s:g:t:d:p:" opt; do
    case $opt in
        s)  # 处理 -s 参数
            genome1=$OPTARG  # 将 -s 参数后面的值存储在变量 genome1 中（软屏蔽后的基因组序列文件）
            ;;
        g)  # 处理 -g 参数
            genome2=$OPTARG  # 将 -g 参数后面的值存储在变量 genome2 中（硬屏蔽后的基因组序列文件）
            ;;
        t)  # 处理 -t 参数
            threads=$OPTARG  # 将 -t 参数后面的值存储在变量 threads 中（线程数）
            ;;
        d)  # 处理 -d 参数
            trans_dir=$OPTARG  # 将 -d 参数后面的值存储在变量 trans_dir 中（转录组测序数据存放路径）
            ;;
        p)  # 处理 -p 参数
            protein=$OPTARG  # 将 -p 参数后面的值存储在变量 protein 中（蛋白质序列路径）
            ;;
        \?)  # 处理无效参数
            echo "无效参数：-$OPTARG"
            exit 1
            ;;
    esac
done

#启动虚拟环境trinity
source /root/miniconda3/bin/activate trinity

# 1. 创建trinity工作目录并进入，HISAT2建立基因组索引
mkdir trinity
cp $genome1 target.fasta
mv target.fasta trinity/
cd trinity
hisat2-build --large-index target.fasta target_index -p $threads

# 2. 循环比对转录组测序数据到基因组
# 遍历指定目录中的所有以 _1.fq.gz 结尾的文件
for fq1 in $trans_dir/*_1.fq.gz; do
     #获取文件前缀（不带后缀）
    base_name=$(basename $fq1 _1.fq.gz)
    
     #构建配对文件的路径，假设配对文件名为 _2.fq.gz
    fq2="$trans_dir/${base_name}_2.fq.gz"

     #如果找到了配对文件 fq2
    if [ -f "$fq2" ]; then
        #执行比对，将比对结果输出到 SAM 文件，文件名基于输入文件的前缀命名
        hisat2 -x target_index -1 $fq1 -2 $fq2 -p $threads -S "${base_name}.sam"
    else
        #如果没有找到配对文件，输出错误信息
        echo "未找到配对的文件: ${base_name}_2.fq.gz"
    fi
done

#3 创建一个文件夹，对比对的sam文件进行过滤，合并
mkdir raw_sam
mv *.sam raw_sam/

#  循环处理 raw_sam 文件夹中的每个 SAM 文件
for sam_file in raw_sam/*.sam; do
    # 获取文件名的前缀（不带扩展名）
    base_name=$(basename $sam_file .sam)
    
    # 使用 samtools 对每个 sam 文件进行排序并输出到 BAM 文件
    samtools sort -@ $threads -o "${base_name}_sorted.bam" "$sam_file"
done

#移动bam文件，然后合并这些bam文件，生成merged.bam
mkdir bam_files
mv *.bam bam_files/
samtools merge -@ $threads merged.bam bam_files/*.bam 

#4. Trinity 有参组装
Trinity --genome_guided_bam merged.bam --genome_guided_max_intron 1600000 --max_memory 120G --CPU $threads --output trinity_output

#切换会base环境
source /root/miniconda3/bin/activate base

#5. PASA回帖
cd ../
mkdir pasa
cp trinity/trinity_output/Trinity-GG.fasta pasa/
cd pasa
cp /home/alignAssembly.config ./

dir=$(pwd)
echo "DATABASE=$dir/t800" >> alignAssembly.config

# PASA预测命令行
/usr/local/src/PASApipeline/Launch_PASA_pipeline.pl -c alignAssembly.config -C -R -g $genome1 -t Trinity-GG.fasta --ALIGNERS blat,gmap --CPU $threads

#预测开放阅读框
/usr/local/src/PASApipeline/scripts/pasa_asmbls_to_training_set.dbi --pasa_transcripts_fasta chr2.assemblies.fasta --pasa_transcripts_gff3 chr2.pasa_assemblies.gff3

#回到主目录
cd ../
mkdir braker
#启动braker虚拟环境，运行命令行
source /root/miniconda3/bin/activate braker3
cp trinity/merged.bam braker/
cp $protein homology_protein.fa
mv homology_protein.fa braker/

cd braker
braker.pl --threads=$threads --species=torreya02 --genome=$genome1 --softmasking --bam=merged.bam --gff3 --prot_seq=homology_protein.fa &
cd ../

#genomethreader预测
mkdir gth
cp $protein homology_protein.fa
mv homology_protein.fa gth/
cd gth
gth -genomic $genome2 -protein homology_protein.fa -gff3out -skipalignmentout -o gth_output.gff

#返回主目录，创建evm文件节哀
cd ../
mkdir evm
cp pasa/chr2.pasa_assemblies.gff3 evm/
cp pasa/chr2.assemblies.fasta.transdecoder.genome.gff3 evm/
cp gth/gth_output.gff evm/
cp braker/braker/braker.gtf evm/

#切换到braker3环境
source /root/miniconda3/bin/activate braker3

#gff文件预处理
cd evm/
/home/soft/EVidenceModeler/EvmUtils/misc/braker_GTF_to_EVM_GFF3.pl braker.gtf > denovo.gff3
/home/soft/EVidenceModeler/EvmUtils/misc/genomeThreader_to_evm_gff3.pl gth_output.gff > evm_pro.gff3
cp chr2.pasa_assemblies.gff3 evm_rna.gff3

cat denovo.gff3 chr2.assemblies.fasta.transdecoder.genome.gff3 > evm_denovo.gff3

EVidenceModeler --sample_id final --genome $genome1 --weights /home/weights.txt --gene_predictions evm_denovo.gff3 --protein_alignments evm_pro.gff3 --transcript_alignments chr2.pasa_assemblies.gff3 --segmentSize 100000 --overlapSize 10000 --CPU $threads

#创建pasa_update，移动evm整合结果
cd ../
mkdir pasa_update
cp evm/chr2.EVM.gff3 pasa_update/
cd pasa_update

#切换回base环境
source /root/miniconda3/bin/activate base

#创建update1文件夹，用于存储第一轮更新数据
mkdir update1
cd update1

#pasa更新，数据导入
/usr/local/src/PASApipeline/scripts/Load_Current_Gene_Annotations.dbi -c ../../pasa/alignAssembly.config -g $genome1 -P ../chr2.EVM.gff3

#第一轮更新
/usr/local/src/PASApipeline/Launch_PASA_pipeline.pl -c ../../pasa/alignAssembly.config -A -g $genome1 -t ../../trinity/Trinity-GG.fasta --CPU $threads

#第二轮更新
cd ..
mkdir update2
cd update2


#导入update1里面更新过1轮的gff3
/usr/local/src/PASApipeline/scripts/Load_Current_Gene_Annotations.dbi -c ../../pasa/alignAssembly.config -g $genome1 -P ../update1/*.gff3

#第二轮更新
/usr/local/src/PASApipeline/Launch_PASA_pipeline.pl -c ../../pasa/alignAssembly.config -A -g $genome1 -t ../../trinity/Trinity-GG.fasta --CPU $threads

#文件重命名
mv *.gff3 final.evm.update2.gff3
cp final.evm.update2.gff3 ../../

cd ../../
mkdir results
mv final.evm.update2.gff3 results/