18 Jan 2025 6483字 22分
CC BY 4.0 （除特别声明或转载文章外）
如果这篇博客帮助到你，可以请我喝一杯咖啡~

在对基因组重复序列和ncRNA进行注释后，接下来是基因预测和功能注释，这也是寻找功能基因的基础和前提。这里主要记录下怎么用Braker3进行基因组的基因预测（也就是结构注释）。

基因预测的方法主要有三种：

基于隐马尔可夫模型的自训练和迭代，获得从头预测的基因结构模型（软件Augustus、GeneMark-ES等）
基于已发表的近缘物种基因序列、蛋白序列的同源预测（软件DIAMOND、GeMoMa等）
基于本物种的RNA-Seq转录组数据，比对基因组内含子结构模型和基因侧翼序列信息（软件Hisat2、STAR等）

一般的流程是将以上三种基因预测结果通过软件EvidenceModeler（EVM）进行整合，最终得到预测结果gff文件。网上对于上面的流程有很多教程，针对不同软件有不同的设置，比如知乎的这篇文章使用AUGSTUS + Geneid + GeneMark + GeMoMa + GenomeThreader + Exonerate 进行基因结构预测 - 知乎 (zhihu.com)

为了简化流程，现在也有越来越多的基因预测pipeline工具得以开发，比较有名的就是Braker和Maker，感兴趣的话可以做两者预测结果的比较，我这里就用发表时间比较近的Braker3为例。

Braker本质上是一个结合了多种基因组注释工具的perl程序，其核心为braker.pl文件。

1. 使用braker3对基因进行预测

1.1. 安装braker3

braker3的安装请查看：

1.2. braker3的使用

基因预测的方法主要有三种：

基于隐马尔可夫模型的自训练和迭代，获得从头预测的基因结构模型（软件Augustus、GeneMark-ES等）
基于已发表的近缘物种基因序列、蛋白序列的同源预测（软件DIAMOND、GeMoMa等）
基于本物种的RNA-Seq转录组数据，比对基因组内含子结构模型和基因侧翼序列信息（软件Hisat2、STAR等）

Braker本质上是一个结合了多种基因组注释工具的perl程序，其核心为braker.pl文件。

上图是整合RNA-Seq数据和蛋白数据跑Braker的流程图，需要注意基因组文件genome.fa在输入前需要需要进行softmasking（重复序列屏蔽为小写字母），官方建议不要用hardmasking（重复序列屏蔽为N），hardmasking后预测的基因数量会偏少，因为重复序列中可能也有功能基因的部分信息，屏蔽为N后就无法检测到了。

对RNA-Seq数据的处理，首先是通过SRA tookit将SRA ID对应的fastq数据下载下来（如果本来就是fastq格式就不需要这一步），用Hisat2比对到softmasking后的参考基因组并生成bam文件，再用stringtie进行转录本组装。

GeneMark-ETP以组装后的转录本和同源蛋白数据库作为输入数据进行训练和预测，之后再用AUGUSTUS软件结合上一步的预测结果进行训练和预测，最后用TSEBRA对预测的基因集进行整合，得到最终的gtf结果文件。

下面的完整步骤：

步骤 1: 文件结构准备

00.ProteinData文件，包含OrthoDB植物蛋白质数据库：Viridiplantae.fa

01.RawData文件，包含以下RNA-seq数据：

步骤2:处理RNAseq文件

#!/bin/bash

# 设置变量
GENOME_FILE="/data/wangxingbin/00.Project/01.Toona_Annotation/hap1/01.repeatMasker/02.TE/02.repeatMasker/repeatmasker_output_softmasking/hap1_T2T.fasta.masked"  # 参考基因组文件路径（如果需要未掩盖的基因组，请使用相应的文件）
GENOME_FASTA=“genome.fasta.masked”
INDEX_DIR="./hista2_idx"  # Hisat2索引文件夹
RAW_DATA_DIR="/data/wangxingbin/00.Project/01.Toona_Annotation/05_RNA_T7/01.RawData"  # RNA-seq原始数据文件夹
THREADS=36  # 使用的线程数
OUTPUT_DIR="./hista2_out"  # 输出结果的文件夹

# 如果输出文件夹不存在，则创建

mkdir -p $OUTPUT_DIR

# 如果Hisat2索引还未建立，执行索引构建
if [ ! -d "$INDEX_DIR" ]; then
echo "正在构建Hisat2索引..."
mkdir -p $INDEX_DIR
ln -s $GENOME_FILE ./$INDEX_DIR/$GENOME_FASTA
cd $INDEX_DIR
hisat2-build -p $THREADS "./$GENOME_FASTA" idx
cd ..
else
echo "Hisat2索引已经存在，跳过索引构建"
fi

# 处理每个样本的RNA-seq数据
for SAMPLE in leaf LY NY root steam; do
# 获取正向和反向读段文件路径
R1="${RAW_DATA_DIR}/${SAMPLE}/${SAMPLE}_1.fq.gz"
R2="${RAW_DATA_DIR}/${SAMPLE}/${SAMPLE}_2.fq.gz"

# 输出BAM文件的路径
BAM_OUTPUT="${OUTPUT_DIR}/${SAMPLE}.bam"

# 执行比对并生成BAM文件
echo "正在处理样本：$SAMPLE"
hisat2 --dta -x $INDEX_DIR/idx -p $THREADS -1 $R1 -2 $R2 | samtools sort -@ 10 -o $BAM_OUTPUT &
done

# 等待所有背景任务完成
wait

echo "所有样本的比对处理已完成！"

步骤3：运行braker3

#!/bin/bash
#export GENEMARK_PATH=$GENEMARK_PATH/gmes
# 设置使用的核数和最大运行时间
threads=36                      # 设置使用的线程数

# 工作目录
wd=braker3_output

# 检查并删除已经存在的工作目录
if [ -d "$wd" ]; then
    rm -r "$wd"                      # 删除已存在的目录，避免冲突
fi


# 设置 BRaKER3 的安装路径 (需要根据实际路径进行修改)
BRAKER_DIR="/data/wangxingbin/Software/braker-3.0.8"

# 设置基因组文件和蛋白质数据库路径 (请替换成你实际的文件路径)
GENOME_PATH="/data/wangxingbin/00.Project/01.Toona_Annotation/hap1/01.repeatMasker/02.TE/02.repeatMasker/repeatmasker_output_softmasking/hap1_T2T.fasta.masked"
PROTEIN_DB="/data/wangxingbin/00.Project/01.Toona_Annotation/hap1/03.functionalGene/01.braker3/00.softmasking_braker3/Viridiplantae.fa"

# 设置 RNA-Seq 数据路径 (按需替换)
RNA_SEQ_DIR="./hista2_out"
#RNA_SEQ_FILES="leaf,LY,NY,root,steam"  # RNA-Seq 文件 IDs，按样本分组

# 直接执行 BRaKER3，不使用容器
echo "开始执行，BRAKER3!"
# 创建工作目录
mkdir -p "$wd"

nohup perl "$BRAKER_DIR/scripts/braker.pl" --genome="$GENOME_PATH" --prot_seq="$PROTEIN_DB" --threads "$threads" --softmasking --workingdir="$wd" --bam $RNA_SEQ_DIR/leaf.bam,$RNA_SEQ_DIR/LY.bam,$RNA_SEQ_DIR/NY.bam,$RNA_SEQ_DIR/root.bam,$RNA_SEQ_DIR/steam.bam --skipOptimize --gff3 > "./run_braker3.log" 2>&1 &
# 提示信息
echo "BRAKER3 gene prediction started in the background."
#nohup bash run_braker3.sh > braker3_output.log" 2>&1 &
# 说明：
# - --genome: 输入的软屏蔽过的基因组文件（即 final_genome_softmasked.fa）。
# - --prot_seq: 输入的蛋白质数据库，这里需要提供适合的参考蛋白质数据库路径。
# - --rnaseq_sets_dirs: RNA-Seq 文件所在的目录路径。包含子目录（例如：leaf, LY, NY, root, steam）。
# - --rnaseq_sets_ids: 用逗号分隔的 RNA-Seq 样本 IDs，确保这些样本数据已经按照指定的格式组织。
# - --workingdir: 用于存放 `BRAKER3` 结果的目录。
# - --threads: 设置计算时使用的线程数。
# - --softmasking: 启用软屏蔽(小写字母屏蔽输入方式)。
# - nohup: 在后台运行作业，即使你退出终端也不会中断任务。

1.3. 结果文件说明

成功运行后结果文件如下所示：

braker.codingseq 基因集的fasta文件
braker.aa 基因集的蛋白序列的fasta文件
braker.gtf Braker预测的基因集最终文件
braker.gff3 Braker预测的基因集的gff3格式文件
Augustus文件夹 AUGUSTUS预测的基因集
GeneMark-ETP文件夹 GeneMark-ETP预测的基因集以及其他中间文件
hintsfile.gff 从RNA-Seq数据和蛋白库数据中提取的外部证据，有内含子等信息

这些结果中是不包含UTR的，需要调用GUSHR且在参数中添加–addUTR=on预测UTR信息

2. 蛋白质数据库下载

同源蛋白库建议用官方推荐的OrthoDB或者找几个模式植物的蛋白数据合并。

2.1. OrthoDB蛋白数据库下载

官方推荐使用OrthoDB数据库作为同源蛋白来源，地址：https://bioinf.uni-greifswald.de/bioinf/partitioned_odb11/。

The following partitions are available for download:

Metazoa.fa.gz (4.6 GB, md5sum: a0e6ff3fe2061cc7a9027616691e6060)
Vertebrata.fa.gz (3.1 GB, md5sum: 27bf52491da63469edc6b323be51bda1)
Viridiplantae.fa.gz (1.3 GB, md5sum: 59bd04dd98102b8cd51b1d0d8d993105)
Arthropoda.fa.gz (1.2 GB, md5sum: 7aa03dbcef38517752be16004ad9cb76)
Eukaryota.fa.gz (8.5 GB, md5sum: ff74dbac56700926b9acc39e17ce323a)
Fungi.fa.gz (2.2 GB, md5sum: 3cd9cb99d70730093f3f64ebe2f15f9b)
Alveolata.fa.gz (201 MB, md5sum: 1c38818489a598843e985e73b45c0ade)
Stramenopiles.fa.gz (108 MB, md5sum: eab4b468d1fcc91c8477d04c797e37aa, should be combined with proteins from another clade for BRAKER, or use Eukaryota.fa.gz)
Amoebozoa.fa.gz (38 MB, md5sum: fa094c9d01fb6bafd5dfa701b1f4efaf, should be combined with proteins from another clade for BRAKER, or use Eukaryota.fa.gz)
Euglenozoa.fa.gz (99 MB, md5sum: 5cb75b7d80b953e1d326b4f863ed6739, should be combined with proteins from another clade for BRAKER, or use Eukaryota.fa.gz)

我这里要分析的物种是香椿(植物)，所以下载绿色植物蛋白库Viridiplantae.fa.gz：

# 下载
wegt -c "https://bioinf.uni-greifswald.de/bioinf/partitioned_odb11/Viridiplantae.fa.gz"
gunzip Viridiplantae.fa.gz #解压

2.2. 自建数据库

在NCBI网站上直接找一些组装注释结果较好的模式生物和近缘物种蛋白：

# 拟南芥（TAIR10.1）
wget -c https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/735/GCF_000001735.4_TAIR10.1/GCF_000001735.4_TAIR10.1_protein.faa.gz
# 栽培烟草（Ntab-TN90）
wget -c https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/715/135/GCF_000715135.1_Ntab-TN90/GCF_000715135.1_Ntab-TN90_protein.faa.gz
# 水稻（IRGSP-1.0）
wget -c https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/001/433/935/GCF_001433935.1_IRGSP-1.0/GCF_001433935.1_IRGSP-1.0_protein.faa.gz
# 近缘物种coffea arabica
wget -c 
# 近缘物种coffea canephora
wget -c

gunzip *.gz
cat *.faa > mydb_proteins.fasta

近缘物种的同源蛋白是如何找到的：

plant Biology - Usadel lab (plabipd.de)这个网站记录了多种已发表的植物基因组文章和数据，点击cladogram view可以直观地看到已测过基因组的植物学名和树状图，比如我要找的物种是夹竹桃科（Apocynaceae），直接ctrl + F 就可以定位到夹竹桃科所处的进化节点。

用Toona近缘且已发表基因组的植物学名，一个一个去搜NCBI网站的Genome库，有protein序列的就可以直接下载。

基因组注释从0开始（4）：使用braker3进行基因预测