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