基因注释软件BRAKER学习笔记

1 BRAKER3的优势

可以使用转录组和蛋白数据
2 基因预测成功的关键
高质量的基因组：short scaffolds太多的话不会得到很准确的结果。
简单的序列名称：如Chr1这种是最好的。
要标记重复序列：the genome should be masked for repeats，避免对重复序列和低复杂度区域预测到基因结构；转录组数据比对是重复序列也会影响；在GeneMark-ES/ET/EP/ETP和AUGUSTUS流程中，softmasking比hardmasking表现更好。
- softmasking：重复区域的核苷酸用小写字母，其他区域的用大写字母。
- hardmasking：重复区域的核苷酸序列用N表示。
注意分支特异性：例如部分真菌会有特殊的分支点等，需要根据物种信息选择合适的参数。
进一步使用预测结果之前需要确定预测的准确性。
3 BRAKER运行模式
3.1 只有基因组
在这种模式下，GeneMark-ES使用基因组进行训练，GeneMark-ES预测的长基因用于训练后续的AUGUSTUS；最终使用AUGUSTUS进行预测。这种模式得到的结果的准确性是最低的。

Pasted image 20250424204809

3.2 基因组+转录组

短读长的RNA-seq数据适合这种预测方法。这种方法要求每个intro被很多的alignments覆盖，但是组装过的转录组比对不适用于这种方法。长读长的RNA-seq数据也是可以的，但是需要注意的是每个进行训练的转录本都需要多次测序并与基因组进行比对。现在BRAKER还不支持短读长和长读长的RNA-seq数据的混合使用。
Pasted image 20250424205041

3.3 基因组+蛋白序列

这种方法尤其适用于没有RNA-seq数据的情况，而且使用的蛋白序列可以是来自和目标物种遗传距离未知的物种，当然，使用近缘物种肯定预测得更准确些。这个方法需要蛋白家族，也就是每个蛋白家族具有代表性的序列应该在数据库中。OrthoDB这个数据库是测试过的。官方提供了OrthoDB v.11和OrthoDB v.12.

Pasted image 20250424205550

3.4 基因组+转录组+蛋白序列

从头预测+转录本+同源蛋白预测。

4 软件安装

4.1 Docker安装

推荐使用docekr进行安装：

1	`docker pull teambraker/braker3`

进入容器：

1	`docker run --rm -it -v $(pwd):/work teambraker/braker3:latest /bin/bash`

4.2 Conda安装

如果没有docker权限，可以尝试使用conda进行安装：

wget https://repo.anaconda.com/archive/Anaconda3-2018.12-Linux-x86_64.sh
bash bin/Anaconda3-2018.12-Linux-x86_64.sh # do not install VS (needs root privileges)
conda install -c anaconda perl
conda install -c anaconda biopython
conda install -c bioconda perl-app-cpanminus
conda install -c bioconda perl-file-spec
conda install -c bioconda perl-hash-merge
conda install -c bioconda perl-list-util
conda install -c bioconda perl-module-load-conditional
conda install -c bioconda perl-posix
conda install -c bioconda perl-file-homedir
conda install -c bioconda perl-parallel-forkmanager
conda install -c bioconda perl-scalar-util-numeric
conda install -c bioconda perl-yaml
conda install -c bioconda perl-class-data-inheritable
conda install -c bioconda perl-exception-class
conda install -c bioconda perl-test-pod
conda install -c bioconda perl-file-which # skip if you are not comparing to reference annotation
conda install -c bioconda perl-mce
conda install -c bioconda perl-threaded
conda install -c bioconda perl-list-util
conda install -c bioconda perl-math-utils
conda install -c bioconda cdbtools
conda install -c eumetsat perl-yaml-xs
conda install -c bioconda perl-data-dumper

BRAKER和AUGUSTUS在conda上也是有的。
个人推荐速度更快的Miniforge替代原生conda。

4.3 软件依赖

AUGUSTUS 3.5.0 F2
GeneMark-ETP (source see Dockerfile)
BAMTOOLS 2.5.1R5
SAMTOOLS 1.7-4-g93586edR6
Spaln 2.3.3d R8, R9, R10
NCBI BLAST+ 2.2.31+ R12, R13
DIAMOND 0.9.24
cdbfasta 0.99
cdbyank 0.981
GUSHR 1.0.0
SRA Toolkit 3.00 R14
HISAT2 2.2.1 R15
BEDTOOLS 2.30 R16
StringTie2 2.2.1 R17
GFFRead 0.12.7 R18
compleasm 0.2.5 R27
4.4 脚本依赖
主要的脚本是braker.pl，在运行braker.pl时会调用以下这些脚本：
`filterGenemark.pl
filterIntronsFindStrand.pl
helpMod_braker.pm
findGenesInIntrons.pl
downsample_traingenes.pl
ensure_n_training_genes.py
get_gc_content.py
get_etp_hints.py
需要确保这些脚本都是可执行状态，也就是有x这个属性；如果不是可执行状态，需要使用chmod a+x *.pl *.py添加可执行属性。
5 运行BRAKER
建议使用softmasked Repeats后的基因组。
5.1 模式选择
5.1.1 转录组
这个流程适用于有转录组但是没有蛋白质数据的情况。BRAKER有多种方式接收传入的转录组数据：
使用SRA数据：有多个SRA数据时使用,进行分隔，BRAKER会自动下载对应的SRA数据：
1
braker.pl --species=yourSpecies --genome=genome.fasta --rnaseq_sets_ids=SRA_ID1,SRA_ID2
使用本地的FASTQ文件：可以是为mapping过的原始数据，同样需要提供FASTQ文件的文件名称和对应的文件的路径，会根据文件的ID自动识别文件，如果是双末端的会识别为SRA_ID1_1.fastq和SRA_ID1_2.fastq这样，如果是单末端的会识别为ID.fastq.
1
2
3
braker.pl --species=yourSpecies --genome=genome.fasta \ --rnaseq_sets_ids=SRA_ID1,SRA_ID2 \ --rnaseq_sets_dirs=/path/to/local/fastq/files/

使用比对后的bam文件

# 类似FASTQ输入方式
braker.pl --species=yourSpecies --genome=genome.fasta \
   --rnaseq_sets_ids=BAM_ID1,BAM_ID2 \
   --rnaseq_sets_dirs=/path/to/local/bam/files/

1 2	`# 直接指定bam文件 braker.pl --species=yourSpecies --genome=genome.fasta --bam=file1.bam,file2.bam`

使用bam文件时，BRAKER默认bam文件是HISAT2比对得到的结果；如果想使用SATR的话，需要加上参数--outSAMstrandField intronMotif.

从比对信息中提取的gff格式信息：

1	`braker.pl --species=yourSpecies --genome=genome.fasta --hints=hints1.gff,hints2.gff`

格式类似于这样：

1 2	`chrName b2h intron 6591 8003 1 + . pri=4;src=E chrName b2h intron 6136 9084 11 + . mult=11;pri=4;src=E`

第二列的b2h和最后一列的src=E是必须要的。

可以将上面这些方法进行组合，但是需要注意的是，如果使用转录组+蛋白组这种方式的话，不建议使用--hints这种方法。

5.1.2 蛋白组

这种方法尤其适用于没有转录组的情况。

执行下面的代码即可运行：

1	`braker.pl --genome=genome.fa --prot_seq=proteins.fa`

推荐使用OrthoDB作为proteins.fa的基础。
也可以使用ProtHint的数据：

1	`braker.pl --genome=genome.fa --hints=prothint_augustus.gff`

对应的文件是这样的：

2R ProtHint intron 11506230 11506648 4 + . src=M;mult=4;pri=4
2R ProtHint intron 9563406  9563473  1 + . grp=69004_0:001de1_702_g;src=C;pri=4;
2R ProtHint intron 8446312  8446371  1 + . grp=43151_0:001cae_473_g;src=C;pri=4;
2R ProtHint intron 8011796  8011865  2 - . src=P;mult=1;pri=4;al_score=0.12;
2R ProtHint start  234524   234526   1 + . src=P;mult=1;pri=4;al_score=0.08;

5.1.2.1 训练和预测UTRs

如果转录组只有bam格式文件，而且基因组是softmasked格式的，那么BRAKER就可以自动训练UTR的参数用于AUGUSTUS预测；预测完成之后就会使用覆盖度等信息进行基因预测，对应的命令是这样的：

1	`braker.pl --species=yourSpecies --genome=genome.fasta --bam=file.bam --UTR=on`

Warnings：

这个流程属于测试性的；
--UTR=on不兼容bamToWig.py这个脚本；
--UTR=on会增加内存需求；
UTR预测有时能优化编码序列的预测准确度，但不是总是能够优化。
5.1.2.2 链特异性的转录组比对结果
对于UTR的训练和预测，链特异性的转录组数据能够提供更多的信息，将bam文件分割成+和-的两个文件，然后运行BRAKER：
1
2
3
braker.pl --species=yourSpecies --genome=genome.fasta \ --bam=plus.bam,minus.bam --stranded=+,- \ --UTR=on
非链特异性的bam文件也可以用，不用于UTR预测，用于基因预测：
1
2
3
braker.pl --species=yourSpecies --genome=genome.fasta \ --bam=plus.bam,minus.bam,unstranded.bam \ --stranded=+,-,. --UTR=on
5.1.3 转录组+蛋白组
首先使用转录组和蛋白组数据训练GeneMark-ETP，输出结果中属于high-confindent的基因用于训练AUGUSTUS，随后使用TSEBRA将这两个软件的结果进行整合。
5.1.3.1 转录组比对
GeneMark-ETP使用Stringtie2对转录本进行组装，也就需要bam文件有XStag，如果比对的时候用的软件是HISAT2，那么就需要加上-dta这个参数；使用STAR的话就要加上--outSAMstrandField intronMotif这个参数选项。有几种代码运行方式：
1
2
3
braker.pl --genome=genome.fa --prot_seq=orthodb.fa \ --rnaseq_sets_ids=SRA_ID1,SRA_ID2 \ --rnaseq_sets_dirs=/path/to/local/RNA-Seq/files/

1 2	`braker.pl --genome=genome.fa --prot_seq=orthodb.fa \ --rnaseq_sets_ids=SRA_ID1,SRA_ID2,SRA_ID3`

1 2	`braker.pl --genome=genome.fa --prot_seq=orthodb.fa \ --bam=/path/to/SRA_ID1.bam,/path/to/SRA_ID2.bam`

5.2 长/短读长转录组

5.3 长读长转录组+蛋白组

如果使用BRAKER的fast-hack模式，那么就需要使用长读长的转录组构造bam文件：

1
2
3

T=48 # 线程数
minimap2 -t${T} -ax splice:hq -uf genome.fa isoseq.fa > isoseq.sam     
samtools view -bS --threads ${T} isoseq.sam -o isoseq.bam

然后拉取适合长读长转录组的BRKAER版本：

1	`singularity build braker3_lr.sif docker://teambraker/braker3:isoseq`

运行：

singularity exec -B ${PWD}:${PWD} braker3_lr.sif braker.pl \
--genome=genome.fa \
--prot_seq=protein_db.fa \
–-bam=isoseq.bam --threads=${T}

Warnings：

不能混合使用短读长和长读长的转录组数据。
isoseq数据的深度严重影响着基因预测的准确度。
5.4 参数选择
5.4.1 —ab_initio
Compute AUGUSTUS ab initio predictions in addition to AUGUSTUS predictions with hints (additional output files: augustus.ab_initio.*. This may be useful for estimating the quality of training gene parameters when inspecting predictions in a Browser.
5.4.2 —augustus_args=”—some_arg=bla”
知道物种适用的AUGUSTUS的话可以把这些参数传入，使用空格进行分隔，如--first_arg=sth --second_arg=sth.
5.4.3 —threads=INT
指定线程数，某些步骤使用单线程，部分步骤使用多线程。设置超过8线程的话，对于多线程步骤会加速计算。
5.4.4 —fungus
是GeneMark-ETP的参数选项，如果是真菌的话使用这个参数，其他不用。
5.4.5 —useexisting
是否使用现成的配置文件，如果使用的话就会覆盖掉软件默认的参数文件。
5.4.6 —crf
是否为AUGUSTUS执行CRF training，生成的参数只是用于最后的预测，使用的话会增加运行时间。
5.4.7 —lambda=int
默认的是参数值是2，如果有大量单个外显子的基因的话，参数可以设置为0.
5.4.8 —UTR=on
见上述步骤。
5.4.9 —addUTR=on
见上述步骤。
5.4.10 —stranded=+,-,.,…
见上述步骤。
6 输出结果
braker.gtf：最终的基因集，不同的模式有不同的结果
- ETPmode：AUGUSTUS和GeneMark-ETP的整合结果
- otherwise：官方不建议使用的模式和结果。
braker.codingseq：最终基因集的编码序列。
braker.aa：最终基因集的蛋白序列。
braker.gff3：最终基因集的注释文件，需要加上--gff3这个参数选项。
Augustus/*：Augustus的输出结果。
GeneMark-E*/genemark.gtf：GeneMark-ES/ET/EP/EP+/ETP的输出结果。
hintsfile.gff：从转录组或者是蛋白质组中输出的结果。
braker_original/*：BRAKER输出的结果，这个结果是运行compleasm之前的结果。
bbc/*：best_by_compleasm.py的输出结果。

如果设置--skipGetAnnoFromFasta参数的话不会生成编码序列和蛋白序列结果。

7 测试练习

7.1 数据下载

基因组数据来自文章Pan-genome analysis of 33 genetically diverse rice accessions reveals hidden genomic variations，选择的是品种是R527，这个基因组不是soft masking的，下载地址https://ricerc.sicau.edu.cn/ 蛋白质组数据来自Rice Resource Center的R498这个水稻品种。转录组数据来自文章的提供的数据。

7.2 运行BRAKER

7.2.1 下载Docker镜像

1	`docker pull teambraker/braker3`

7.2.2 进入容器

1
2
3

docker run --rm -it -v $(pwd):/home/jovyan/braker3 teambraker/braker3:latest /bin/bash

cd braker3 #进入挂载的目录

7.2.3 运行BRAKER

nohup braker.pl --genome=data/R527.fa --prot_seq=data/R498.pep.fa \
        --rnaseq_sets_ids=ERR13336195,ERR13336317,ERR13336454,ERR13336603 \
        --rnaseq_sets_dirs=data/ \
        --threads=65 \
        --crf >> braker3.log 2>&1 &

7.2.4 后台运行容器

上面先进入容器再运行的方法有个问题就想终端关闭或者是exit退出容器之后，进程会被杀死。改用下面的方式运行：

nohup docker run --rm -v $(pwd):/home/jovyan/ teambraker/braker3:latest \
	braker.pl --genome=data/R527.fa --prot_seq=data/R498.pep.fa \
        --rnaseq_sets_ids=CRR279535,CRR279536,CRR279537,CRR279538,CRR279539,CRR279540,CRR279541,CRR279542,CRR279543,CRR279544,CRR279545,CRR279546,CRR279547,CRR279548,CRR279549,CRR279550,CRR279551,CRR279552,CRR279553,CRR279554,CRR279555,CRR279556,CRR279557,CRR279558,CRR279559,CRR279560,CRR279561,CRR279562,CRR279563 \
        --rnaseq_sets_dirs=data/ \
        --threads=65 \
        --crf >> braker.log 2>&1 &

7.2.5 报错处理

遇到报错：

Error: Encountered one or more spaces while parsing the quality string for read ERR13336195.10655764 D00738:173:H2TGTBCX3:2:1201:21091:100335#0/2.  If this 
is a FASTQ file with integer (non-ASCII-encoded) qualities, try re-running with the --integer-quals option.                                                 
terminate called after throwing an instance of 'int'                                               
Aborted (core dumped)                                                                              
(ERR): hisat2-align exited with value 134

先对转录组数据进行处理：

# 创建输出目录                                                                                       
mkdir -p data/clean_data                                                                                    

# 遍历所有 R1 文件（自动匹配对应的 R2 文件）                                                            
for r1_file in data/*_1.fastq.gz; do                                                               
    # 生成 R2 文件名（替换 "R1" 为 "R2"）                                                             
    r2_file="${r1_file/_1/_2}"                                                                                             
                                     
    # 生成输出文件名                                                                                 
    out_r1="data/clean_data/$(basename "$r1_file")"                                                
    out_r2="data/clean_data/$(basename "$r2_file")"   
                                                                                                          
    # 运行 fastp，处理双端数据                                                                        
    fastp -i "$r1_file" -I "$r2_file" \                                                            
          -o "$out_r1" -O "$out_r2" \                                                              
          --html "data/clean_data/$(basename "$r1_file").html" \                                   
          --json "data/clean_data/$(basename "$r1_file").json"                                    
done

最终确定是文件下载的问题，最后只使用ERR13336195这个转录组数据集。

8 比较注释结果

使用下面的脚本评估新转录组数据和注释结果的匹配率（未经过验证）。

#!/bin/bash

# 参数设置
GENOME=$1 # 基因组fasta文件
GFF=$2 # 注释GFF文件
TRANSCRIPTOME=$3 # 转录组fastq文件
THREADS=8 # 并行线程数

# Step 1: 转录组比对到基因组
hisat2-build $GENOME genome_index # 构建索引
hisat2 -p $THREADS -x genome_index -U $TRANSCRIPTOME | samtools sort -O BAM -o aligned.bam
samtools index aligned.bam # 生成索引


# Step 2: 计算覆盖区域并合并连续区间
samtools depth -aa aligned.bam | awk '$3 > 0 {print $1"\t"$2-1"\t"$2}' > coverage.bed
bedtools merge -i coverage.bad > merged_coverage.bed

# Step 3: 提取GFF中的基因坐标
awk '$3=="gene" {print $1"\t"$4-1"\t"$5}' $GFF > genes.bed # 提取基因区域
bedtools sort -i genes.bed > sorted_genes.bed

# Step 4: 比较覆盖区域与注释基因
bedtools intersect -a merged_coverage.bed -b sorted_genes.bed -wa -wb > overlap_results.txt

# Step 5: 统计准确性指标
TOTAL_COV=$(wc -l < merged_coverage.bed)
TOTAL_GENE=$(wc -l < sorted_genes.bed)
MATCHED=$(cut -f1-3 overlap_results.txt | sort -u | wc -l)

# 计算匹配率
MATCH_RATE=$(echo "scale=2; $MATCHED*100/$TOTAL_GENE" | bc)
echo "基因注释匹配率: ${MATCH_RATE}%"

1	`./validate_annotation.sh genome.fa annotations.gff rnaseq.fastq`

9 长读长转录组

9.1 实验和测序设计

9.1.1 测序平台和文库构建

关注准确性和read长度就用PacBio
关注RNA修饰就用ONT
9.1.2 测序深度
通常建议大于2千万条reads或者是深度超过40x，或者根据公式进行计算： $Depth_{sample}=\frac{\sum Depth_{isoforms}}{N_{expressed-isoforms-in-sample}}$
9.1.3 样品重复
尽可能多一点。
9.2 数据处理
先关注用户层面不够用到的。
9.2.1 Mapping
作者推荐的是Minimap2，单是面对小的外显子或者复杂的可变剪切时还是有局限性。如果面对复杂的可变剪切或者是小的外显子，可以使用uLTRA、Graphmap2、deSALT、2passtools、Splam和Magic-BLAST.

9.2.2 转录本鉴定和定量

长读长转录组的主要目的是鉴定和定量转录本。从头预测转录可以发现新的转录本，但是耗时耗计算资源；使用参考基因组对转录本定量的话就可能会丢失新的转录本信息。根据不同的目的使用不同的方法和策略。

对转录本进行鉴定的方法可以分为基于聚类、基于graph和基于class的方法。

9.3 下游分析

9.3.1 差异表达和可变剪切

9.3.2 基因组注释

可以使用的软件有Augustus、SNAP、MAKER和BRAKER.

生物信息学

#生物信息学

基因注释软件BRAKER学习笔记

https://lixiang117423.github.io/article/braker/

作者

李详【Xiang LI】

发布于

2025年5月9日

许可协议

泛基因组时代中的序列到图映射算法综述上一篇

结构变异可视化工具SVbyEye 下一篇

基因注释软件BRAKER学习笔记

1 BRAKER3的优势

2 基因预测成功的关键

3 BRAKER运行模式

3.1 只有基因组

3.2 基因组+转录组

3.3 基因组+蛋白序列

3.4 基因组+转录组+蛋白序列

4 软件安装

4.1 Docker安装

4.2 Conda安装

4.3 软件依赖

4.4 脚本依赖

5 运行BRAKER

5.1 模式选择

5.1.1 转录组

5.1.2 蛋白组

5.1.2.1 训练和预测UTRs

5.1.2.2 链特异性的转录组比对结果

5.1.3 转录组+蛋白组

5.1.3.1 转录组比对

5.2 长/短读长转录组

5.3 长读长转录组+蛋白组

5.4 参数选择

5.4.1 —ab_initio

5.4.2 —augustus_args=”—some_arg=bla”

5.4.3 —threads=INT

5.4.4 —fungus

5.4.5 —useexisting

5.4.6 —crf

5.4.7 —lambda=int

5.4.8 —UTR=on

5.4.9 —addUTR=on

5.4.10 —stranded=+,-,.,…

6 输出结果

7 测试练习

7.1 数据下载

7.2 运行BRAKER

7.2.1 下载Docker镜像

7.2.2 进入容器

7.2.3 运行BRAKER

7.2.4 后台运行容器

7.2.5 报错处理

8 比较注释结果

9 长读长转录组

9.1 实验和测序设计

9.1.1 测序平台和文库构建

9.1.2 测序深度

9.1.3 样品重复

9.2 数据处理

9.2.1 Mapping

9.2.2 转录本鉴定和定量

9.3 下游分析

9.3.1 差异表达和可变剪切

9.3.2 基因组注释