0基础学习基因组三代测序组装（7）——基因组组装质量评估（BUSCO、LAI指数）

BUSCO（Benchmarking Universal Single-Copy Orthologs）评估是在基因含量层面上评估基因组完整性。简单来说，通过已有的直系同源数据库进行基因组比对，同源的生物之间有保守基因序列，能比对上的基因数越多说明组装的结果越靠谱。

安装过程：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

# 1. 源码安装（需要安装前置软件）
git clone https://gitlab.com/ezlab/busco.git
cd busco/
python setup.py install

# 前置软件：
https://biopython.org/
https://pandas.pydata.org/
https://jgi.doe.gov/data-and-tools/software-tools/bbtools/
https://ftp.ncbi.nlm.nih.gov/blast/executables/blast+/LATEST
http://bioinf.uni-greifswald.de/augustus/
https://github.com/soedinglab/metaeuk
https://github.com/hyattpd/Prodigal
http://hmmer.org/
https://github.com/smirarab/sepp/
https://www.r-project.org/

# 2. conda安装（推荐）
conda install -c conda-forge -c bioconda busco=5.3.2

conda安装可能会比较慢，需要多试几次。实在不行就源码下载编译，不过需要下载非常多的前置软件，不同软件可能会有环境冲突问题、gcc版本问题等等（我花了大半天时间在折腾环境）。安装之后通过busco -h查看是否安装成功，如果提示缺什么软件就用conda补上（我当前环境中没有安装pandas就会有提示）。

通过busco --list-datasets可以查看当前有哪些物种的数据库，我的植物是双子叶龙胆目，这里的数据库只有真双子叶植物（eudicots）分支离的最近，因此选择这个数据库，v5版本所有单拷贝直系同源数据库网址https://busco-data.ezlab.org/v5/data/lineages/

下载的数据库放在busco_downloads文件夹中，解压即可使用：

1
2

nohup wget https://busco-data.ezlab.org/v5/data/lineages/eudicots_odb10.2020-09-10.tar.gz &
tar -zxvf eudicots_odb10.2020-09-10.tar.gz

busco的详细参数可以看官网的user guide User guide BUSCO v5.4.4 (ezlab.org)

简单讲一讲格式和能用到的参数：

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

busco -i [SEQUENCE_FILE] -l [LINEAGE] -o [OUTPUT_NAME] -m [MODE] [OTHER OPTIONS]
'''
主要参数：
 -i		序列文件位置
 -l		下载的同源物种保守基因数据库位置
 -o		输出文件名
 -m		模式，分为genome,proteins,transcriptome三种
 其他参数：
 --cpu		设置cpu数量
 --download		在线下载数据库，根据分类有"all"、"prokaryota"、"eukaryota"和"virus" （不推荐，速度慢）
 --offline		离线模式，不会更新数据库
'''

以下是我跑的程序，大约用了1个小时：

1

busco -i /public/home/wlxie/NextPolish/01_rundir/genome.nextpolish.fasta -l /public/home/wlxie/busco_soft/busco/test_data/eukaryota/busco_downloads/lineages/eudicots_odb10 -o baima -m genome --cpu 8 --offline

截至2023/02/28，真双子叶植物库有2326个保守BUSCO基因序列，比对结果文件在short_summary.specific.xxx.xxx.txt中，如下：

• Complete BUSCOs (C) 多少个基因完全比对上BUSCOs
• Complete and single-copy BUSCOs (S) 多少个基因比对上单拷贝的BUSCOs
• Complete and duplicated BUSCOs (D) 多少个基因比对上多拷贝的BUSCOs
• Fragmented BUSCOs (F) 多少个基因部分比对上BUSCOs，可能基因只是部分存在
• Missing BUSCOs (M) 多少个基因没有比对上BUSCOs，可能这些直系同源基因是缺失的

从上面的数据看，组装结果还是不错的。从中也可以看到BUSCO运行的两个步骤：用metaeuk进行基因预测（真核生物可以用tBLASTn与对应的BUSCO数据库序列进行比对从而确定候选区域，然后使用 Augustus 软件进行基因结构预测，两个软件可以替代metaeuk，详细参数见官网），以及HMMER进行同源基因的比对，从而评估基因组组装的完整性。

官方还提供了相应的python程序绘制结果图（调用了R包ggplot2），先将BUSCO结果文件放到新建的文件夹，运行相应的py程序，指定工作目录即可：

1
2
3
4
5

mkdir summaries

cp baima/short_summary.specific.eudicots_odb10.baima.txt summaries

generate_plot.py -wd summaries

结果图如下：

当然，有结果数据就可以自己做更好看的图了，不一定要用官方的。

2018年发表在Nucleic Acids Research上的一篇文章Assessing genome assembly quality using the LTR Assembly Index (LAI)，研究者提出了一种对长末端重复序列（long terminal repeats,LTRs）评估从而评价基因组完整度的方法，并且开发了对应的分析工具LTR_retriever

具体的LTR注释我会在后续的基因组注释笔记中更新，这里暂时跳过原理和背景部分，介绍下文章中提出的评估核心——LAI指数（LTR Assembly Index，LAI），也就是长末端重复序列组装指数。raw LAI = (完整LTR-RTs长度/总LTR长度)*100，修正后，LAI = raw LAI + 2.8138 × (94 – 整个基因组LTR identity)。

以下是一个完整的LTR-RTs的结构示意图：

文章还阐明LAI独立于基因组大小、LTR-RT含量以及基因空间评估指标（如BUSCO和CEGMA）等参数，可以用于鉴定低质量的基因组区域。使用这个指标要求**基因组中完整的LTR-RTs应至少占基因组0.1%且总LTR-RTs长度至少占5%**。

文章最后给出了LAI评价基因组完整度的三个指标：

2.1 LTR序列预测

LTR_retriever需要以LTR_finder和/或ltrharvest的LTR预测结果文件为输入，也可以整合两个软件的预测结果作为输入（或者其他符合格式的LTR结果文件），因此需要先安装并运行以上软件，我这里以文章中提到的软件和参数运行。

1
2
3
4
5

# LTR_finder、ltrharvest和LTR_retriever安装（ltrharvest是genometools软件的一部分）

conda install -c bioconda ltr_finder
conda install -c bioconda genometools-genometools
conda install -c bioconda ltr_retriever

LTR_finder预测LTR序列（参数均由作者给出，只有文件是自己的）：

1
2
3
4

#!/bin/bash
#SBATCH -n 10

ltr_finder -D 15000 -d 1000 -L 7000 -l 100 -p 20 -C -M 0.85 /public/home/wlxie/NextPolish/01_rundir/genome.nextpolish.fasta > baima_ltrfinder.scn

参数解释：

-D NUM Max distance between 5’&3’LTR, default is 20000 # 5’和3’LTR之间的最大距离

-d NUM Min distance between 5’&3’LTR, default is 1000

-L NUM Max length of 5’&3’LTR, default is 3500 # 5’和3’LTR最大长度

-l NUM Min length of 5’&3’LTR, default is 100

-p NUM min length of exact match pair, default is 20 # 完全匹配最小长度

-C detect Centriole, delete highly repeat regions # 检测中心粒，删除高度重复区域

-M NUM min LTR similarity threshold, default is 0.00, [0,1] #最小LTR相似度

ltrharvest预测LTR序列（ltrharvest参数均由作者给出，只有文件是自己的）：

1
2
3
4
5
6
7
8

#!/bin/bash
#SBATCH -n 10

mkdir index

gt suffixerator -db /public/home/wlxie/NextPolish/01_rundir/genome.nextpolish.fasta -indexname index/baima -tis -suf -lcp -des -ssp -sds -dna

gt ltrharvest -index index/baima -minlenltr 100 -maxlenltr 7000 -mintsd 4 -maxtsd 6 -motif TGCA -motifmis 1 -similar 85 -vic 10 -seed 20 -seqids yes > baima_ltrharvest.scn

这里要注意下要先使用genome tools里的suffixerator创建基因组索引文件，然后才可以使用ltrharvest进行LTR预测。

创建基因组索引的参数不做解释了，可以 gt suffixerator -help 查看。稍微记录下ltrharvest参数：

-minlenltr specify minimum length for each LTR，default: 100

-mintsd specify minimum length for each TSD，default: 4

-motif specify 2 nucleotides startmotif + 2 nucleotides endmotif: ****

-motifmis specify maximum number of mismatches in motif [0,3]，default: 4

-similar specify similaritythreshold in range [1..100%]，default: 85.00

-vic specify the number of nucleotides (to the left and to the right) that will be searched for TSDs and/or motifs around 5’ and 3’boundary of predicted LTR retrotransposons, default: 60

-seed specify minimum seed length for exact repeats，default: 30

-seqids use sequence descriptions instead of sequence numbers in GFF3 output，default: no

以上两个软件以同样的LTR最小相似度0.85预测LTR，得到两个结果文件baima_ltrfinder.scn和baima_ltrharvest.scn。

2.2 LAI指数计算

用上一步的输出的两个结果文件，运行LTR_retriever鉴定LTR和计算LAI指数。

1
2
3
4

#!/bin/bash
#SBATCH -n 20

LTR_retriever -genome /public/home/wlxie/NextPolish/01_rundir/genome.nextpolish.fasta -inharvest baima_ltrharvest.scn -infinder baima_ltrfinder.scn -threads 20

这一步运行结果如下：

其他文件可以后续做LTR分析用到，这里我们只要看最后一个LAI的计算结果文件：

可以看到这个结果文件中包含了整个genome和各个contig的raw LAI和LAI指数，这里就只看整个genome的LAI指数15.37，根据上面文章作者提到的分类，属于Reference级别，也就是说可以认为达到参考基因组级别。