基于模式物種的快速同源搜索軟件基準(zhǔn)測(cè)試

王殷偉,武晶菁,張宸寧,華宜家,李 鵬,嚴(yán) 潔

(南京師范大學(xué)生命科學(xué)學(xué)院,江蘇 南京 210023)

同源搜索對(duì)比較基因組學(xué)分析十分重要,是后續(xù)諸如基因組注釋、基因鑒定、基因家族聚類、GO[1]、KEGG[2]富集等一系列比較基因組學(xué)流程分析的基礎(chǔ). 以最流行的兩個(gè)直系同源推斷和基因家族聚類軟件OrthoMCL[3]和OrthoFinder[4]為例:OrthoMCL通常需要基于blastp結(jié)果進(jìn)行分析,而OrthoFinder則是內(nèi)置了blastp、MMseqs和Diamond三種搜索軟件供用戶選擇. 直系同源被定義是一種在物種形成事件后分化的特征[5],一對(duì)在不同物種中相同的基因通常會(huì)被認(rèn)定為直系同源基因. Reciprocal Best Hits(RBH)的方法因具有更少假陽(yáng)性而被廣泛用于直系同源基因的推斷[6-8]. RBH指兩個(gè)物種中通過(guò)比對(duì)搜索軟件得到最佳匹配或者最高打分的一對(duì)基因. 隨著基因組時(shí)代的來(lái)臨,蛋白數(shù)據(jù)呈現(xiàn)爆炸式的指數(shù)增長(zhǎng),而傳統(tǒng)blastp搜索越來(lái)越難以應(yīng)對(duì)快速搜索的需求,帶來(lái)了對(duì)于快速同源搜索算法和軟件開(kāi)發(fā)的需求. 各類算法的快速同源搜索軟件應(yīng)運(yùn)而生,但速度的提升往往會(huì)帶來(lái)精度上的下降. 因此,基于時(shí)間消耗、同源對(duì)、RBH以及軟件錯(cuò)誤率等指標(biāo)對(duì)于各種同源搜索軟件的評(píng)估十分重要,它決定能否找到正確的直系同源基因,決定比較基因組學(xué)分析的正確與否.

近幾年來(lái)不乏有一些針對(duì)同源搜索軟件的比較分析. 2014年Ward等[9]將blast+[10]中的blastp與last[11]、ublast[12]和blat[13]進(jìn)行了比較,結(jié)果發(fā)現(xiàn)雖然blat速度最快,但在RBH、同源對(duì)搜索上有著最低的預(yù)測(cè)數(shù),并且在RBH推斷上有著最高的錯(cuò)誤率,而ublast和last則相比blat來(lái)說(shuō),有著更多的RBH和同源對(duì)預(yù)測(cè)數(shù)以及更低的錯(cuò)誤率;隨著衡量物種親緣的指標(biāo)基因組相似度GSS(Genomic Similarity Scores)[14]的降低,各個(gè)軟件在RBH和同源對(duì)預(yù)測(cè)數(shù)上都有著隨之下降的趨勢(shì),而錯(cuò)誤率則有著隨之升高的趨勢(shì). 2016年Saripella等[15]基于16個(gè)模式物種以及結(jié)合蛋白結(jié)構(gòu)數(shù)據(jù)庫(kù)信息,對(duì)基于譜(profile-based)搜索的軟件cs-blast[16]、hhsearch[17]、phmmer[18]以及非基于譜搜索的軟件blast+中的blastp、usearch[12]、ublast和fasta[19]進(jìn)行了綜合評(píng)估,結(jié)果發(fā)現(xiàn)基于譜搜索的軟件相比非基于譜的,有著更高的AUC值,表明其精度更高,但這同樣帶來(lái)了時(shí)間消耗的巨大增加. 2020年,Hernndez等[20]對(duì)last、blast、diamond[21]和MMseqs2[22]進(jìn)行了類似于Ward等[9]的研究,在不同GSS下得到的結(jié)果和趨勢(shì)也是類似的,并且還發(fā)現(xiàn) diamond 的“very”模式在速度和RBH結(jié)果上有著良好的平衡,diamond是綜合來(lái)說(shuō)最好的軟件.

盡管有了上述的一些研究,對(duì)更多優(yōu)秀、快速同源搜索的軟件進(jìn)行比較仍然是必要的. 首先,上述的大部分軟件,都會(huì)不斷地進(jìn)行定期更新,其運(yùn)算的速率、精度是會(huì)改變的,需要進(jìn)行重新評(píng)估;其次,上述的一些研究,選擇的軟件并不全面,一些主流和新開(kāi)發(fā)的快速搜索軟件并未參與評(píng)估. 縱觀近十年的同源搜索軟件和算法的開(kāi)發(fā)和進(jìn)展,選取了具有代表性的快速搜索軟件,包括usearch/ublast、last、lambda[23]、ghostx[24]、diamond、MMseqs2以及blast在內(nèi),共8種非基于譜的軟件或算法進(jìn)行評(píng)估. 為何僅選擇非基于譜的,這很大程度上是因?yàn)?基于譜的算法和軟件,雖然帶來(lái)了精度上的提升,速度卻相比blastp有著大幅下降[15],更難實(shí)現(xiàn)大規(guī)模數(shù)據(jù)情況下的同源蛋白搜索. 本文旨在篩選出相比blastp來(lái)說(shuō)更快的算法或軟件,與此同時(shí)精度上有著更少下降或更高的替代品,以應(yīng)對(duì)大規(guī)模數(shù)據(jù)下的同源蛋白搜索.

1 材料與方法

1.1 蛋白序列收集

研究部分參照Saripella等[15]的做法,選取了15個(gè)具有代表性的、有一定跨度的模式物種(表1),涵蓋原核與真核生物,它們分別為,屬于細(xì)菌的Escherichiacoli和Staphylococcusaureus,原生動(dòng)物的Chlamydomonasreinhardtii和Dictyosteliumdiscoideum,真菌的Aspergillusnidulans和Saccharomycescerevisiae,植物的Arabidopsisthaliana和Zeamays,無(wú)脊椎動(dòng)物的Drosophilamelanogaster和Caenorhabditiselegans以及脊椎動(dòng)物的Homosapiens、Daniorerio、Xenopustropicalis、Gallusgallus和Musmusculus,各自從NCBI基因組數(shù)據(jù)庫(kù)中下載對(duì)應(yīng)基因組序列和gff注釋,結(jié)合基因組序列和注釋信息,提取蛋白序列,而對(duì)于有著不同可變剪切轉(zhuǎn)錄本的基因,則保留最長(zhǎng)的蛋白序列作為該基因的代表,因而得到對(duì)于每個(gè)物種來(lái)說(shuō)都是非冗余的蛋白序列.

1.2 實(shí)驗(yàn)機(jī)器、軟件運(yùn)行以及時(shí)間與內(nèi)存消耗

本研究采用軟件usearch/ublast v11.0.667_i86linux32、last 1256、lambda2 v1.9.5、ghostx v1.3.6、diamond v2.0.6.144、MMseqs2 Release 13-45111和blast v2.5.0,以及部分軟件的不同精度,即diamond-fast、diamond-sensitive、diamond-more、diamond-very、diamond-ultra、MMseqs-s3、MMseqs-s5和MMseqs-s7.5,基于一臺(tái)系統(tǒng)為Centos8的中小型服務(wù)器以4線程以及1e-6的閾值設(shè)定運(yùn)行(表2),物種的蛋白集合兩兩比對(duì),并且包括自身比對(duì). 服務(wù)器CPU型號(hào)為AMD Ryzen Threadripper 3970X 32-Core Processor,一共32核,每核兩線程,內(nèi)存總大小為120GB,使用unix命令“time”計(jì)算并記錄每次運(yùn)行的真實(shí)時(shí)間(real times),在unix下的“while”循環(huán)中檢測(cè)監(jiān)測(cè)私有內(nèi)存與共有內(nèi)存的總即時(shí)消耗,程序運(yùn)行完成后計(jì)算平均內(nèi)存消耗,每次運(yùn)行時(shí)確保無(wú)其他任務(wù)運(yùn)行占用計(jì)算額外資源導(dǎo)致時(shí)間計(jì)算出現(xiàn)偏差.

表1 選擇的15個(gè)模式物種基因組序列信息Table 1 Genome sequence information for candidate 15 model species

表2 用于運(yùn)行軟件的命令行Table 2 Command lines used to run each program

1.3 同源對(duì)和RBH的統(tǒng)計(jì)計(jì)算

同源對(duì)數(shù)目即每個(gè)軟件運(yùn)行得到的結(jié)果數(shù)目即hits數(shù),而RBH則是每次運(yùn)行結(jié)果中每個(gè)Query對(duì)應(yīng)最佳的匹配. 首先繪制了每個(gè)軟件或不同精度對(duì)應(yīng)的箱線圖,然后根據(jù)blastp結(jié)果計(jì)算基因組相似度得分,即GSS. 在這里,計(jì)算了每對(duì)物種包括物種本身的GSSa[14],并且以相對(duì)于blastp的時(shí)間、同源對(duì)數(shù)目、RBH分別與GSS作誤差棒圖. 還計(jì)算了每種具有共線性順序支持的可信RBH數(shù)目,具體來(lái)說(shuō),如果一個(gè)物種相鄰的兩個(gè)基因a和b,與另一個(gè)物種相鄰的基因a′、b′相比,a和a′為RBH,b和b′為RBH,那么這兩個(gè)RBH都為真陽(yáng)性結(jié)果[7,25-26],在這里稱之為可信RBH,在本文暫稱之為CRBH(Credible RBH),同樣,也繪制了CRBH與GSS的誤差棒圖,最后繪制了6種有較好表現(xiàn)軟件結(jié)果的韋恩圖,所有繪圖、統(tǒng)計(jì)計(jì)算均在Python 3.6下進(jìn)行,以Matplotlib模塊繪制圖片.

1.4 計(jì)算錯(cuò)誤發(fā)現(xiàn)率

為了評(píng)估不同模型或程序得到結(jié)果的準(zhǔn)確性,首先用InterproScan v5.8[27]軟件的superfamily v1.75[28]蛋白家族結(jié)構(gòu)數(shù)據(jù)庫(kù)對(duì)所有蛋白進(jìn)行了注釋,接下來(lái)對(duì)所有程序運(yùn)行結(jié)果在不同期望閾值下的錯(cuò)誤發(fā)現(xiàn)率進(jìn)行統(tǒng)計(jì)計(jì)算. 具體來(lái)說(shuō),如果程序搜索匹配得到的一對(duì)蛋白,有著完全相同的superfamily注釋結(jié)果,則為陽(yáng)性匹配;有著部分相同的superfamily注釋結(jié)果,則為模糊匹配,不參與后續(xù)計(jì)算;有著完全不同的superfamily注釋,則為陰性匹配,錯(cuò)誤發(fā)現(xiàn)率則等于在不同期望閾值下的假陽(yáng)性結(jié)果總數(shù)除以所有陽(yáng)性結(jié)果的總數(shù).

2 結(jié)果與討論

2.1 運(yùn)行時(shí)間

將所有程序運(yùn)行計(jì)算得到的時(shí)間與blastp相除,來(lái)觀測(cè)其各自相對(duì)blastp所節(jié)省的時(shí)間. 結(jié)果發(fā)現(xiàn),大部分程序運(yùn)行時(shí)間平均數(shù)不到blastp時(shí)間消耗的5%,而last、usearch和diamond-fast運(yùn)行時(shí)間平均數(shù)則不到blastp時(shí)間消耗的2.5%,是速度最快的3個(gè)程序(圖1). 雖然從圖中看似并無(wú)太大差異,但配對(duì)樣本t檢驗(yàn)表明,3個(gè)程序在時(shí)間消耗節(jié)省上來(lái)說(shuō),diamond-fast

除了關(guān)注不同軟件間的比較,還關(guān)注了相同軟件不同精度下的運(yùn)行效率,即diamond和MMseqs2的不同精度模式下的表現(xiàn).

對(duì)于diamond來(lái)說(shuō),研究發(fā)現(xiàn),與之前的研究類似,“sensitive”、“more”和“very”的運(yùn)行時(shí)間幾乎無(wú)太大差異[20],因此,在選擇這三種選項(xiàng)時(shí),理論上來(lái)說(shuō)一定是精度更高的“very”模式更好,而“fast”相比這三種精度下的運(yùn)行速度明顯更快,“ultra”相比則更慢.

對(duì)于MMseqs2來(lái)說(shuō),不同的精度下,速度有著明顯的差異,尤其是s7.5精度模式,時(shí)間消耗已經(jīng)高于除了blastp外的所有程序,并且在一些運(yùn)行中,速度慢于blastp,這體現(xiàn)了其運(yùn)行時(shí)間的巨大變異性和不穩(wěn)定性,如果該精度下后續(xù)的評(píng)估中并沒(méi)有發(fā)現(xiàn)隨之帶來(lái)的各種指標(biāo)評(píng)估下的良好改進(jìn),那么該精度模式則是個(gè)不被推薦的選項(xiàng).

本研究并未發(fā)現(xiàn)不同程序時(shí)間消耗與GSS之間的關(guān)聯(lián)性(圖2),但可以發(fā)現(xiàn)的是,MMseqs的s7.5精度模式的相對(duì)耗時(shí)在不同GSS下波動(dòng)較大,呈現(xiàn)“中間低,兩頭高”的模式,diamond的“ultra”模式也呈現(xiàn)類似模式,但較前者更加平穩(wěn)些,其他軟件之間的波動(dòng)以及趨勢(shì)并不完全一致,但總體呈現(xiàn)在不同GSS下的平穩(wěn)均勻分布,這說(shuō)明大部分軟件每次實(shí)驗(yàn)相對(duì)于blastp節(jié)省時(shí)間的比例是一個(gè)穩(wěn)定的小區(qū)間,換言之,最不穩(wěn)定的MMseqs的s7.5精度模式,如果在精度上沒(méi)有表現(xiàn)出巨大優(yōu)越性,那么其在大規(guī)模同源搜索時(shí)就不是一個(gè)良好的選擇.

圖1 不同程序每次同源搜索相對(duì)blastp運(yùn)行速度的差異Fig.1 Differences in the speed of homologous protein searchrelative to blastp by different programs

圖2 不同程序在不同GSS下的相對(duì)于blastp的時(shí)間消耗的誤差棒圖Fig.2 Error bar graph of time consumption relative to blastpfor different programs at different GSS

2.2 運(yùn)行內(nèi)存

同樣如上,將所有程序計(jì)算得到的平均內(nèi)存消耗與blastp相除,得到相對(duì)blastp的內(nèi)存消耗. 結(jié)果顯示(圖3),除了MMseqs2,其余所有程序的相對(duì)運(yùn)行內(nèi)存都顯著高于blastp,其中g(shù)hostx最為顯著,其平均運(yùn)行內(nèi)存消耗約為blastp的22倍,這表明在運(yùn)行g(shù)hostx進(jìn)行同源搜索時(shí),尤為需要注意可用內(nèi)存空間的大小. 其次比較高的是diamond的“ultra”精度模式,而usearch、ublast、last和lambda則有著比blastp較高但接近的內(nèi)存占用. 值得注意的是,這些結(jié)果在大體上與時(shí)間消耗表現(xiàn)一致,這可能表明有些內(nèi)存占用較少的軟件可能不需要太大的開(kāi)銷進(jìn)行更多的搜索,從而在時(shí)間消耗上表現(xiàn)較低.

當(dāng)關(guān)注到相同軟件不同精度下的內(nèi)存占用時(shí),MMseqs2三種精度模式下相對(duì)平均內(nèi)存消耗的分布并無(wú)太大差異,但明顯低于其他程序,包括blastp(配對(duì)樣本t檢驗(yàn),p-value<1e-5),而diamond則隨著精度的提高,呈現(xiàn)出明顯的內(nèi)存消耗提升.

2.3 同源蛋白對(duì)

與時(shí)間消耗誤差棒圖(圖2)一致的是,研究發(fā)現(xiàn),不同軟件在不同GSS下鑒定出的相對(duì)于blastp的同源蛋白數(shù)的數(shù)目分布也呈現(xiàn)“中間低,兩頭高”的趨勢(shì)(圖4),并且總體分布趨勢(shì)與時(shí)間消耗圖吻合,這表明程序能夠鑒定出的結(jié)果數(shù)目的大小和時(shí)間是有一定關(guān)聯(lián)性的,這也與預(yù)期和直覺(jué)一致,結(jié)果數(shù)目越多,時(shí)間消耗越大.

幾乎所有軟件在不同GSS下鑒定出的同源蛋白對(duì)都要少于blastp的結(jié)果. 唯一的例外是ghostx在低GSS的情況下獲得的結(jié)果大大增加,最多可達(dá)blastp結(jié)果的3倍以上. 而ublast則在很高GSS的情況下,鑒定的結(jié)果數(shù)目越為接近blastp,但大部分情況下,MMseqs的s7.5和s5模式有更多的結(jié)果數(shù). 可以看到MMseqs的s7.5精度模式在時(shí)間消耗上的提升確實(shí)帶來(lái)了鑒定同源蛋白對(duì)數(shù)目上的提升,但提升并不明顯:在低GSS下顯著低于ghostx,在高GSS下的則略少于ublast. 如果僅從同源蛋白結(jié)果數(shù)目上來(lái)看,ghostx適用遠(yuǎn)緣搜索,ublast適用于近緣搜索. 在中等GSS的情況下,MMseqs的s5精度模式獲得結(jié)果的數(shù)目?jī)H次于s7.5精度模式,但考慮到時(shí)間消耗,并且s5與s7.5的數(shù)目差異并不太大,MMseqs的s5精度模式是更好的選擇.

圖3 不同程序每次同源搜索相對(duì)blastp運(yùn)行內(nèi)存消耗的差異Fig.3 Differences in memory consumption of homologousprotein search relative to blastp by different programs

圖4 不同程序在不同GSS下鑒定出的相對(duì)于blastp的同源蛋白數(shù)的誤差棒圖Fig.4 Error bars of the number of homologous proteins identifiedby different programs relative to blastp at different GSS

2.4 RBH

除了考慮同源蛋白數(shù)目上的評(píng)估,對(duì)于RBH的評(píng)估也是非常重要的,因?yàn)榍罢哧P(guān)聯(lián)著基因鑒定,而后者則關(guān)聯(lián)著直系同源推斷,兩者都是比較基因組分析的重要步驟. 不同于相對(duì)時(shí)間消耗以及同源蛋白與GSS不明顯的關(guān)系模式,不同程序鑒定出的RBH相對(duì)于blastp鑒定出的數(shù)量,在大部分軟件中都呈現(xiàn)出了隨著GSS降低而降低的趨勢(shì)(圖5),而MMseqs的s7.5精度模式,則在所有GSS下與blastp數(shù)目保持一致和穩(wěn)定. 令人驚訝的是,diamond的“ultra”、“sensitive”、“more”和“very”的RBH數(shù)目,盡管在低GSS下有所波動(dòng)下降,但整體都穩(wěn)定在blastp結(jié)果數(shù)目的約90%左右,而ghostx則呈現(xiàn)出隨著GSS降低相對(duì)RBH數(shù)顯著增高的趨勢(shì),這同樣也表現(xiàn)出了ghostx在遠(yuǎn)緣搜索的相對(duì)優(yōu)勢(shì). 結(jié)合相對(duì)時(shí)間消耗來(lái)考慮,這些結(jié)果表明,在進(jìn)行遠(yuǎn)緣搜索進(jìn)行直系同源推斷的時(shí)候,在不考慮錯(cuò)誤率的情況下,ghostx是一個(gè)良好的選擇,可以獲得最多的結(jié)果數(shù)目,而在其他情況下,考慮到MMseqs7.5的耗時(shí)之多,以及diamond除了“fast”外其他精度結(jié)果數(shù)目的接近以及“sensitive”、“more”和“very”時(shí)間消耗的接近,diamond的“very”仍然是一個(gè)綜合來(lái)說(shuō)非常優(yōu)秀的運(yùn)行模式和優(yōu)先考慮的選擇.

然而,單純從數(shù)量上來(lái)評(píng)估RBH鑒定及直系同源推斷能力的優(yōu)秀與否是不可行的,因?yàn)楦嗟腞BH有可能引入更多的假陽(yáng)性結(jié)果,因此需要對(duì)RBH評(píng)估的錯(cuò)誤率進(jìn)行推斷,之前的兩個(gè)研究都是考慮共線性關(guān)系,以及旁系同源關(guān)系,來(lái)計(jì)算錯(cuò)誤率[9,20].

然而,不同軟件鑒定的旁系同源基因,仍然有假陽(yáng)性的可能,這里考慮以CRBH進(jìn)行比較,即僅比較相對(duì)blastp來(lái)說(shuō),具有相鄰共線性位置關(guān)系支持的可信RBH數(shù)目,因而避免其他噪聲. 結(jié)果發(fā)現(xiàn),在高GSS的情況下,各程序CRBH數(shù)目與blastp結(jié)果數(shù)目差別不大,而隨著GSS的降低,各程序的相對(duì)CRBH出現(xiàn)了明顯的分歧(圖6). ghostx同樣在CRBH上,體現(xiàn)著其在遠(yuǎn)緣搜索的優(yōu)勢(shì),這表明,ghostx得到的RBH、同源蛋白數(shù)目在遠(yuǎn)緣中的增加,同時(shí)確實(shí)會(huì)帶來(lái)真實(shí)的、可靠的RBH結(jié)果的增加,而其次優(yōu)秀的就是diamond的“ultra”、“sensitive”、“more”和“very”模式,在GSS降低的情況下,也展現(xiàn)出了比blastp更多的結(jié)果,并且彼此之間的差異不大,這說(shuō)明,從CRBH上來(lái)考慮,diamond的“very”精度模式仍然是一個(gè)良好的選擇.

圖5 不同程序在不同GSS下鑒定出的相對(duì)于blastp的RBH數(shù)的誤差棒圖Fig.5 Error bars of the number of RBH identified by differentprograms relative to blastp at different GSS

圖6 不同程序在不同GSS下鑒定出的相對(duì)于blastp的CRBH數(shù)的誤差棒圖Fig.6 Error bars of the number of CRBH identified by differentprograms relative to blastp at different GSS

2.5 代表性軟件交集評(píng)估

對(duì)于各個(gè)軟件的結(jié)果重合情況評(píng)估,也是非常重要的,越多的重合也同樣表明了結(jié)果的可靠性. 為了使得繪圖具有可讀性,并且考慮到usearch和last在除了時(shí)間消耗之外其他指標(biāo)上的評(píng)估表現(xiàn)較差以及MMseqs2和diamond的更高精度理論上應(yīng)有更好的結(jié)果表現(xiàn),僅選取了blastp、diamond的“ultra”精度模式、MMseqs2的s7.5精度模式、ublast和lambda作為代表,進(jìn)行評(píng)估.

首先統(tǒng)計(jì)、繪制同源蛋白對(duì)結(jié)果的韋恩圖(圖7),如圖所示,大部分軟件之間所鑒定出的同源蛋白對(duì)之間都有著或多或少的重合,所有軟件的重合僅有4.8%,大部分軟件都與blastp之間有著重合,這些重合結(jié)果具有可靠性. ghostx,則有著22.0%的最多特有搜索結(jié)果,這表明其與其他軟件搜索結(jié)果交集比較少,原因正如前文所述,其鑒定出了更多的同源蛋白對(duì)、更多的RBH、更多的CRBH;而diamond的“ultra”精度模式僅有0.1%的特有搜索結(jié)果,這表明其大量的搜索結(jié)果與其余軟件都有交集,表明了其結(jié)果的可靠性,diamond軟件的優(yōu)秀之處又在此處展現(xiàn). ublast和lambda這些更加快速但精度相較來(lái)說(shuō)更低、搜索結(jié)果更少的軟件,特有的結(jié)果分別僅有5.6%和3.0%,這可能表明,這些快速搜索軟件能夠快速搜索出一些與其他軟件重合的、可靠的同源搜索結(jié)果,然后快速結(jié)束搜索,保留下少部分可靠的結(jié)果,至少在其鑒定出的結(jié)果上來(lái)說(shuō),是比較可靠的.

接下來(lái)繪制了不同代表性軟件的RBH結(jié)果的韋恩圖(圖8),可以明顯看出,相比同源蛋白鑒定結(jié)果在各個(gè)軟件中的差異性,RBH結(jié)果的差異性明顯更少,所有軟件的共有的RBH占到了所有結(jié)果的 21.7%,這同樣也是可以預(yù)期的,因?yàn)檎鎸?shí)的直系同源對(duì)相較于其他的同源基因?qū)?往往會(huì)有著最高的相似性,從而在軟件搜索中獲得最高得分而被保留. diamond的“ultra”精度模式和MMseqs2的s7.5模式所擁有的特有RBH最少,約為1.2%,但是仍然需要速度上的考量,diamond在快速預(yù)測(cè)出可靠的、準(zhǔn)確的RBH上具有優(yōu)勢(shì). 同樣符合預(yù)期的是,在遠(yuǎn)緣預(yù)測(cè)以及RBH數(shù)量、CRBH數(shù)量上具有明顯優(yōu)勢(shì)的ghostx,具有顯著最多的特有RBH預(yù)測(cè)占比,達(dá)到了約43.4%,lambda和ublast在特有的搜索結(jié)果上占比同樣較低,分別有5.4%和3.8%.

總而言之,6種代表性軟件統(tǒng)計(jì)、繪制的韋恩圖與預(yù)期相符,ghostx在遠(yuǎn)緣搜索上能夠鑒定出更多的同源蛋白對(duì)、RBH和CRBH,因此會(huì)有著更多的特有鑒定結(jié)果,而diamond軟件鑒定結(jié)果可靠性,也體現(xiàn)在其與各個(gè)軟件結(jié)果都互有交集,特有結(jié)果數(shù)目較少上.

圖7 代表性軟件鑒定的同源蛋白韋恩圖Fig.7 Venn diagram of homologous proteins identified byrepresentative software

圖8 代表性軟件鑒定的RBH韋恩圖Fig.8 Venn diagram of RBH identified byrepresentative software

圖9 不同程序搜索結(jié)果的錯(cuò)誤發(fā)現(xiàn)率Fig.9 False discovery rate of search results for different programs

2.6 錯(cuò)誤發(fā)現(xiàn)率

為了對(duì)不同軟件程序的準(zhǔn)確性進(jìn)行更進(jìn)一步的評(píng)估,根據(jù)superfamily數(shù)據(jù)庫(kù)對(duì)所有搜索蛋白的注釋結(jié)果,計(jì)算了不同期望閾值下的錯(cuò)誤發(fā)現(xiàn)率(圖9). 結(jié)果發(fā)現(xiàn),ghostx、MMseqs2的s5、s7.5以及diamond所有精度下的錯(cuò)誤發(fā)現(xiàn)率均低于0.01,blastp的則低于0.02,lambda、ublast和usearch在低期望閾值(小于 1e-10)下能夠?qū)㈠e(cuò)誤發(fā)現(xiàn)率總體控制在0.05以下,而last和MMseqs2的s3的錯(cuò)誤發(fā)現(xiàn)率比較高,last總體在0.1以下但高于0.08,MMseqs2的s3則表現(xiàn)最差,總體在0.14以下但高于0.12,這可能表明s3模式相對(duì)來(lái)說(shuō)要盡量避免使用.

值得注意的是,雖然從直覺(jué)上來(lái)說(shuō)時(shí)間消耗較少、精度模式選擇較低的軟件應(yīng)該有著很高的錯(cuò)誤率,但實(shí)際上并非如此,diamond和MMseqs2的不同精度模式就是很好的例子. diamond的“fast”精度模式下的錯(cuò)誤發(fā)現(xiàn)率是最低的,這表明其雖然進(jìn)行快速搜索后,在數(shù)量上要少于高精度的搜索,但其至少能保證得到的結(jié)果是準(zhǔn)確的,這與MMseqs的s3模式是相對(duì)的. blastp、ghostx、MMseqs2和diamond高精度的搜索,在獲取的結(jié)果數(shù)量上高于其他快速搜索軟件的同時(shí),也能將錯(cuò)誤發(fā)現(xiàn)率控制在較低的值,可見(jiàn)這些軟件的可靠性和優(yōu)秀性.

3 結(jié)論

本研究選取了共7種快速同源搜索軟件或程序usearch/ublast、last、lambda、ghostx、diamond、MMseqs2,包括diamond的5種不同精度模式“fast”、“sensitive”、“more”、“very”、“ultra”以及MMseqs2的3種不同精度模式s3、s5、s7.5,與blastp在時(shí)間消耗、同源蛋白對(duì)、RBH、CRBH、重合狀況以及錯(cuò)誤發(fā)現(xiàn)率上進(jìn)行綜合比較,來(lái)選擇在不同狀況下對(duì)大數(shù)據(jù)進(jìn)行搜索的blastp的替代品. 結(jié)果表明,如果追求速度同時(shí)保證準(zhǔn)確性,diamond的“very”精度模式是最佳選擇,因?yàn)槠溆兄畹偷腻e(cuò)誤發(fā)現(xiàn)率以及最快的搜索速度,而MMseqs的s3精度模式則有著最高的錯(cuò)誤發(fā)現(xiàn)率,可能需要避免使用;在進(jìn)行遠(yuǎn)緣物種同源搜索、直系同源推斷時(shí),ghostx由于能夠得到更多的直系同源對(duì)、RBH、CRBH以及適中的時(shí)間消耗節(jié)省和非常低的錯(cuò)誤發(fā)現(xiàn)率,成為最佳的選擇,盡管其有著最高的內(nèi)存消耗;而只有在進(jìn)行近緣物種搜索時(shí),快速搜索軟件ublast能得到更多的同源蛋白對(duì)結(jié)果,在大部分的GSS下,對(duì)于同源蛋白的搜索與鑒定,MMseqs7.5與MMseqs5差異不大,且兩者錯(cuò)誤發(fā)現(xiàn)率均低于0.01,都是良好的兩個(gè)選擇,但考慮到時(shí)間消耗,MMseqs5應(yīng)是更好的選擇;如果是應(yīng)對(duì)于進(jìn)行直系同源推斷的研究目的需求,diamond是綜合來(lái)說(shuō)最佳的軟件,并且其最適合以“very”精度選項(xiàng)運(yùn)行,能夠得到速度與精度良好的權(quán)衡. 本研究為不同目的下選擇和使用不同快速搜索軟件提供了參考和指南.