基因组注释（4）——基因预测

目前为止（2023年4月3日）Braker的最新版本为3.0.2，conda上能搜到的最新版本只有2.1.6，因此不建议用conda安装，尤其是最新的版本Braker可以直接使用RNA-seq和蛋白数据，整合GeneMark-ETP和AUGUSTUS训练和预测基因，对于预测结果有较高的支持度。

因为整个pipeline包含了十几个注释用的软件，用到的perl模块也非常多（数了一下配置环境需要安装20个perl模块），还是推荐用给官方给的**container**。

1.1 申请和下载GeneMark-ETP密钥

在Braker3中使用RNA-seq数据和蛋白数据预测基因，都要用到GeneMark-ETP这个软件。但是这个软件不能直接用，需要到GeneMark网站申请和下载对应的密钥文件放在集群用户的家目录中。

申请完成之后获得名称为gm_key_64.gz的密钥文件，解压之后命名为.gm_key（注意点号）并上传到集群用户的家目录下即可。

1
2

gunzip gm_key_64.gz
mv gm_key_64 .gm_key

1.2 创建Braker3镜像

这一步在dockerhub网站的Braker3仓库中有详细说明，我这里选择创建singularity镜像：

1

singularity build braker3.sif docker://teambraker/braker3:latest

得到的braker3.sif就是Braker3的singularity image

创建braker3镜像文件的环境变量：

1

export BRAKER_SIF=/your/path/to/braker3.sif

可以复制三个示例脚本到当前目录：

1
2
3

singularity exec -B $PWD:$PWD braker3.sif cp /opt/BRAKER/example/singularity-tests/test1.sh .
singularity exec -B $PWD:$PWD braker3.sif cp /opt/BRAKER/example/singularity-tests/test2.sh .
singularity exec -B $PWD:$PWD braker3.sif cp /opt/BRAKER/example/singularity-tests/test3.sh .

在本地申请计算资源并跑一下三个示例脚本：

1
2
3
4
5

salloc -n 50	# 申请50个核跑test.sh，注意不要在登录节点直接运行计算程序
bash test1.sh # tests BRAKER1
bash test2.sh # tests BRAKER2
bash test3.sh # tests BRAKER3
exit	# 退出并释放计算资源

官方提供了4种BRAKER pipeline 模式：

• RNA-Seq数据跑BRAKER
• 蛋白数据跑BRAKER
• 整合RNA-Seq数据以及蛋白数据跑BRAKER
• 整合短读长与长读长的RNA-Seq数据以及蛋白数据跑BRAKER

4种pipeline模式在调用软件的方法上有区别，根据自己手上有的数据选择用哪种，我这里选择第三种。

上图是整合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结果文件。

barker.sh脚本可以如下编写:

1
2
3
4
5
6
7
8
9
10
11
12

#!/bin/bash
#SBATCH -n 48
#SBATCH -t 7200

wd=baima_pre

if [ -d $wd ]; then
    rm -r $wd
fi

singularity exec -B ${PWD}:${PWD} ${BRAKER_SIF} braker.pl --genome=/public/home/wlxie/biosoft/db_data/baima/RepeatMasker_soft/genome.nextpolish.fasta.masked --prot_seq=/public/home/wlxie/busco_soft/busco/test_data/eukaryota/busco_downloads/lineages/eudicots_odb10/refseq_db.faa --softmasking --threads 48 --workingdir=${wd} --rnaseq_sets_dirs=/public/home/wlxie/RNAseq/BYT2022020901/rnaseq/baima --rnaseq_sets_ids=4-216031965_raw

几个参数的解释：

genome softmasking后的基因组文件位置

prot_seq 同源蛋白库的文件位置

–softmasking mask的方式

–threads 跑程序用的核数

–workingdir 工作目录位置

rnaseq_sets_dirs RNA-Seq数据所在目录

–rnaseq_sets_ids 双端测序数据文件前缀（比如我这里是4-216031965_raw_1.fq和4-216031965_raw_2.fq）

说明一下同源蛋白来源于前面做BUSCO评估的真双子叶植物单拷贝直系同源库，怎么来的详情可见这篇博客（同源蛋白库建议用官方推荐的OrthoDB或者找几个模式植物的蛋白数据合并，见博客最下方的更新）

前面说过塔大集群的计算节点没有安装singularity，所以在运行该容器的时候要在申请核在本地跑程序，并且用screen维持当前会话：

1
2
3
4
5
6
7

screen -S singularity	# 创建singularity会话
salloc -n 48 -t 7200	# 申请计算资源
bash barker.sh

screen -r singularity	# 进入singularity会话
exit	# 退出会话
exit	# 运行结束释放计算资源

正常跑完花费了9个小时时间（200Mbp大小的基因组），如果中途不幸出bug，braker支持有限度的断点重新运行，主要分为以下三个阶段：

只要有中间文件存在，就可以在这三个阶段继续加入其他参数，跳过已经运行的阶段继续运行，详情可以看官方文档https://github.com/Gaius-Augustus/BRAKER#starting-braker-on-the-basis-of-previously-existing-braker-runs

可以在前面给定的工作目录中看到如下结果文件:

• braker.gtf——Braker预测的基因集，包括了各种不同的基因结构预测结果
• braker.codingseq——fasta格式的编码序列基因集（基因序列）
• braker.aa——fasta格式的蛋白序列基因集（蛋白序列）
• braker.gff3——需要–gff3参数指定，这里我没有，就是基因集的gff3格式
• Augustus/*——AUGUSTUS预测的基因集（包括gtf文件、基因序列和蛋白序列）
• GeneMark-ETP/*——GeneMark-ETP预测的基因集以及其他中间文件
• hintsfile.gff——从RNA-Seq数据和蛋白库数据中提取的外部证据数据

可以通过awk命令查看gft文件的第三列，查看预测的编码蛋白基因数量和转录本数量：

1
2
3
4
5

$ awk '$3=="gene"' braker.gtf | wc -l
17869

$ awk '$3=="transcript"' braker.gtf | wc -l
20832

后续就可以对这些预测的基因做质量评估，然后比对各数据库做功能注释。

1. OrthoDB蛋白数据库下载

官方推荐使用OrthoDB数据库作为同源蛋白来源。

1
2
3
4
5
6
7
8
9
10
11
12

# 真菌
Fungi: https://v100.orthodb.org/download/odb10_fungi_fasta.tar.gz
# 后生动物
Metazoa: https://v100.orthodb.org/download/odb10_metazoa_fasta.tar.gz
	# 节肢动物
	- Arthropoda: https://v100.orthodb.org/download/odb10_arthropoda_fasta.tar.gz
	# 脊椎动物
	- Vertebrata: https://v100.orthodb.org/download/odb10_vertebrata_fasta.tar.gz
# 单细胞生物
Protozoa: https://v100.orthodb.org/download/odb10_protozoa_fasta.tar.gz
# 绿色植物
Viridiplantae: https://v100.orthodb.org/download/odb10_plants_fasta.tar.gz

我这里要分析的物种是植物，所以下载最后一个绿色植物蛋白库：

1
2
3

nohup wget https://v100.orthodb.org/download/odb10_plants_fasta.tar.gz &
tar zxvf odb10_plants_fasta.tar.gz
cat plants/Rawdata/* > plant_proteins.fasta

绿色植物蛋白库约780Mb大小，建议用wget下载，如果windows下载再ftp拖拽上传到集群可能会损坏文件（并且没有md5效验码没办法确认是否真的损坏）。

解压并将所有数据合并到一个文件plant_proteins.fasta中，截至目前2023年4月7日为止，这个数据库共有3510742条蛋白序列，总文件大小为1.4Gb，比原来我比对的蛋白库大了1500倍。而蛋白比对是一个很缓慢的过程，因此这一步预测基因的时间将会很长，可以根据自己要做的物种确定用哪些注释比较完善的模式生物的蛋白库。

2. 自建蛋白数据库

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

1
2
3
4
5
6
7
8
9
10
11
12
13

# 拟南芥（TAIR10.1）
wget 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 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 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 https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/003/713/225/GCF_003713225.1_Cara_1.0/GCF_003713225.1_Cara_1.0_protein.faa.gz
# 近缘物种coffea canephora
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/900/059/795/GCA_900059795.1_AUK_PRJEB4211_v1/GCA_900059795.1_AUK_PRJEB4211_v1_protein.faa.gz

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

顺便记录一下近缘物种的同源蛋白是如何找到的：

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

然后用Apocynaceae祖先节点和子节点的已发表基因组的植物学名，一个一个去搜NCBI网站的Genome库，有protein序列的就可以直接下载。

两种同源蛋白建库方式预测的基因数量和花费的时间：

用自建蛋白数据库跑braker预测的基因数更多，且花费时间更短。后续以该预测结果继续分析。

BRAKER v3.0.2 中存在bug，最后生成的braker.gff3文件中缺少GeneMark.hmm3预测结果。该bug在BRAKER v3.0.3版本中修复，请使用最新版本。