Easy353使用,快速簡便提取被子植物353基因集
置頂注釋:建議先閱讀我最下面的?24.1.14 加筆,看完后再考慮看不看其他正文內(nèi)容
被子植物353基因集(Angiosperms353 gene set, AGS) 是這樣一組通用的低拷貝基因集合。川大余巖老師團(tuán)隊開發(fā)的easy353程序可以高效地從高通量測序數(shù)據(jù)中過濾reads并從頭組裝,提取出AGS353序列。
運行原理簡要如圖,詳細(xì)內(nèi)容可閱讀原文獻(xiàn)進(jìn)行了解:
https://doi.org/10.1093/molbev/msac261

了解、下載easy353:https://github.com/plant720/Easy353
用相同方法挖掘更多基因,可以使用GeneMiner擴展版: https://github.com/happywithxpl/GeneMiner
Easy353有兩種版本:命令行操作版本和gui交互窗口版本。
窗口界面易于操作,本著能不麻煩絕不去找事的原則,選擇使用這個版本。提供了win和macOS,根據(jù)自己的操作系統(tǒng)按需下載,解壓即可使用。https://github.com/plant720/Easy353/release

參數(shù)界面分了三個部分:Basic、General、Advance三個部分,下面逐一介紹。

Basic
Unpaired fq file:輸入單端測序的結(jié)果。
Paired fq file 1/2:分別輸入雙端測序的結(jié)果文件。
(上面兩個根據(jù)測序方式選擇輸入)
Taxonomy:參考序列類群范圍。(選擇關(guān)系最近的類群,否則很難得到結(jié)果;參考序列會自動下載,數(shù)據(jù)來自邱園的項目https://treeoflife.kew.org/,沒有網(wǎng)絡(luò)的話會報錯)
Output dir:設(shè)置輸出文件夾。

General
Exclude file和Exclude:用于指定不使用哪些參考物種序列。(一般不用)
Filtering Kmer:設(shè)定過濾時的Kmer長度。越長就越難匹配到參考序列上,但是太短了匹配上的reads就會非常多,會出現(xiàn)錯誤。默認(rèn)為31。
Assembly Kmer:設(shè)定組裝的Kmer長度。默認(rèn)為41。
(測的乘數(shù)較高的話默認(rèn)參數(shù)應(yīng)該足夠了,推薦的乘數(shù)為20×以上,乘數(shù)不夠或者沒有較近關(guān)系的類群參考,可以適當(dāng)降低Kmer設(shè)定,過濾Kmer不應(yīng)低于21,組裝Kmer不應(yīng)短于31,越低結(jié)果錯誤的風(fēng)險越大)
Filtering thread:過濾運行線程,gui版本為1,電腦能雙線程的話可以嘗試,不過要么沒用要么很卡,反正我沒試過。
Assembly thread:組裝運行線程,默認(rèn)為4。
Advanced
Step length:獲取Kmer的步長,比如一個reads為AATTCCGG,我設(shè)定Kmer長度為5,步長為1的話得到以下Kmer:AATTC、ATTCC、TTCCG、TCCGG,設(shè)定步長為2的話就是:AATTC、TTCCG。舉的例子可能不好,但大致就是這個意思,步長短獲取的Kmer就多,約容易獲取匹配,如果測序覆蓋度夠高可以適當(dāng)調(diào)高來減少運行時間。
Ref number:當(dāng)參考物種較多時,設(shè)定轉(zhuǎn)化為哈希表的最大參考數(shù),用于減少運行時間。
Change seed: The setting for the number of seed changes, default value is 32. Actually, change seed is the amount of times the assembly's beginning point can be changed. The seeds are high-abundance K-mers selected from filtered reads that serve as the beginning point for de novo assembly. When the assembled gene's length is less than the set value, Easy353 will alter the assembly beginning point.(不能理解,但感覺對結(jié)果影響不是很大,默認(rèn)值應(yīng)該就OK了)
Kmer limit:用于刪除豐度低或錯誤的Kmer,只有出現(xiàn)次數(shù)大于限定設(shè)置的Kmer才回用于組裝,如果測的數(shù)據(jù)集覆蓋率很高,可以使用更大的值來保證質(zhì)量。
Min/Max length ratio:組裝基因與參考基因長度比的上下限設(shè)定。默認(rèn)為0.8-2.0范圍。
輸出會有三個文件夾
353_genes_of_species:從邱園的官網(wǎng)上下載的參考物種序列。
easy353_result:其中filtered_reads_for_genes為過濾后的每個基因?qū)?yīng)的Kmer;而target_genes包含組裝結(jié)果,組裝不成功的基因會放置在單獨的文件夾中,assemble_log.csv文件包含結(jié)果記錄。
reference_of_353_genes:將參考物種序列按照基因編號進(jìn)行分組。
?
easy353_result文件重命名備份,另兩個文件夾作是參考序列,可以保留在工作目錄中,以免同樣的內(nèi)容反復(fù)下載浪費流量。

23.11.15加筆
GeneMiner和其gui模式也可以從github獲取:https://github.com/happywithxpl/GeneMiner/releases
GeneMiner運行需要Python環(huán)境,首先要保證你的服務(wù)器安裝了Python,安裝后要確認(rèn)一下是不是加入了自己的環(huán)境變量或全局環(huán)境變量。
GeneMiner安裝和具體使用說明參考:https://github.com/happywithxpl/GeneMiner
下面僅做示例
輸出如下

使用的時候可能會報錯找不到 'Bio',這時可以安裝一個Biopython包
23.11.21 加筆
GeneMiner的內(nèi)存利用率似乎偏低,在組裝那一步會占用大量內(nèi)存,同時跑過多線程容易導(dǎo)致服務(wù)器死機(快被機房管理老師記熟了……)
雖然可以控制數(shù)據(jù)使用量來降低內(nèi)存消耗,但是可能會降低輸出質(zhì)量。
服務(wù)器提前安裝earlyoom,防止死機
bootstrap可以用來篩一遍,內(nèi)存占用極大,可以用,但沒必要。篩完剩的可能每個樣品又有點差別,最后矩陣?yán)锞筒皇讉€了,剩下的未必就對的齊。

23.11.23 加筆
GeneMiner現(xiàn)在暴露出一個缺陷,具有內(nèi)含子的基因挖出來會帶著內(nèi)含子,因此以轉(zhuǎn)錄組做參考得到的目標(biāo)序列會在中間多出一段或多段。
內(nèi)含子跟外顯子的進(jìn)化速率一般來說是不同的,然而GeneMiner得到的結(jié)果也無法直接用于分區(qū)分析,只能通過人工校對來手動剪切,對于大數(shù)據(jù)量來說是十分痛苦的過程,但是不去校正有可能對最終樹的結(jié)構(gòu)造成影響。
尋求最佳解決辦法中。

23.12.23 加筆
easyminer似乎有一定的改進(jìn),另外又集成了一些別的功能。內(nèi)含子的問題好像有解決,但是效果有待驗證,其次對步驟似乎也有精簡,省去了組裝文件與結(jié)果文件的重復(fù),所以最后結(jié)果中的組裝輸出文件就是最終的結(jié)果。但是過濾文件實在是太大了,如果在結(jié)束后不能自動刪除的話,對電腦容量提出了挑戰(zhàn)。
此外還發(fā)現(xiàn)GeneMiner的運行占用內(nèi)存會隨著基因數(shù)量增加而比較平緩的線性增加,但隨著參考數(shù)量的增加內(nèi)存占會用劇增,因而選擇少量、近緣、最佳的序列作為參考是最優(yōu)的。

24.1.14 加筆
GeneMiner新版本上線(就是上面提到的easyminer),目前僅支持window版本。新版本集合了很多功能,包括OrthoFinder的同源基因篩選、系統(tǒng)發(fā)育重建(建樹)、時間標(biāo)定等,可以對獲取的序列快速篩選(包括根據(jù)與參考序列比對后,序列差異較大的可以直接設(shè)定閾值篩掉)。新版本由于是窗口界面,使用的門檻可以說是非常低了,再加上加入了中英版本一鍵切換功能,完全是為國人設(shè)計的軟件,這一點十分贊。
優(yōu)點在于window,明顯的一個缺點也在window。
首先是算力不夠的問題:一般人的計算機算力都是有限的,無論如何都不可能跟單個服務(wù)器甚至集成去比,隨著參考和目標(biāo)序列的增加而增加的內(nèi)存占用量可不能小看,跑著一個GeneMiner,如果用上的線程過多的話,連Excel都打不開。
再是本地空間不足的問題:一般植物基因組大小有0.5G很正常吧——我想要高質(zhì)量的數(shù)據(jù),至少得測個20X,那么我測個20G數(shù)據(jù)不過分吧(老標(biāo)本破碎的嚴(yán)重我測個50G保個底也不過分吧)——測序公司好心給我都測了些數(shù)據(jù),最后給了我70G雙端數(shù)據(jù)的cleandata還不多收我錢,我不能說這個公司做得不對吧。GeneMiner推薦20G以上內(nèi)存,我D盤還剩27G,結(jié)果跑了一晚上D盤紅了,寫入不了后GeneMiner跑斷了,斷了還只能重跑。如果我的電腦只有256G內(nèi)存,這個批量操作的功能給我也只能喂狗。(這就去裝個2T的好吧,裝照片也不裝測序數(shù)據(jù))
還有個缺點前面的加筆中也反復(fù)提及了,GeneMiner在設(shè)計的時候似乎就沒有考慮基因內(nèi)含子這回事
我又稍微仔細(xì)看了一下工作原理,雖然我沒太看懂,但我大概估摸這個軟件對有內(nèi)含子的基因就是沒法,從原理上來看,只要內(nèi)含子長度過長,就不可能把完整的基因拼出來。GeneMiner是先過濾后拼接的,內(nèi)含子長度長過reads后就有可能有那么一段跟哪個參考都匹配不上(我一個內(nèi)含子長度超過150bp很難嗎),被過濾掉后就沒了。拼接摒棄了SPAdes等成熟的拼接軟件,創(chuàng)新了拼接方法,會邊跟參考比對邊拼接,也就是說基本上完全忠于參考的,那么參考不含內(nèi)含子,那么拼出來的序列就不可能會包含,完全異于Hybpiper的先拼整段再去除內(nèi)含子的思路。這樣造成的結(jié)果就是,有的樣拼出來是只包含這邊一段的外顯子,有的只包含另一邊的外顯子,剩下的既包含部分內(nèi)含子又有部分外顯子,MAFFT后成分復(fù)雜慘不忍睹。如果目標(biāo)基因多,篩掉一些還剩一些,走量的話總是有能用的,但是最后得到的序列是個啥成分很難去評判,只能說位點夠多,足夠保證節(jié)點支持率是100/1.00。