Genome assembly assessment first identifies candidate regions to be assessed with tBLASTn searches using BUSCO consensus sequences. Gene structures are then predicted using Augustus with BUSCO block profiles. These predicted genes, or all genes from an annotated gene set or transcriptome, are then assessed using HMMER and lineage specific BUSCO profiles to classify matches. The recovered matches are classified as ‘complete’ if their lengths are within the expectation of the BUSCO profile match lengths. If these are found more than once they are classified as ‘duplicated’. The matches that are only partially recovered are classified as ‘fragmented’, and BUSCO groups for which there are no matches that pass the tests of orthology are classified as ‘missing’.                    --BUSCO manual

BUSCO流程框架如下图：

开始时BUSCO载入程序配置文件config.ini，这个文件记录一些默认参数和软件路径。如果run_BUSCO.py指定的是genome模式，则运行如下：

1 tblastn

按照源代码，比对参数如下：

makeblastdb -in xx.fasta -dbtype nucl -out temp.file

tblastn

    -query 1e-03

    -num_threads 30

    -query xx.fasta

    -db temp.file

    -out blast_output/tblastn_xx.tsv

    -outfmt 7

#注意：其中-query的参数可以通过run_BUSCO.py的--evalue指定；-num_threads可以通过run_BUSCO.py --cpu指定；-out中会用到run_BUSCO.py的--out指定；

生成的文件如下：

注：BUSCO会读取这个文件最后一行是否有"processed"字样，以自动评估这个文件是否正常生成

再解析tblastn output表格，先按run_BUSCO.py --limit选项(默认为3，后面讲解--limit=3)载入每个BUSCO比对到的前3条记录。当载入第4个记录且这个记录的contig不在上述的列表中，而且当这个记录中evalue更小时，则这个记录进入前三(通常tblastn是按bitscore排序的，然而BUSCO会按evalue值排序)。因为blast是局部比对，如果这个记录中contig在前3中出现过，会尽量做一个merge操作，具体如下fragement1已经在前三，fragment2是新的记录，它们会按最长50kb合并(更多的时候是fragment1与fragment2之间有gap区)

完成后比对长度还会做一个筛选。先迭代每个比对区间，检测最大的比对长度max_size。再迭代一次，当某个区间的比对长度小于max_size*0.7时则去除这个比对区间。

运行完结果会记录在blast_output/coordinates_xx.tsv，这个文件中每个BUSCO对应的contigs数最多为3。格式如下，第一列是库里的BUSCO ID，第二列是contig名，第三列是比对区间的起始位点，第四列是比对区间的终止位点。这个文件作为augustus的输入

这一步是采用BuscoAnalysis类中的_get_coordinates()方法完成的

2 augustus

先用augustus找出tblastn比对出的contigs片段上可能的编码区，并将其转化为protein序列，这个序列作为hmmersearch的输入。按照源码参数如下：

/path/to/augustus-3.3.1/bin/augustus

    --codingseq=1

    --proteinprofile=xx.protein.fasta

    --predictionStart=xx

    --predictionEnd=xx

    --species=xx

augustus的运行参数不多，species可以通过以下查看可选的物种

/path/to/augustus-3.3.1/bin/augustus --species=help

运行完，BUSCO清掉文件大小为0的空文件，再没有做其它清理工作

3 hmmersearch

按照源码，具体扫描参数如下：

/path/to/hmmer-3.1b1/bin/hmmsearch

     -o ./tmp/temp_all_ctg_1852078540

     --domtblout /path/to/workdir/hmmer_output/EOG09370VLM.out.3

     --cpu 1

     /path/to/eukaryota_odb9/hmms/EOG09370VLM.hmm

     /path/to/workdir/augustus_output/extracted_proteins/EOG09370VLM.faa.3

#其中extracted_proteins faa序列即augustus生成的蛋白序列。

输出的--domtblout结果格式如下。

获取到hmmersearch结果后，程序会进行一些过滤，主要是在通过BuscoAnalysis.py文件中_parse_hmmer()方法完成的。具体步骤如下：

(1) run_BUSCO.py需要指定lineage信息，在这个lineage目录下有个scores_cutoff文件，它是每个BUSCO算法训练出的最小的bitscore阈值。首先按照这个值对hmmersearch结果进行过滤。

(2) 构建一个字典，这一步看源码很重要，如果不看源码可以跳过，这里没过滤步骤

先构建一个如下图的字典，其中黄色线第一个是BUSCO ID，第二个是ctg id(augustus转protein的)，第三个是hmmersearch生成的文件；后面的蓝色框内的部分是list-of-list，第1个元素是[hmm_start1, hmm_end1, bit_score1, total_length, sigma, hmmersearch_outfile]，后面的元素是[hmm_start2, hmm_end2, bit_score2]，接下来只会用到第1个列表元素

#注意：其中，

=> total_length是这段protein上的search到这个BUSCO的总长度，如上，

(68-3) + (101-66) + (115-97) + (179-139) = 158

=> sigma是BUSCO最核心的地方，它算出来的这个sigma值用于评估是否完整(即complete还是fragment)

(3) 过滤掉识别为两个或多个BUSCOs 的基因 (identify genes belonging to two or more buscos, keep the best score)  => 一个基因对多个BUSCOs时

源码中称augustus检测出来的序列为"gene"，例如g4[000010F_arrow_pilon:804424-805324]。这里是先遍历每个基因，找到bitscore最高的那个gene保留下来，所以一个基因只会对应一个BUSCO (注意：因为contigs上可能有多个基因，因此一个contigs可能对应多个BUSCO)。

# 这里似乎有点怪异，只是用了bitscore值没有涉及total length！因为一方面只用了上图中字典的第1个元素获取bitscore，另一方面按照代码本意应该是会涉及total_length，否则专门计算出来，但是又完全不使用。

这步是采用BuscoAnalysis类中的_filter_multi_match_genes()方法完成的，源码中的注释可以帮助理解:

This function identifies genes that match the same BUSCO, and keeps the one with the best score only

(4) 比对到同一个BUSCO的基因中，过滤掉比对得分低的基因。=> 一个BUSCO对多个基因时

如果有多个gene通过hmmersearch对应到相同的BUSCO，如下。这时先找出最高的bitscore值max_score，再迭代一次，找出bitscore值不小于max_score*0.85的这些hmmersearch结果。如果其数量大于2，就产生了Duplicated BUSCOs，否则即为single-copy BUSCOs。

'EOG09370YC4': {

'g5[000048F_arrow_pilon:78932-80517]-run_all_ctg/hmmer_output/EOG09370YC4.out.2': 

           [28,160, 60.8,132,3, 

            'run_all_ctg/hmmer_output/EOG09370YC4.out.2', 2,  

            'run_all_ctg/hmmer_output/EOG09370YC4.out.3', -1, 

            'run_all_ctg/hmmer_output/EOG09370YC4.out.1'],

'g6[000007F_arrow_pilon:176764-178195]-run_all_ctg/hmmer_output/EOG09370YC4.out.1': 

           [13, 161, 66.5, 148, 1, 

            'run_all_ctg/hmmer_output/EOG09370YC4.out.1']

}, ...

这一步采用BuscoAnalysis类中的_remove_bad_ratio_genes()方法完成的，源码的注释可以帮助理解：

This function removes duplicate positive results if the score is above a 0.85 threshold compared to the top scoring match.

参考材料：

[1] https://busco.ezlab.org/busco_userguide.html

[2] Mathieu Seppey, Mosè Manni, Evgeny M Zdobnov. BUSCO: Assessing Genome Assembly and Annotation Completeness. 2019. Gene Prediction: Methods and Protocols, Methods in Molecular Biology.