開頭引入可以先從“生信工具的更迭速度其實不慢但真正能沉淀下來變成行業(yè)基礎(chǔ)設(shè)施的很少”切入引出EMBOSS這個老牌工具包的不可替代性。前100字內(nèi)自然融入核心關(guān)鍵詞“生信分析”、“Linux”、“EMBOSS”、“安裝”。寫生信腳本的人多少都有這樣的經(jīng)歷師兄的腳本里赫然寫著transeq、getorf、restrict你往PATH里一查發(fā)現(xiàn)根本沒有這個命令然后開始翻文檔、裝環(huán)境一折騰就是半天。EMBOSSEuropean Molecular Biology Open Software Suite就是這樣一個你遲早會撞上的名字——它不算新但幾乎所有搞序列分析和基因組注釋的生信人都繞不開它。這篇文章就圍繞EMBOSS在Linux環(huán)境下的安裝和核心功能展開既講清楚為什么要用它也把從編譯安裝到常用命令的實操細(xì)節(jié)完整過一遍幫你把這個老牌工具包真正用起來。1. 一個能活二十年的生信工具包EMBOSS到底解決了什么問題1.1 為什么如今還要學(xué)“老古董”EMBOSS最早可以追溯到上世紀(jì)九十年代末歐洲幾個生物信息學(xué)實驗室為了統(tǒng)一序列分析程序而發(fā)起的項目。二十多年過去它的界面談不上炫酷命令行的交互方式也不算友好但在生信分析里它依然占據(jù)一個特殊位置當(dāng)你需要快速完成格式轉(zhuǎn)換、序列比對、引物設(shè)計、限制性酶切位點分析、開放閱讀框預(yù)測這些基礎(chǔ)但高頻的序列操作時EMBOSS往往是最順手的那個工具箱。有人可能會問現(xiàn)在Python生態(tài)里有BiopythonR語言里有Bioconductor為什么還要折騰一個老掉牙的命令行工具我的看法是很多高頻小任務(wù)用命令行工具反而是最高效的。比如你想快速知道一個FASTA文件里有多少條序列、每條序列多長、堿基組成如何跑一條infoseq就完事了用Python寫腳本反而要好幾行代碼。更關(guān)鍵的是EMBOSS的命令行接口設(shè)計得相當(dāng)統(tǒng)一熟悉了一兩個命令之后剩下的都能觸類旁通。1.2 EMBOSS能干的六類核心事情為了讓你對這套工具包有一個整體認(rèn)知我先把EMBOSS最常用的功能類別列出來。這套工具包含五百多個程序覆蓋面很廣但日常使用頻率最高的主要集中在以下幾類功能類別代表性命令典型用途序列檢索與格式轉(zhuǎn)換seqret、infoseqFASTA轉(zhuǎn)GenBank、讀取任一格式序列序列比對needle、water、stretcher全局比對、局部比對、大規(guī)?？焖俦葘Ψg與ORF預(yù)測transeq、getorf、sixpack六框翻譯、開放閱讀框提取限制性酶切分析restrict、remap、recutter查找酶切位點、繪制酶切圖譜引物設(shè)計eprimer3、primersearchPCR引物設(shè)計、引物結(jié)合位點搜索序列特征搜索fuzznuc、cpgplot、einverted核酸序列模式查找、CpG島預(yù)測、反向重復(fù)檢測從這張表基本可以總結(jié)出EMBOSS的定位它不是用來做復(fù)雜統(tǒng)計建?；驒C(jī)器學(xué)習(xí)的平臺而是把日常序列操作中那些“重復(fù)勞動”固化成標(biāo)準(zhǔn)命令讓你在終端里一行命令解決問題輸出結(jié)果穩(wěn)定可復(fù)現(xiàn)。1.3 適合誰看、需要什么基礎(chǔ)這篇文章主要面向三類讀者剛開始接觸生信分析、需要在Linux服務(wù)器上搭建基礎(chǔ)工具鏈的初學(xué)者已經(jīng)用其他語言處理序列、但對命令行工具查漏補(bǔ)缺的從業(yè)者以及做流程搭建、需要在腳本里調(diào)用成熟命令行工具的工程師。前置要求不高只需要你有一個Linux環(huán)境物理機(jī)、虛擬機(jī)、云服務(wù)器都可以會基本的終端操作比如切換目錄、編輯文件。如果這些還不熟也不妨礙閱讀我會把每一步命令寫清楚你照著執(zhí)行就行。2. 動手前先理清環(huán)境依賴、目錄與版本選擇2.1 需要準(zhǔn)備的Linux環(huán)境和基礎(chǔ)依賴我第一次在全新的Ubuntu服務(wù)器上編譯EMBOSS時踩了不少坑回頭總結(jié)發(fā)現(xiàn)核心問題就一個依賴沒裝全。EMBOSS的源碼編譯依賴于幾個基礎(chǔ)組件最主要的是C語言編譯器gcc和make工具。大部分Linux系統(tǒng)默認(rèn)不會安裝完整的編譯工具鏈所以需要提前準(zhǔn)備好。以Ubuntu/Debian系為例編譯前建議先把這些包裝上sudo apt update sudo apt install -y build-essential zlib1g-dev libpng-dev libxpm-dev其中build-essential包含了gcc和makezlib1g-dev是處理壓縮文件所需libpng-dev和libxpm-dev是EMBOSS繪圖和圖形界面功能依賴的庫。如果你用的不是Ubuntu而是CentOS/RHEL系對應(yīng)命令是yum install gcc make zlib-devel libpng-devel libXpm-devel。2.2 選哪種版本源碼發(fā)布版、發(fā)行版打包版、conda版的區(qū)別EMBOSS當(dāng)前的穩(wěn)定主線版本是6.6.0雖然發(fā)布有一段時間了但作者團(tuán)隊仍在維護(hù)后續(xù)的小版本更多是修復(fù)bug和完善文檔。在選擇安裝方式時你有三條路徑從官網(wǎng)下載源碼自己編譯、用系統(tǒng)包管理器直接安裝、通過conda在獨(dú)立環(huán)境里安裝。當(dāng)時我圖省事直接用apt install emboss裝過一版結(jié)果發(fā)現(xiàn)版本偏舊部分新命令缺失而且系統(tǒng)已經(jīng)有一個老版本的Perl庫和它存在兼容性隱患。從那以后我就養(yǎng)成了習(xí)慣在生信項目中只要是能通過conda解決的工具一律優(yōu)先用conda把版本和依賴都鎖在環(huán)境里。但如果追求極致可控比如你想自定義安裝目錄、添加自己寫的模塊源碼編譯還是最穩(wěn)的路徑。2.3 一個容易踩的坑EMBOSS_ACD_ROOT與數(shù)據(jù)目錄EMBOSS的設(shè)計里有一個比較獨(dú)特的概念——ACDAjax Command Definition文件。每次運(yùn)行EMBOSS命令時它都需要讀取對應(yīng)的ACD定義文件來知道參數(shù)類型和默認(rèn)值。源碼編譯安裝完成后ACD文件默認(rèn)放在安裝目錄的share/EMBOSS/acd下數(shù)據(jù)文件則放在share/EMBOSS/data下。問題就出現(xiàn)在這里如果你用非標(biāo)準(zhǔn)路徑安裝或者系統(tǒng)里有多個版本的EMBOSS運(yùn)行時經(jīng)常報錯說找不到.acd文件。解決方式是在環(huán)境變量里顯式指定這兩個路徑。后面我會給出完整配置先記住這個關(guān)鍵點就行因為類似“命令找不到ACD文件”的問題折騰過EMBOSS的人幾乎都遇到過。3. 三種安裝路徑橫向?qū)Ρ仍创a編譯、apt、conda3.1 源碼編譯安裝可控但別忽略configure參數(shù)源碼編譯是最正統(tǒng)的安裝方式也是理解EMBOSS內(nèi)部結(jié)構(gòu)最好的途徑。建議從官網(wǎng)的FTP下載穩(wěn)定版壓縮包wget ftp://emboss.open-bio.org/pub/EMBOSS/EMBOSS-6.6.0.tar.gz tar -zxvf EMBOSS-6.6.0.tar.gz cd EMBOSS-6.6.0接下來是配置、編譯、安裝三步./configure --prefix/opt/emboss make -j4 make install這里強(qiáng)烈建議你在configure時指定--prefix也就是安裝目錄。默認(rèn)情況下EMBOSS會裝到/usr/local下雖然也能用但后續(xù)升級和卸載都不太方便我習(xí)慣把生物信息學(xué)相關(guān)的軟件統(tǒng)一裝在/opt或?qū)ｉT的數(shù)據(jù)盤下比如/opt/bioinfo/emboss這樣管理起來一目了然。make -j4中的-j4表示用4個線程并行編譯能明顯加快速度。如果你的服務(wù)器核心數(shù)多可以改成-j8甚至更高。編譯過程中如果報錯絕大多數(shù)情況都是缺依賴回到2.1節(jié)把依賴補(bǔ)上再重新執(zhí)行make即可。安裝完成后還需要把命令路徑和ACD路徑寫入環(huán)境變量echo export PATH/opt/emboss/bin:$PATH ~/.bashrc echo export EMBOSS_ACD_ROOT/opt/emboss/share/EMBOSS/acd ~/.bashrc echo export EMBOSS_DATA/opt/emboss/share/EMBOSS/data ~/.bashrc source ~/.bashrc3.2 發(fā)行版?zhèn)}庫安裝最快但不一定最省心用包管理器安裝是真的快sudo apt update sudo apt install emboss裝完直接在終端輸入embossversion就能看到版本號。不過我有一次在一些相對精簡的鏡像源上安裝時發(fā)現(xiàn)命令被拆分到了多個包中某些命令缺失需要手動補(bǔ)裝emboss-lib等相關(guān)包。而且發(fā)行版的版本更新往往滯后缺少新特性。這種方式適合臨時應(yīng)急不適合正經(jīng)的生信分析項目。3.3 conda環(huán)境安裝生信工作流里最穩(wěn)妥的選擇如果你已經(jīng)裝了Miniconda或Anaconda強(qiáng)烈推薦用conda來管理EMBOSSconda create -n emboss -c bioconda emboss conda activate embossconda安裝的好處是會自動解決所有依賴并且把EMBOSS放在一個獨(dú)立環(huán)境里不影響系統(tǒng)的全局環(huán)境。我在實際項目中通常會把EMBOSS和其他序列分析工具、Python包放在同一個conda環(huán)境里這樣整個分析流程的環(huán)境是可復(fù)現(xiàn)的。不同項目之間需要不同版本的工具時直接創(chuàng)建多個conda環(huán)境即可。另外一個細(xì)節(jié)conda安裝完成后命令路徑和ACD路徑都被conda環(huán)境管理好了默認(rèn)就能運(yùn)行不用手動設(shè)置環(huán)境變量這對新手非常友好。3.4 安裝完成后的自檢清單不管用哪種方式裝完以后建議從以下三個方面做一次自檢確認(rèn)安裝真的成功了。我在調(diào)試環(huán)境時通常不會只跑一個embossversion而是會實際運(yùn)行幾個命令驗證庫依賴是否完整。示例如下embossversion infoseq -auto -filter test.fasta geecee test.fasta -outfile gc.txt第一個命令確認(rèn)主程序能啟動第二個命令會讀取一個FASTA文件并輸出每條序列的基本信息如果運(yùn)行不報錯說明輸入輸出模塊工作正常第三個命令計算GC含量能進(jìn)一步驗證數(shù)據(jù)處理和結(jié)果寫入沒有問題。跑完這幾個命令后如果沒有異常基本可以判定安裝是完整的可以進(jìn)行下一步實戰(zhàn)了。4. 核心命令逐個拆解從序列讀取到比對翻譯4.1 seqret與infoseq序列格式轉(zhuǎn)換和信息預(yù)覽生信分析中常年伴隨的一個小麻煩就是序列文件格式不統(tǒng)一。你從NCBI下載的是一個GenBank格式的完整記錄從測序公司拿回來的可能是ABI格式的峰圖數(shù)據(jù)而后續(xù)分析軟件又只認(rèn)FASTA格式。seqret就是來解決這個問題的它可以把幾乎任意常見的序列格式轉(zhuǎn)換成你需要的格式。比如把GenBank格式的文件轉(zhuǎn)成FASTAseqret input.gb -outseq output.fastaseqret還支持直接讀取EMBL或GenBank的Accession號通過互聯(lián)網(wǎng)檢索序列并轉(zhuǎn)換比如seqret embl:AB001abc -outseq myseq.fasta當(dāng)然這個功能依賴網(wǎng)絡(luò)環(huán)境如果服務(wù)器不能聯(lián)網(wǎng)直接用本地文件就好。infoseq則像是序列文件的“體檢報告”infoseq myseq.fasta它會輸出一個表格包含每條序列的名稱、長度、GC含量、是否含有未知堿基等信息。你可以在拿到一個新序列文件時先用infoseq快速確認(rèn)文件有沒有問題比如發(fā)現(xiàn)某個樣本的序列長度異常短可能就意味著測序或組裝出了問題。4.2 needle、water、stretcher兩套經(jīng)典算法解決不同比對需求EMBOSS最被人稱道的功能之一是序列比對這里有必要把三個命令的算法背景講清楚因為很多人在選擇時容易弄混。needle是基于Needleman-Wunsch算法的全局比對核心思想是把兩條序列從頭到尾進(jìn)行比對允許引入gap來盡可能匹配所有位置。它適合比對長度相近、結(jié)構(gòu)類似的序列比如兩條來自不同物種的同源蛋白序列needle seq1.fasta seq2.fasta -gapopen 10 -gapextend 0.5 -outfile alignment.txtwater則是基于Smith-Waterman算法的局部比對。它不會強(qiáng)迫整條序列參與比對而是找到兩條序列中相似度最高的區(qū)域適合分析結(jié)構(gòu)域的保守片段、短序列的匹配情況。需要說明的是這里“water”這個名字是EMBOSS作者對自己實現(xiàn)該算法程序的命名與一些文獻(xiàn)中提到的水溶性蛋白無關(guān)初學(xué)者不用把它和結(jié)構(gòu)生物學(xué)的概念聯(lián)想到一起。使用方式water seq1.fasta seq2.fasta -gapopen 10 -gapextend 0.5 -outfile local_aln.txtstretcher是全局比對的“極速版”使用更快的啟發(fā)式算法適合對大片段序列先做粗篩。比如全基因組范圍內(nèi)快速比對、判斷兩條序列有沒有同源關(guān)系用stretcher跑一遍就能得到結(jié)論。實測下來同樣的兩條全基因組序列needle可能要跑幾分鐘stretcher幾秒鐘就出結(jié)果了。在比對命令的參數(shù)中g(shù)apopen和gapextend也非常重要gapopen是引入一個gap的罰分gapextend是每延長一個gap堿基的額外罰分。具體數(shù)值需要根據(jù)序列的生物學(xué)特性來調(diào)整。蛋白序列比對的gap罰分通常低于核酸序列因為密碼子的存在使得核酸比對中插入缺失事件的代價更高所以如果比對核酸序列可以把gapopen設(shè)得稍高一些比如12到15。4.3 transeq、getorf、sixpack翻譯與ORF識別核酸序列拿到手之后最常見的需求就是翻譯成蛋白質(zhì)序列或者找出開放閱讀框ORF。這也是EMBOSS被我使用頻率最高的場景之一。transeq做六框翻譯輸入一條DNA序列輸出三個正鏈和三個負(fù)鏈共六種可能的蛋白序列transeq dna.fasta -outseq protein.fasta -table 1-table 1是遺傳密碼表的編號表示標(biāo)準(zhǔn)密碼子表。如果你做的是線粒體基因組密碼子表就需要換成對應(yīng)物種的編號比如脊椎動物線粒體是-table 2無脊椎動物線粒體是-table 5。密碼子表編號搞錯的后果很嚴(yán)重翻譯出來全是錯位蛋白這個細(xì)節(jié)要特別注意。getorf的作用是從一段核酸序列中提取所有可能的開放閱讀框getorf dna.fasta -outseq orfs.fasta -minsize 300 -find 1-minsize 300表示只保留長度大于等于300個核苷酸的ORF對應(yīng)100個氨基酸-find 1的含義是只輸出正鏈上的ORF如果你想同時找正負(fù)鏈把它改成-find 3即可。我在實際分析中通常先跑一次getorf再用transeq把提取出來的ORF翻譯成蛋白序列然后對蛋白序列做數(shù)據(jù)庫比對判斷這些ORF的真實性。sixpack則把翻譯和ORF標(biāo)注合到了一個程序里輸出一個帶注釋的可視化文本文件適合快速瀏覽序列的編碼潛力。第一次使用sixpack時我發(fā)現(xiàn)它默認(rèn)輸出的文件名是六個字母的組合這個命名規(guī)則確實容易讓人摸不著頭腦。建議在命令里用-outfile顯式指定輸出文件名避免生成一堆難以辨識的文件。4.4 restrict、remap、recutter限制性酶切分析做分子克隆的人對限制性內(nèi)切酶應(yīng)該再熟悉不過。restrict命令內(nèi)置了REBASE限制酶數(shù)據(jù)庫可以在序列上查找所有限制性內(nèi)切酶的識別位點并給出切割后的片段長度restrict plasmid.fasta -enzymes EcoRI -outfile cuts.txt-enzymes參數(shù)是大小寫敏感的指定酶名時必須用標(biāo)準(zhǔn)的命名方式比如 EcoRI、BamHI。如果不指定-enzymesrestrict會調(diào)用全部酶庫輸出結(jié)果可能會非常冗長所以日常使用中建議把它加進(jìn)去只分析關(guān)心的酶。remap和restrict的區(qū)別在于它會生成一個序列圖譜文件把每個酶切位點標(biāo)注在序列的對應(yīng)位置上。你可以用文本編輯器打開這個圖譜文件借助它的位置信息判斷某個酶切位點附近有沒有其他特征序列。recutter是另一個比較有意思的小工具它允許你模擬“改變一個堿基”之后酶切位點的變化情況比如想知道一個SNP突變是否會導(dǎo)致某個酶切位點消失recutter就能直接給出前后對比。這個功能在分子標(biāo)記開發(fā)里非常實用。4.5 eprimer3與primersearch引物設(shè)計的實際用法引物設(shè)計是EMBOSS中最常用的功能之一。eprimer3包裝了經(jīng)典的Primer3算法但你不需要直接去安裝Primer3EMBOSS會調(diào)用內(nèi)置的邏輯完成計算。它需要輸入一段模板序列然后搜索合適的引物對eprimer3 template.fasta -outfile primers.txt輸出文件里會列出候選引物對包括每對引物的序列、長度、Tm值、GC含量以及它們可能形成發(fā)夾結(jié)構(gòu)或引物二聚體的風(fēng)險提示。不過要提醒的是eprimer3的輸出結(jié)果比較質(zhì)樸它會把所有候選信息堆在一個文本文件里。如果你追求圖形界面建議把篩選出的候選區(qū)段轉(zhuǎn)移到Primer-BLAST或者其他在線工具進(jìn)行二次驗證。我自己的使用習(xí)慣是先用eprimer3批量預(yù)篩再用在線工具復(fù)核關(guān)鍵引物對。primersearch則是另一個方向的工具用于在一段長序列中搜索已知引物的結(jié)合位點primersearch -seqall query.fasta -primer primers.txt -mismatchpercent 20-mismatchpercent允許的錯配率需要根據(jù)實驗?zāi)康膩矶?。做普通PCR驗證時20%的容錯率通常夠用但如果引物需要高特異性地結(jié)合靶標(biāo)序列比如做qPCR實驗我會把這個值調(diào)到10%以內(nèi)因為qPCR對引物特異性的要求遠(yuǎn)高于普通PCR。4.6 fuzznuc、cpgplot、einverted特征搜索與基因組結(jié)構(gòu)分析除了上面這些高頻功能EMBOSS還有一批用于序列特征搜索的小工具同樣很值得掌握。fuzznuc用于在一個核酸序列中查找特定的序列模式支持IUPAC簡并密碼子符號fuzznuc genome.fasta -pattern TATAA[AT]A -outfile motif_hits.txt這條命令會搜索類似TATA盒的啟動子基序并且允許最后一位是A或T。fuzzpro是它的蛋白版可以用來搜索蛋白序列中的功能位點。cpgplot則是預(yù)測CpG島的工具輸入一個基因組區(qū)段它就知道這些區(qū)域有沒有CpG島以及落在哪里。啟動子區(qū)甲基化研究、基因組注釋中經(jīng)常用到cpgplot genome_seq.fasta -outfile cpg.txt -outfeat cpg.gffeinverted用于查找序列中的反向重復(fù)結(jié)構(gòu)這類結(jié)構(gòu)在基因調(diào)控、轉(zhuǎn)座子研究中意義重大。核酸序列中方向相反的重復(fù)片段能夠形成發(fā)夾結(jié)構(gòu)影響基因的表達(dá)調(diào)控很多關(guān)鍵調(diào)控元件的注釋都離不開這一功能einverted genome_seq.fasta -outfile inverted.txt -gap 20 -match 3-match和-gap參數(shù)分別控制匹配得分和允許的最大gap大小數(shù)值設(shè)定不同會直接影響重復(fù)序列的檢出靈敏度建議根據(jù)序列復(fù)雜度適當(dāng)調(diào)整。用一個生活類比來收束本節(jié)如果把EMBOSS比作一個綜合工具箱seqret就是改錐套件把“零件”換成你想要的型號needle和water是一把精細(xì)測量的游標(biāo)卡尺測量相似度getorf是把毛坯原料切割成標(biāo)準(zhǔn)件截取ORFeprimer3有點像帶精度檢驗的夾具固定擴(kuò)增位置。工具本身不神秘關(guān)鍵是知道什么場景該拿哪一把。5. 實戰(zhàn)串講一條未知序列從檢查到產(chǎn)出引物的完整流程5.1 場景設(shè)定與輸入數(shù)據(jù)為了讓你真正看到這些命令是怎么串聯(lián)起來的我這里設(shè)定一個典型場景手頭有一段從公共數(shù)據(jù)庫中獲取的未知功能基因序列存在gene_unknown.fasta文件中。任務(wù)目標(biāo)是通過EMBOSS完成四個環(huán)節(jié)的分析先確認(rèn)序列本身沒問題再預(yù)測它是否編碼蛋白、找到ORF區(qū)域然后排查這段序列上有沒有適合做克隆實驗的限制性酶切位點最后設(shè)計一組PCR引物用于擴(kuò)增驗證。5.2 從格式清洗到基本信息統(tǒng)計拿到序列第一件事不是急著分析而是先確認(rèn)文件格式和序列質(zhì)量。我遇到過太多次這樣的情況表面上是FASTA文件其實里面混入了Excel導(dǎo)出的換行符、空格、甚至Tab制表符導(dǎo)致很多下游工具解析失敗。先用infoseq看一下文件概況infoseq gene_unknown.fasta正常情況下輸出會有一個表格列出序列ID、長度、GC含量等信息。如果這一步報錯或者輸出異常就用seqret清洗一遍seqret gene_unknown.fasta -outseq gene_clean.fasta -filter -osformat fasta-filter參數(shù)讓結(jié)果直接輸出到屏幕而不是寫入文件如果不指定它就生成一個默認(rèn)文件這個參數(shù)也可以配合重定向使用。清洗后重新用infoseq查看確認(rèn)序列長度是3的倍數(shù)這暗示它可能是一條完整的編碼區(qū)。同時計算GC含量geecee gene_clean.fastaGC含量異常比如超過65%或低于35%就要警惕了可能意味著樣品污染、測序錯誤或序列拼接有問題后續(xù)分析前需要確認(rèn)。5.3 ORF預(yù)測與蛋白翻譯驗證確認(rèn)序列沒有明顯問題后用getorf找開放閱讀框getorf gene_clean.fasta -outseq gene_orfs.fasta -minsize 300 -find 3-find 3表示同時搜索正鏈和負(fù)鏈上的ORF這樣不會遺漏潛在反向編碼的基因。查看輸出文件里的ORF數(shù)量和長度分布如果只有一個較長的ORF位于序列起始位置附近那么這段序列很可能就是一個完整的CDS。用transeq翻譯這條ORFtranseq gene_orfs.fasta -outseq gene_proteins.fasta -table 1說句實在話這一步對初學(xué)者特別友好——你不需要懂密碼子表命令跑完直接得到蛋白序列肉眼就能看序列里有沒有提前終止的星號標(biāo)注。檢查這段蛋白序列的N端有沒有起始氨基酸MethionineC端是否自然終止如果兩端都符合預(yù)期就說明ORF預(yù)測比較可靠。5.4 酶切位點排查與引物設(shè)計很多克隆實驗設(shè)計的第一步就是確認(rèn)目標(biāo)序列上有沒有關(guān)鍵酶切位點。比如你想用EcoRI和BamHI這兩個酶做克隆就先用restrict檢查它們會不會切斷目標(biāo)序列restrict gene_clean.fasta -enzymes EcoRI,BamHI -outfile gene_cuts.txt查看輸出文件時如果發(fā)現(xiàn)EcoRI或BamHI在序列中出現(xiàn)了酶切位點就說明這兩個酶不適合這個克隆構(gòu)建需要換酶或者引入突變。模擬一個堿基突變來消除酶切位點這種事可以交給recutter去評估。引物設(shè)計用eprimer3eprimer3 gene_clean.fasta -outfile gene_primers.txt打開輸出文件重點看排名靠前的引物對確認(rèn)它們的Tm值落在50到65攝氏度之間、GC含量在40%到60%區(qū)間并且盡量避免3末端出現(xiàn)連續(xù)GC序列。所選引物的擴(kuò)增產(chǎn)物長度也要留意如果只是為了做表達(dá)驗證產(chǎn)物在100到300bp之間比較合適如果需要克隆到載體做后續(xù)實驗產(chǎn)物長度就要覆蓋完整編碼區(qū)。5.5 把常用步驟固化成腳本的想法上面這一套流程如果每次都一條條敲命令效率很低。我的個人習(xí)慣是把這類分析流程寫成shell腳本方便后續(xù)重復(fù)使用。腳本內(nèi)容就是把這些命令按順序排列并加好日志輸出#!/bin/bash SEQ_FILE$1 BASE_NAME$(basename $SEQ_FILE .fasta) infoseq $SEQ_FILE ${BASE_NAME}_info.txt getorf $SEQ_FILE -outseq ${BASE_NAME}_orfs.fasta -minsize 300 -find 3 transeq ${BASE_NAME}_orfs.fasta -outseq ${BASE_NAME}_proteins.fasta -table 1 restrict $SEQ_FILE -enzymes EcoRI,BamHI -outfile ${BASE_NAME}_cuts.txt eprimer3 $SEQ_FILE -outfile ${BASE_NAME}_primers.txt腳本寫完后加上執(zhí)行權(quán)限chmod x run_analysis.sh ./run_analysis.sh gene_clean.fasta這樣每拿到一個新序列只需改文件名就能跑完整套初步分析。人機(jī)交互時間被壓縮到最低出錯率也顯著下降。這一步在效率上帶來的提升遠(yuǎn)比想象中大。6. 安裝和運(yùn)行中常見的報錯與排查經(jīng)驗6.1 編譯階段缺編譯器與缺少開發(fā)頭文件源碼編譯階段最常見的報錯是configure: error: no acceptable C compiler found in $PATH這個錯誤其實不是EMBOSS本身的問題而是系統(tǒng)里壓根沒有g(shù)cc。解決方法是安裝完整的編譯工具鏈sudo apt install -y build-essential安裝后重新運(yùn)行./configure。有時候還會遇到configure: error: libpng not found configure: error: X11/Xlib.h: No such file or directory這就對應(yīng)前面提到的libpng-dev和libxpm-dev沒裝全。這些開發(fā)包的命名習(xí)慣是在庫名后面加-devUbuntu或-develCentOS編譯任何C/C軟件遇到“找不到某個頭文件”的報錯基本都是因為對應(yīng)的-dev包沒裝。順便說一句技巧在報錯提示中搜索No such file or directory看它缺的到底是.h頭文件還是.so動態(tài)庫再去搜索包名比盲目亂試要高效得多。6.2 運(yùn)行階段找不到共享庫與ACD目錄編譯安裝完第一次運(yùn)行命令可能遇到emboss: error while loading shared libraries: libpng.so: cannot open shared object file根本原因是動態(tài)庫的路徑?jīng)]被系統(tǒng)搜索到。源碼編譯安裝的軟件默認(rèn)裝到/usr/local/lib但部分Linux發(fā)行版的新版本系統(tǒng)默認(rèn)不掃描這個目錄。解決辦法是把庫路徑加入ldconfig配置sudo sh -c echo /usr/local/lib /etc/ld.so.conf.d/emboss.conf sudo ldconfig還有一個高頻報錯Unable to open acd file: infoseq.acd這個跟ACD文件路徑配置有關(guān)。檢查你的EMBOSS_ACD_ROOT環(huán)境變量是否指向了正確的目錄echo $EMBOSS_ACD_ROOT find / -name infoseq.acd 2/dev/null如果EMBOSS_ACD_ROOT為空或者指向錯誤目錄重新設(shè)置環(huán)境變量后再執(zhí)行命令。export EMBOSS_ACD_ROOT/opt/emboss/share/EMBOSS/acd同樣地如果運(yùn)行涉及數(shù)據(jù)庫檢索的命令時報錯找不到數(shù)據(jù)庫還需要確認(rèn)EMBOSS_DATA指向了正確的data目錄。6.3 數(shù)據(jù)階段格式識別錯誤與文件路徑問題EMBOSS對輸入文件的格式相當(dāng)嚴(yán)格。運(yùn)行g(shù)etorf或restrict時報錯Unable to read sequence: bad trace or format not recognised這通常是文件格式不標(biāo)準(zhǔn)或內(nèi)容有誤。解決辦法是先通過seqret把文件統(tǒng)一轉(zhuǎn)換成標(biāo)準(zhǔn)FASTA格式再進(jìn)行分析。還有一個常見問題序列ID中如果包含空格、|等特殊符號EMBOSS有時會解析異常。建議在處理前先清理序列標(biāo)題行只保留簡潔的ID。如果你拿到的文件路徑包含中文或空格也要先處理為純英文路徑否則終端解析可能會出意外。6.4 一個排查思路先看命令的ACD定義最后分享一個我自己調(diào)試EMBOSS命令時的通用思路當(dāng)你不確定某個命令能接受哪些參數(shù)、參數(shù)類型是什么時不用去翻文檔直接查看它的ACD定義文件cat $EMBOSS_ACD_ROOT/restrict.acdACD文件的結(jié)構(gòu)很像一份機(jī)器可讀的迷你文檔里面列出了每個參數(shù)的名稱、類型、默認(rèn)值和說明。我多次靠這個方式在離線條件下解決了參數(shù)格式猜不對的問題。比如發(fā)現(xiàn)某個參數(shù)默認(rèn)值寫的是0.5但自己傳一個字符串進(jìn)去就一直報錯一看ACD才知道那個位置需要數(shù)字。從安裝配置到日常使用EMBOSS給我的整體感覺是學(xué)習(xí)曲線確實有一點陡界面設(shè)計也不符合現(xiàn)代審美但它用二十多年積累下來的功能穩(wěn)定性和覆蓋面在生信分析里依然有不可替代的位置。我自己在寫分析流程時到現(xiàn)在還經(jīng)常打開EMBOSS的命令參考找靈感因為它已經(jīng)把序列分析的邊界摸得很清楚了。你不需要背下來所有命令掌握本文提到的那十幾個高頻工具在實戰(zhàn)中遇到具體需求再按ACD文件查參數(shù)就能應(yīng)付絕大部分日常序列分析任務(wù)了。