EMBOSS軟件包序列分析程序應(yīng)用實(shí)例

羅 靜 初

(北京大學(xué) 生命科學(xué)學(xué)院，北京大學(xué)生物信息中心,北京 100871)

1 概述

1.1 簡(jiǎn)介

歐洲分子生物學(xué)開(kāi)放軟件包(European Molecular Biology Open Software Suite, EMBOSS)誕生于上個(gè)世紀(jì)九十年代末，是較早投入使用的大型生物信息學(xué)開(kāi)放軟件包，包括300多個(gè)程序，主要用于核酸和蛋白質(zhì)序列分析[1-3]。EMBOSS是歐洲分子生物學(xué)網(wǎng)絡(luò)組織(European Molecular Biology Network, EMBnet)啟動(dòng)的以歐洲國(guó)家為主的國(guó)際合作項(xiàng)目，主要發(fā)起人和開(kāi)發(fā)者為Peter Rice和Alan Bleasby。EMBOSS是開(kāi)源軟件包，源代碼完全公開(kāi)，任何人均可免費(fèi)獲取、安裝、使用和修改，并可進(jìn)行二次開(kāi)發(fā)，例如開(kāi)發(fā)瀏覽器用戶界面等。EMBOSS官方網(wǎng)站除提供軟件下載外，還提供用戶文檔、使用教程、開(kāi)發(fā)指南和常見(jiàn)問(wèn)題解答等相關(guān)資料。

2001年起，筆者在北京大學(xué)開(kāi)設(shè)“實(shí)用生物信息技術(shù)”研究生課程[4]，曾介紹EMBOSS軟件包中部分常用程序。2020年新冠病毒疫情期間，為北京大學(xué)生物信息學(xué)專(zhuān)業(yè)本科生開(kāi)設(shè)“Linux基礎(chǔ)及其在生物信息領(lǐng)域中的應(yīng)用”線上課程，較為系統(tǒng)地介紹了EMBOSS軟件包主要程序。

為推廣EMBOSS軟件包在生物信息學(xué)研究中的應(yīng)用，本文基于教學(xué)中的一些實(shí)例，介紹EMBOSS軟件包中主要程序的用途、用法，及其運(yùn)行結(jié)果的生物學(xué)意義。文中所舉實(shí)例的蛋白質(zhì)序列源自國(guó)際蛋白質(zhì)數(shù)據(jù)庫(kù)UniProt[5]，核酸序列源自美國(guó)國(guó)家生物技術(shù)信息中心(National Center for Biotechnology Information, NCBI)的核酸序列數(shù)據(jù)庫(kù)GenBank和參考序列數(shù)據(jù)庫(kù)RefSeq。文中對(duì)所舉實(shí)例的生物學(xué)背景作了簡(jiǎn)要說(shuō)明，并給出相關(guān)文獻(xiàn)，具體細(xì)節(jié)可參閱“實(shí)用生物信息技術(shù)”課程(Applied Bioinformatics Course, ABC)教學(xué)網(wǎng)站(https://bigd.big.ac.cn/education/)。獲中國(guó)科學(xué)院北京基因組研究所大數(shù)據(jù)中心支持，本文補(bǔ)充材料已上載到該所國(guó)家生物信息中心網(wǎng)站(https://bigd.big.ac.cn/education/)。文中所涉及的網(wǎng)站請(qǐng)參閱補(bǔ)充材料1，所舉實(shí)例中的輸入數(shù)據(jù)可通過(guò)補(bǔ)充材料2中的鏈接下載，運(yùn)行結(jié)果可通過(guò)補(bǔ)充材料3中的鏈接查看。計(jì)劃將文中主要程序?qū)嵗木幊删W(wǎng)絡(luò)教程，并不斷進(jìn)行擴(kuò)充和更新。

1.2 用戶界面

EMBOSS序列分析程序的主要使用方式有以下三種，讀者可根據(jù)自己的實(shí)際情況，適當(dāng)采用以下一種或幾種方式。

1.2.1 命令行方式

EMBOSS基于UNIX操作系統(tǒng)開(kāi)發(fā)，所有程序均可在Linux系統(tǒng)上用命令行方式執(zhí)行。國(guó)內(nèi)外不少生物信息學(xué)相關(guān)科研機(jī)構(gòu)、高等學(xué)校和公司企業(yè)的Linux服務(wù)器上裝有EMBOSS軟件包。讀者可向服務(wù)器管理員申請(qǐng)Linux系統(tǒng)賬號(hào)，登錄后即可通過(guò)命令行方式運(yùn)行該軟件包中的程序。MacOS系統(tǒng)用戶可下載EMBOSS源代碼并編譯、安裝該軟件包?；赪indows系統(tǒng)版本mEMBOSS可在磁盤(pán)操作系統(tǒng)(Disk Operation System, DOS)環(huán)境下以命令行方式運(yùn)行。

1.2.2 圖形界面

命令行方式操作快速高效，是熟悉Linux系統(tǒng)用戶的首選。不熟悉Linux系統(tǒng)的用戶可下載基于Windows系統(tǒng)的mEMBOSS軟件包，在裝有Java運(yùn)行環(huán)境的Windows系統(tǒng)中啟動(dòng)Jemboss圖形界面，運(yùn)行EMBOSS程序[6]。

1.2.3 瀏覽器界面

除上述兩種用戶界面外，還可通過(guò)瀏覽器訪問(wèn)裝有EMBOSS軟件包的服務(wù)器。這類(lèi)服務(wù)器提供基于EMBOSS軟件包開(kāi)發(fā)的用戶界面，如北京大學(xué)生物信息中心開(kāi)發(fā)的生物信息網(wǎng)上實(shí)驗(yàn)室WebLab，荷蘭瓦格寧根大學(xué)(Wageninen University)開(kāi)發(fā)的EMBOSS Explorer，以及歐洲生物信息學(xué)研究所(European Bioinformatics Institute, EBI)安裝的部分EMBOSS程序(補(bǔ)充材料)。

1.3 運(yùn)行模式

上述三種用戶界面中，命令行方式是最基本的運(yùn)行方式，本文僅介紹這種運(yùn)行方式。具體運(yùn)行時(shí)，按參數(shù)選擇方式的不同，可分為以下三種模式。

1.3.1 交互式

第一種為交互式，運(yùn)行時(shí)只輸入程序名，系統(tǒng)給出提示符后再逐項(xiàng)輸入需要處理的序列文件名、輸出結(jié)果文件名和程序參數(shù)。參數(shù)可以有1個(gè)或多個(gè)，用戶可采用系統(tǒng)提供的默認(rèn)參數(shù)，也可自行輸入。下面，我們以血紅蛋白(Hemoglobin)為例，說(shuō)明交互式程序運(yùn)行模式。

血紅蛋白是重要生物大分子，在生命科學(xué)研究歷史中具有特殊地位。國(guó)際蛋白質(zhì)結(jié)構(gòu)數(shù)據(jù)庫(kù)(Protein Data Bank, PDB)分子月報(bào)(Molecule of the Month)網(wǎng)站科普短文介紹了它的結(jié)構(gòu)功能關(guān)系。UniProt數(shù)據(jù)庫(kù)蛋白質(zhì)分子精選(Protein Spotlight)網(wǎng)站介紹了血紅蛋白研究歷史。血紅蛋白是第一個(gè)被確定生理功能的生物大分子，于1849年獲得純化晶體；也是第一個(gè)準(zhǔn)確測(cè)得分子量的蛋白質(zhì)、第一個(gè)在體外非細(xì)胞體系中人工合成的真核生物大分子，其編碼基因的mRNA最先被分離并測(cè)序。人類(lèi)基因組中有alpha和beta兩大類(lèi)血紅蛋白編碼基因，共編碼9種不同的血紅蛋白[4]。成人血液中的血紅蛋白由兩個(gè)alpha血紅蛋白亞基和兩個(gè)beta血紅蛋白亞基組成。研究表明，小鼠alpha血紅蛋白和人的alpha血紅蛋白是直系同源蛋白，起源于共同祖先，其序列相似性較高。

從蛋白質(zhì)數(shù)據(jù)庫(kù)UniProt中下載FASTA格式的人和小鼠alpha血紅蛋白氨基酸序列，分別保存為HBA_HUMAN.FASTA和HBA_MOUSE.FASTA，用EMBOSS軟件包中的程序needle進(jìn)行序列比對(duì)。

needle

Needleman-Wunsch global alignment of two sequences

Input sequence:HBA_HUMAN.FASTA

Second sequence(s):HBA_MOUSE.FASTA

Gap opening penalty [10.0]:

Gap extension penalty [0.5]:

Output alignment [hba_human.needle]:HUMAN-MOUSE.NEEDLE

上述運(yùn)行過(guò)程說(shuō)明如下。

1)用戶輸入程序名needle，按回車(chē)鍵后開(kāi)始運(yùn)行。

2)系統(tǒng)給出所運(yùn)行程序的簡(jiǎn)單說(shuō)明Needleman-Wunsch global alignment of two sequences，并提示用戶輸入用于比對(duì)的序列文件名Input sequence。

3)輸入第一個(gè)序列文件名HBA_HUMAN.FASTA后系統(tǒng)提示輸入第二個(gè)序列文件名。

4)輸入第二個(gè)序列文件名HBA_MOUSE.FASTA后系統(tǒng)提示輸入起始空位罰分值Gap opening penalty，并給出默認(rèn)值10.0，按回車(chē)鍵Enter則采用默認(rèn)值。

5)系統(tǒng)提示輸入延伸空位罰分值Gap extension penalty，并給出默認(rèn)值0.5，按回車(chē)鍵采用默認(rèn)值。

6)系統(tǒng)提示比對(duì)結(jié)果輸出文件名Output alignment，并將第一個(gè)序列FASTA格式第1行注釋信息中的序列名作為默認(rèn)輸出結(jié)果文件名hba_human.needle(默認(rèn)輸出文件名用小寫(xiě)字母)。也可由用戶輸入指定的文件名，如此處的HUMAN-MOUSE.NEEDLE。

1.3.2 參數(shù)式

上述交互式運(yùn)行模式適用于初學(xué)者，用戶可根據(jù)提示信息確定需要輸入的文件名和參數(shù)，許多情況下可先使用默認(rèn)值，在分析運(yùn)行結(jié)果后再次運(yùn)行時(shí)則可適當(dāng)調(diào)節(jié)參數(shù)，以得到更好的結(jié)果。對(duì)于熟練用戶，則可采用參數(shù)模式運(yùn)行程序，即在命令行中直接給出輸入文件名、輸出文件名和程序參數(shù)。例如，用needle程序?qū)θ?、小鼠和大鼠三種哺乳動(dòng)物alpha血紅蛋白進(jìn)行兩兩比對(duì)，則可分別采用以下命令。

needle HBA_HUMAN.FASTA HBA_MOUSE.FASTA -gapo 10.0 -gape 0.5 -out HUMAN-MOUSE.NEEDLE

needle HBA_HUMAN.FASTA HBA_RAT.FASTA -gapo 10.0 -gape 0.5 -out HUMAN-RAT.NEEDLE

needle HBA_MOUSE.FASTA HBA_RAT.FASTA -gapo 10.0 -gape 0.5 -out MOUSE-RAT.NEEDLE

上述needle程序運(yùn)行過(guò)程說(shuō)明如下。

1)第一次對(duì)人和小鼠alpha血紅蛋白進(jìn)行比對(duì)，采用默認(rèn)起始空位罰分-gapo 10.0和默認(rèn)延伸空位罰分-gape 0.5，運(yùn)行結(jié)果存放到輸出文件HUMAN-MOUSE.NEEDLE中。參數(shù)-gapo是-gapopen的簡(jiǎn)略形式，-gape是-gapextension的簡(jiǎn)略形式。

2)第二次對(duì)人和大鼠alpha血紅蛋白進(jìn)行比對(duì)，空位罰分參數(shù)同上，運(yùn)行結(jié)果存放到輸出文件HUMAN-RAT.NEEDLE中。

3)第三次對(duì)小鼠和大鼠alpha血紅蛋白進(jìn)行比對(duì)，空位罰分參數(shù)同上，運(yùn)行結(jié)果存放到輸出文件MOUSE-RAT.NEEDLE中。

上述命令的運(yùn)行結(jié)果包括比對(duì)分值、相同位點(diǎn)數(shù)及百分比、相同及相似位點(diǎn)數(shù)及百分比(見(jiàn)表1)。

表1 人、小鼠和大鼠alpha血紅蛋白兩兩比對(duì)結(jié)果

采用參數(shù)式運(yùn)行模式時(shí)，輸入文件、輸出文件和程序參數(shù)均在命令行中給出，運(yùn)行過(guò)程中不必逐個(gè)輸入；與交互模式相比，運(yùn)行效率有所提高，特別適合批量處理。

1.3.3 菜單式

采用交互式運(yùn)行時(shí)，除個(gè)別參數(shù)可由用戶輸入外，大部分參數(shù)由系統(tǒng)默認(rèn)給定。若用戶需要改變這些默認(rèn)參數(shù)，則可采用菜單式運(yùn)行。即在輸入程序名時(shí)，同時(shí)輸入選項(xiàng)參數(shù)-options，程序運(yùn)行時(shí)則列出所有可選參數(shù)。這種運(yùn)行方式，對(duì)具有較多選擇參數(shù)的程序十分便利。下面，我們以豌豆開(kāi)花后特異表達(dá)基因(Peapost-floral-specificgene,PPF-1)為例，說(shuō)明如何按菜單式運(yùn)行程序getorf，從mRNA中提取編碼區(qū)(Coding sequence, CDS)序列。

1997年，朱玉賢等利用差異表達(dá)方法，從豌豆天然突變體G2株系中分離到開(kāi)花后特異表達(dá)基因。為探索該基因是否與衰老相關(guān)，對(duì)其進(jìn)行了序列分析，初步推斷該基因的表達(dá)產(chǎn)物為內(nèi)膜蛋白(inner-membrane protein)，定位于葉綠體中，與某些細(xì)菌內(nèi)膜蛋白有相同的親疏水性模式[7]。

上述豌豆開(kāi)花后特異表達(dá)基因序列分析的第一步，需要提取mRNA序列中自起始密碼子ATG到終止密碼子之間的編碼區(qū)序列。該序列已遞交NCBI核酸序列數(shù)據(jù)庫(kù)GenBank(登錄號(hào)Y12618)，并作了初步注釋?zhuān)蛄腥L(zhǎng)1 523個(gè)核苷酸，包括第48-1 373位編碼區(qū)、5’端和3’端非翻譯區(qū)(Untranslated region, UTR)。

從核酸序列數(shù)據(jù)庫(kù)GenBank中下載FASTA格式序列文件Y12618.FASTA，采用菜單式運(yùn)行編碼區(qū)提取程序getorf，步驟如下。

getorf-options

Finds and extracts open reading frames (ORFs)

Input nucleotide sequence(s):Y12618.FASTA

Genetic codes

0 : Standard

1 : Standard (with alternative initiation codons)

…… (選項(xiàng)2-23省略)

Code to use [0]:

Minimum nucleotide size of ORF to report [30]:1000

Maximum nucleotide size of ORF to report [1000000]:

Type of sequence to output

0 : Translation of regions between STOP codons

1 : Translation of regions between START and STOP codons

…… (選項(xiàng)2-6省略)

Type of output [0]: 1

protein output sequence(s) [y12618.orf]:Y12618_AA.FASTA

上述運(yùn)行過(guò)程說(shuō)明如下。

1)輸入程序名并啟用菜單式選項(xiàng)getorf -option。

2)輸入豌豆開(kāi)花后特異表達(dá)基因mRNA序列Y12618.FASTA。

3)系統(tǒng)顯示0-23種可選遺傳密碼表，按回車(chē)鍵選擇默認(rèn)通用密碼表。

4)系統(tǒng)顯示最小讀碼框長(zhǎng)度，默認(rèn)為30，輸入1 000，僅獲取長(zhǎng)度大于1 000 bp的編碼區(qū)序列。

5)系統(tǒng)顯示最大讀碼框長(zhǎng)度，按回車(chē)鍵選擇默認(rèn)值1 000 000。

6)系統(tǒng)顯示輸出序列種類(lèi)，共有7種不同選擇，輸入1則提取編碼區(qū)序列并翻譯成氨基酸。

1.4 三個(gè)幫助程序

EMBOSS軟件包整合了300多個(gè)程序，可通過(guò)三個(gè)幫助程序wossname, tfm和seealso了解某個(gè)程序的用途和用法。

1.4.1 wossname

第一個(gè)程序?yàn)閣ossname，其含義為What’s the name，可通過(guò)輸入關(guān)鍵詞查找特定用途的程序名稱(chēng)。例如，可用以下命令找到所有點(diǎn)陣圖(Dot-plot)序列比對(duì)程序。

wossnamedotplot

SEARCH FOR 'DOTPLOT'

dotmatcher Draw a threshold dotplot of two sequences

dotpath Draw a non-overlapping wordmatch dotplot of two sequences

dottup Displays a wordmatch dotplot of two sequences

polydot Draw dotplots for all-against-all comparison of a sequence set

1.4.2 tfm

第二個(gè)幫助程序?yàn)閠fm，其含義為T(mén)he file manual，可用來(lái)顯示某個(gè)程序的使用方法和可選參數(shù)，內(nèi)容十分詳盡，命令參數(shù)部分列出可供用戶選擇的所有參數(shù)及其數(shù)據(jù)類(lèi)型，并給出默認(rèn)值和可選范圍(見(jiàn)表2)。

表2 EMBOSS軟件包中tfm程序幫助信息Table 2 Help information of the tfm program in EMBOSS

命令tfm的功能與Linux系統(tǒng)man命令的功能類(lèi)似，詳細(xì)列出某程序所有幫助信息。此外，也可以在程序運(yùn)行時(shí)用-help參數(shù)列出該程序簡(jiǎn)略信息，包括EMBOSS軟件包的版本。

1.4.3 seealso

第三個(gè)幫助程序?yàn)閟eealso，即英文See also，其含義為列出EMBOSS軟件包中與某程序相關(guān)的其它程序。例如，以下命令列出與點(diǎn)陣圖程序dotmatcher相關(guān)的其它點(diǎn)陣圖程序。

seealsodotmatcher

Finds programs with similar function to a specified program

SEE ALSO

dotpath Draw a non-overlapping wordmatch dotplot of two sequences

dottup Displays a wordmatch dotplot of two sequences

polydot Draw dotplots for all-against-all comparison of a sequence set

2 序列處理

EMBOSS軟件包整合的程序幾乎涵蓋了序列分析的所有方面。本文按功能分類(lèi)，簡(jiǎn)要介紹部分常用程序，包括序列處理程序、序列比對(duì)程序、核酸序列分析程序和蛋白質(zhì)序列分析程序，對(duì)其中一些具有代表性的程序結(jié)合實(shí)例給出具體操作步驟和運(yùn)行結(jié)果，并對(duì)運(yùn)行結(jié)果的生物學(xué)意義略加說(shuō)明。

序列處理是序列分析的基礎(chǔ)，下面我們分別介紹最為常用的格式轉(zhuǎn)換、序列提取和序列變換三類(lèi)序列處理程序。

2.1 格式轉(zhuǎn)換

核酸和蛋白質(zhì)序列格式有多種，不同格式之間的轉(zhuǎn)換在序列分析中經(jīng)常遇到。EMBOSS程序多以FASTA作為輸入序列格式。從GenBank, EMBL等核酸序列數(shù)據(jù)庫(kù)和UniProt等蛋白質(zhì)序列數(shù)據(jù)庫(kù)下載的原始格式序列條目，可轉(zhuǎn)換成FASTA格式。

EMBOSS包括多個(gè)序列格式轉(zhuǎn)換程序，此處介紹最為常用的seqret。該程序是EMBOSS軟件包開(kāi)發(fā)的第一個(gè)程序，除了格式轉(zhuǎn)換外，還有其它許多功能。

2.1.1 seqret用法實(shí)例

從NCBI核酸序列數(shù)據(jù)庫(kù)下載GenBank格式豌豆開(kāi)花后特異表達(dá)基因mRNA序列(登錄號(hào)Y12618)，以Y12618.GB為文件名保存?？捎胹eqret程序?qū)⑵滢D(zhuǎn)換成FASTA格式。采用參數(shù)式方法，直接在命令行指定輸入文件Y12618.GB和輸出文件Y12618.FASTA。

seqret Y12618.GB Y12618.FASTA

seqret - 格式轉(zhuǎn)換程序

Y12618.GB - GenBank格式豌豆開(kāi)花后特異表達(dá)基因mRNA序列文件

Y12618.FASTA - FASTA格式輸出結(jié)果文件

上述豌豆開(kāi)花后特異表達(dá)基因的表達(dá)產(chǎn)物為內(nèi)膜蛋白，已在UniProt蛋白質(zhì)數(shù)據(jù)庫(kù)Swiss-Prot子庫(kù)中收錄。UniProt數(shù)據(jù)庫(kù)中每個(gè)序列都有特定的序列條目名(Entry name)。下載UniProt/Swiss-Prot格式豌豆內(nèi)膜蛋白序列(序列條目名PPF1_PEA)，保存為PPF1_PEA.SW，可用以下命令轉(zhuǎn)換成FASTA格式序列文件。

seqret PPF1_PEA.SW PPF1_PEA.FASTA

seqret - 格式轉(zhuǎn)換程序

PPF1_PEA.SW - UniProt/Swiss-Prot格式豌豆內(nèi)膜蛋白序列文件

PPF1_PEA.FASTA - FASTA格式輸出結(jié)果文件

2.2 序列提取

EMBOSS軟件包中的序列提取程序，包括以下兩大類(lèi)。第一類(lèi)用于FASTA格式輸入文件。這類(lèi)程序操作比較簡(jiǎn)單，常用的有seqretsplit和extractseq，前者可將一個(gè)FASTA格式多序列文件拆分成多個(gè)FASTA格式單序列文件，后者可根據(jù)用戶指定的區(qū)域，提取序列中的子序列，合并成新序列，或按多個(gè)序列保存。

另一類(lèi)程序則用于GenBank和RefSeq格式的核酸序列或UniProt格式的蛋白質(zhì)序列。這些數(shù)據(jù)庫(kù)中的序列條目通常均包含序列特征注釋信息，也稱(chēng)序列特征表(Feature Table)。這里程序則可根據(jù)序列特征表中提供的注釋信息，提取其中的子序列。核酸序列數(shù)據(jù)庫(kù)的序列特征表包括mRNA、編碼序列、翻譯產(chǎn)物蛋白質(zhì)、外顯子(Exon)、內(nèi)含子(Intron)、非翻譯區(qū)、序列標(biāo)簽位點(diǎn)(Sequence Tag Site, STS)等。蛋白質(zhì)序列數(shù)據(jù)庫(kù)的序列特征表包括二級(jí)結(jié)構(gòu)(HELIX, STRAND, TURN)、跨膜螺旋(TRANSMEM)、變異位點(diǎn)(VARIANT)、活性位點(diǎn)(ACT_SITE)、金屬結(jié)合位點(diǎn)(METAL)、糖基化位點(diǎn)(CARBOHYDR)、DNA結(jié)合區(qū)域(DNA_BIND)、二硫鍵(DISULFID)、信號(hào)肽(SIGNAL)、序列模體(MOTIF)和結(jié)構(gòu)域(DOMAIN)等[5]。

這類(lèi)程序中最為常用的有coderet和extractfeat。coderet用法比較簡(jiǎn)單，可用于提取mRNA、編碼序列和所編碼的蛋白質(zhì)序列。extractfeat功能十分強(qiáng)大，用法也比較靈活，可提取更多種類(lèi)的子序列，包括外顯子、內(nèi)含子、重復(fù)序列和多聚腺苷酸信號(hào)等。此外，通過(guò)設(shè)定特征表類(lèi)型(Type)、標(biāo)簽(Tag)和標(biāo)簽值(Value)等參數(shù)，也可提取用戶指定的一個(gè)或幾個(gè)特定子序列。

2.2.1 coderet用法實(shí)例

以小鼠alpha血紅蛋白編碼基因(GenBank登錄號(hào)V00714)DNA序列為例，根據(jù)序列注釋信息，該基因包括三個(gè)外顯子，其mRNA序列位于372-500, 623-826和961-1 191，編碼區(qū)位于405-500, 623-826和 961-1 089。運(yùn)行coderet程序，可分別提取mRNA序列、編碼區(qū)序列、翻譯產(chǎn)物蛋白質(zhì)序列和非編碼區(qū)序列。

coderet V00714.GB

Extract CDS, mRNA and translations from feature tables

Output file [v00714.coderet]:

Coding nucleotide output sequence(s) (optional) [v00714.cds]:

Messenger RNA nucleotide output sequence(s) (optional) [v00714.mrna]:

Translated coding protein output sequence(s) (optional) [v00714.prot]:

Non-coding nucleotide output sequence(s) (optional) [v00714.noncoding]:

調(diào)用coderet程序，輸入小鼠alpha血紅蛋白GenBank格式文件V00714.GB，程序以GenBank序列條目中基因座(LOCUS)名稱(chēng)v00714為默認(rèn)輸出文件名，將輸出結(jié)果分別保存到5個(gè)文件中，分別為列表文件(v00714.coderet)、編碼區(qū)序列文件(v00714.cds)、mRNA序列文件(v00714.mrna)、所編碼的蛋白質(zhì)序列文件(v00714.prot)和非編碼區(qū)序列文件(v00714.noncoding)。按EMBOSS軟件包習(xí)慣，默認(rèn)輸出文件名用小寫(xiě)字母表示。

2.2.2 extractfeat用法實(shí)例

以上述小鼠alpha血紅蛋白編碼基因?yàn)槔?，根?jù)序列注釋信息，該基因包括兩個(gè)內(nèi)含子，分別位于501-622和827-960，用以下命令可提取這兩個(gè)內(nèi)含子的序列：

extractfeat V00714.GB V00714.INTRON -type"intron"

extractfeat - 子序列提取程序

V00714.GB - GenBank格式小鼠alpha血紅蛋白編碼基因序列

V00714.INTRON - 輸出結(jié)果內(nèi)含子序列

-type "intron" - 指定提取注釋信息中的內(nèi)含子intron

2.3 序列變換

EMBOSS軟件包中的序列變換程序包括reveseq, msbar和 shuffleseq等。程序revseq將已知序列轉(zhuǎn)換成反向互補(bǔ)序列，msbar對(duì)已知序列進(jìn)行突變，shuffleseq則用于產(chǎn)生隨機(jī)序列。

程序msbar可用于對(duì)已知序列進(jìn)行單點(diǎn)或多點(diǎn)隨機(jī)突變，突變方式可以是替換、插入或刪除，突變位點(diǎn)可以是單個(gè)核苷酸點(diǎn)突變(Point mutation)，也可插入或刪除一個(gè)序列片段(Block mutation)，還可插入或刪除一個(gè)密碼子(Codon mutation)。程序運(yùn)行默認(rèn)方式為交互式，即屏幕顯示交互菜單，用戶可以自行選擇突變次數(shù)(Number of times)，也可按上述三種不同突變時(shí)，選擇插入(Insertion)、刪除(Deletion)、替換(Substitution)、復(fù)制(Duplication)和移動(dòng)(Move)等不同突變方式。以豌豆開(kāi)花后特異表達(dá)基因編碼區(qū)序列為例，以下命令對(duì)該序列進(jìn)行1次單點(diǎn)插入突變、1次片段刪除突變和1次密碼子替換突變。

msbar Y12618_CDS.FASTA Y12618_CDS_NEW.FASTA

Mutate a sequence

Number of times to perform the mutation operations [1]: 1

Point mutation operations

0 : None

1 : Any of the following

2 : Insertions

3 : Deletions

4 : Changes

5 : Duplications

6 : Moves

Types of point mutations to perform [0]: 2

Block mutation operations

(6個(gè)選項(xiàng)與點(diǎn)突變相同，略)

Types of block mutations to perform [0]: 3

Codon mutation operations

(6個(gè)選項(xiàng)與點(diǎn)突變相同，略)

Types of codon mutations to perform [0]: 4

熟練用戶也可在命令行中直接設(shè)定參數(shù)，即：

msbar Y12618_CDS.FASTA Y12618_CDS_NEW.FASTA -count 1 -point 2 -block 3 -codon 4

msbar - 序列變換程序

Y12618_CDS.FASTA - FASTA格式豌豆開(kāi)花后特異表達(dá)基因編碼區(qū)序列

Y12618_CDS_NEW.FASTA - 隨機(jī)突變輸出結(jié)果

-count 1 - 突變次數(shù)1次

-point 2 - 單點(diǎn)突變方式選擇插入(2: Insertions)

-block 3 - 片段突變方式選擇刪除(3: Deletions)

-codon 4 - 密碼子突變方式選擇改變(4: Change)

突變結(jié)果可用雙序列比對(duì)程序needle驗(yàn)證。由于程序msbar基于隨機(jī)突變，即使運(yùn)行時(shí)設(shè)定的參數(shù)完全一致，兩次運(yùn)行結(jié)果并不相同。

2.4 序列顯示

EMBOSS軟件包中的序列顯示程序包括infoseq、showseq和showfeat等。程序infoseq顯示序列簡(jiǎn)單信息，包括序列名稱(chēng)、長(zhǎng)度和GC含量等，showseq按不同格式輸出序列，而showfeat則根據(jù)序列特征表輸出序列注釋信息。

2.4.1 infoseq用法實(shí)例

程序infoseq簡(jiǎn)單實(shí)用，可用于顯示序列名稱(chēng)、格式及登錄號(hào)等基本信息，并可統(tǒng)計(jì)序列長(zhǎng)度。對(duì)于核酸序列，還能統(tǒng)計(jì)GC含量。可用以下命令顯示豌豆開(kāi)花后特異表達(dá)基因的基本信息。

infoseq Y12618.FASTA -outfile Y12618.INFO

infoseq - 序列信息顯示程序

Y12618.FASTA - FASTA格式豌豆開(kāi)花后特異表達(dá)基因mRNA序列

Y12618.INFO - 輸出結(jié)果文件

程序infoseq既可用于GenBank和Swiss-Prot等格式，也可用于FASTA格式；既可用于單個(gè)序列，也可用于多序列文件?？捎靡韵旅铒@示12個(gè)人源癌胚抗原(Carcinoembryonic Antigen, CEA)蛋白質(zhì)分子的基本信息[8]。

infoseq 12HUMAN_CEA.FASTA -outfile 12HUMAN_CEA.INFO

infoseq - 序列信息顯示程序

12HUMAN_CEA.FASTA - 12個(gè)FASTA格式人癌胚抗原蛋白質(zhì)序列

12HUMAN_CEA.INFO - 輸出結(jié)果文件

2.4.2 showseq用法實(shí)例

程序showseq可用不同方式顯示核酸序列，也可顯示按不同讀碼框翻譯得到的氨基酸序列、反向互補(bǔ)序列，以及酶切位點(diǎn)等。以下命令顯示豌豆開(kāi)花后特異表達(dá)基因第1-120位序列，并用標(biāo)尺顯示序列位點(diǎn)，第48-120位用大寫(xiě)字母顯示。

showseq Y12618.GB Y12618.SHOWSEQ -sbegin 1 -send 120 -format 3 -upper"48-120"

showseq - 核酸子序列顯示程序

Y12618.GB - GenBank格式豌豆開(kāi)花后特異表達(dá)基因序列

Y12618.SHOWSEQ - 輸出結(jié)果文件

-sbegin 1 - 指定子序列開(kāi)始位點(diǎn)為1

-send 120 - 指定子序列終止位點(diǎn)為120

-format 3 - 指定輸出格式

-upper "48-120" - 指定48-120位用大寫(xiě)字母表示

2.4.3 showfeat用法實(shí)例

程序showfeat可用于顯示GenBank和Swiss-Prot等序列的特征信息。以下命令以圖形方式顯示GenBank格式小鼠alpha血紅蛋白編碼基因V00714.GB中外顯子和內(nèi)含子位置。

showfeat V00714.GB V00714.SHOWFEAT

showfeat - 序列特征信息顯示程序

V00714.GB - GenBank格式小鼠血紅蛋白編碼基因序列

V00714.SHOWFEAT - 輸出結(jié)果文件

3 序列比對(duì)

序列比對(duì)在生物信息學(xué)中占有重要地位，是核酸和蛋白質(zhì)序列分析的基礎(chǔ)。EMBOSS軟件包整合了十多個(gè)序列比對(duì)程序，包括雙序列比對(duì)、多序列比對(duì)、數(shù)據(jù)庫(kù)搜索，以及基于點(diǎn)陣圖的可視化序列比對(duì)等。

3.1 雙序列比對(duì)

EMBOSS中整合的雙序列比對(duì)程序包括needle, water, stretcher, matcher, seqmatcherall, supermatcher和esim4等，其中needle和stretcher為基于全局相似性的序列比對(duì)程序，其余為基于局部相似性的序列比對(duì)程序。needle和water最為常用，廣泛用于核酸和蛋白質(zhì)序列比對(duì)。它們均基于動(dòng)態(tài)規(guī)劃算法，在給定計(jì)分矩陣和空位罰分前提下，能夠得到最佳比對(duì)結(jié)果，即最優(yōu)解。程序needle所采用的算法由Needleman和Wunsch于1970年提出，而程序water所采用的算法由Smith和Waterman于1981年提出。程序stretcher是在needle基礎(chǔ)上稍作修改，運(yùn)行時(shí)所需內(nèi)存大為降低，而運(yùn)行時(shí)間稍長(zhǎng)。而程序matcher則是water的改進(jìn)版，可由用戶指定輸出一個(gè)或多個(gè)最佳局部比對(duì)結(jié)果。程序seqmatcherall和supermatcher用于多條序列比對(duì)或數(shù)據(jù)庫(kù)搜索，運(yùn)行時(shí)間較長(zhǎng)。此外，esim4是將mRNA序列定位于基因組序列的程序，而est2genome則是將表達(dá)序列標(biāo)簽(Expressed Sequence Tag, EST)定位于基因組序列的程序。

3.1.1 needle用法實(shí)例

下面，我們以人癌胚抗原為例，說(shuō)明全局比對(duì)程序needle和局部比對(duì)程序water的用途和用法。

人癌胚抗原是一種細(xì)胞表面糖蛋白，多在直腸癌、胃癌等惡性腫瘤中表達(dá)[8]。CEA基因家族分CEA和妊娠特異性beta-1糖蛋白(Pregnancy-specific beta-1-glycoprotein, PSG)兩個(gè)亞家族，其中CEA亞家族包括12個(gè)不同成員(見(jiàn)表3)。CEA蛋白質(zhì)分子屬免疫球蛋白超家族，N-端含長(zhǎng)度為34個(gè)氨基酸的信號(hào)肽，第35位開(kāi)始則為免疫球蛋白可變結(jié)構(gòu)域(Immunoglobin Variable Domain, IgV)，長(zhǎng)度約為110個(gè)氨基酸。除可變結(jié)構(gòu)域外，有的CEA分子還含一個(gè)或多個(gè)免疫球蛋白恒定結(jié)構(gòu)域(Immunoglobin Constant Domain, IgC)，分不同亞型。

表3 UniProt/Swiss-Prot中收錄的12個(gè)人源癌胚抗原蛋白質(zhì)分子Table 3 Twelve human CEA proteins in UniProt/Swiss-Prot

UniProt/Swiss-Prot蛋白質(zhì)數(shù)據(jù)庫(kù)中收錄了12個(gè)人源癌胚抗原蛋白質(zhì)CEA(見(jiàn)圖1，根據(jù)德國(guó)Ludwig-Maximilians大學(xué)Zimmermann教授CEA網(wǎng)站改編)。圖中標(biāo)有N的結(jié)構(gòu)域?yàn)槊庖咔虻鞍卓勺兘Y(jié)構(gòu)域IgV；標(biāo)有A和B的結(jié)構(gòu)域?yàn)槊庖咔虻鞍缀愣ńY(jié)構(gòu)域IgC，各分3種亞型(A1, A2, A3和B1, B2, B3)。嵌入磷脂雙層膜的箭頭表示糖基磷脂酰肌醇(Glycosylphosphatidylinositol, GPI)膜結(jié)合位點(diǎn)，穿過(guò)磷脂雙層膜的螺旋表示跨膜螺旋(Transmembrane helix, TMH)。

人的III型(CEAM3_HUMAN)和IV型(CEAM4_HUMAN)CEA分子長(zhǎng)度接近，各含1個(gè)可變結(jié)構(gòu)域、1個(gè)跨膜螺旋區(qū)和1個(gè)膜內(nèi)區(qū)。用全局比對(duì)程序needle對(duì)其進(jìn)行序列比對(duì)，命令如下。

needle CEAM3_HUMAN.FASTA CEAM4_HUMAN.FASTA CEAM3-CEAM4.NEEDLE -gapo 20 -gape 2

needle-EMBOSS程序，用于全局序列比對(duì)

CEAM3_HUMAN.FASTA-FASTA格式III型癌胚抗原分子序列

CEAM4_HUMAN.FASTA-FASTA格式IV型癌胚抗原分子序列

CEAM3-CEAM4.NEEDLE-輸出結(jié)果文件

-gapo 20- 起始空位罰分

-gape 2-延伸空位罰分

分析比對(duì)結(jié)果可以發(fā)現(xiàn)，這兩個(gè)亞型的CEA分子具有較高的相似性，僅在C-端有1個(gè)長(zhǎng)度為8個(gè)氨基酸殘基的插入。上述命令中，起始空位罰分設(shè)為20，而不用默認(rèn)值10，延伸空位罰分設(shè)為2，而不用默認(rèn)值0.5，可用來(lái)避免不必要的插入或刪除。

圖1 UniProt/Swiss-Prot數(shù)據(jù)庫(kù)中12個(gè)人源癌胚抗原蛋白質(zhì)分子Fig.1 Twelve human CEA proteins in UniProt/Swiss-Prot

3.1.2 water用法實(shí)例

全局比對(duì)程序needle適用于長(zhǎng)度相差不大的兩個(gè)序列，如上述CEAM3_HUMAN和CEAM4_HUMAN，而CEAM5_HUMAN含7個(gè)結(jié)構(gòu)域，除可變結(jié)構(gòu)域IgV外，還包括6個(gè)恒定結(jié)構(gòu)域IgC。若需比較其可變結(jié)構(gòu)域與CEAM3_HUMAN的可變結(jié)構(gòu)域序列相似性，則可用局部比對(duì)程序water進(jìn)行比對(duì)。這兩個(gè)分子N-端均有長(zhǎng)度為34個(gè)氨基酸的信號(hào)肽，C-端有跨膜螺旋，免疫球蛋白結(jié)構(gòu)域位于35-142區(qū)域，比對(duì)時(shí)可用參數(shù)指定。

water CEAM3_HUMAN.FASTA -sbegin 35 -send 142 CEAM5_HUMAN.FASTA -sbegin 35 -send 142 CEAM3-CEAM5.WATER -gapo 20 -gape 2

water-雙序列局部比對(duì)程序

CEAM3_HUMAN.FASTA-FASTA格式III型癌胚抗原分子序列

CEAM5_HUMAN.FASTA-FASTA格式V型癌胚抗原分子序列

CEAM3-CEAM5.WATER-輸出結(jié)果文件

-sbegin 35-指定比對(duì)序列起始位點(diǎn)35

-send 142-指定比對(duì)序列終止位點(diǎn)142

-gapo 20-起始空位罰分

-gape 2-延伸空位罰分

比對(duì)結(jié)果表明，這兩個(gè)CEA分子N-端可變結(jié)構(gòu)域序列具有較高相似性。

3.2 多序列比對(duì)

和雙序列比對(duì)一樣，多序列比對(duì)在核酸和蛋白質(zhì)序列分析中也廣泛使用。EMBOSS軟件包中整合了多序列比對(duì)程序emma和edialign。程序emma是EMBOSS整合的基于全局相似性多序列比對(duì)程序ClustalW，而edialign則是EMBOSS整合的基于全局加局部相似性的多序列比對(duì)程序Dialign。

3.2.1 emma用法實(shí)例

下面以血紅蛋白為例，說(shuō)明emma用法。UniProt/Swiss-Prot中收錄了9個(gè)已經(jīng)審閱的人源血紅蛋白[4]，分屬于alpha和beta兩個(gè)亞家族。alpha亞家族有5個(gè)基因，編碼4種蛋白質(zhì)，其中HBA1和HBA2兩個(gè)基因編碼的蛋白質(zhì)序列完全相同；beta亞家族也有5個(gè)基因，編碼5種蛋白質(zhì)(見(jiàn)表4)。

表4 UniProt/Swiss-Prot中收錄的人的9種不同血紅蛋白Table 4 Nine human hemoglobins in UniProt/Swiss-Prot

可用以下命令對(duì)上述9個(gè)血紅蛋白進(jìn)行多序列比對(duì)，除輸出FASTA格式序列比對(duì)文件外，同時(shí)輸出Newick格式分支圖文件，可用MEGA軟件顯示其樹(shù)形分支結(jié)構(gòu)。

emma 9HUMAN_HB.FASTA 9HUMAN_HB.ALN 9HUMAN_HB.DND

emma - 多序列比對(duì)程序

9HUMAN_HB.FASTA - 9個(gè)FASTA格式人源血紅蛋白序列

9HUMAN_HB.ALN - FASTA格式輸出結(jié)果文件

9HUMAN_HB.DND - Newick格式輸出結(jié)果文件

程序emma有許多可調(diào)參數(shù)，包括計(jì)分矩陣、空位罰分、比對(duì)方式和輸出格式等，可用菜單運(yùn)行模式，即運(yùn)行時(shí)加-options參數(shù)，即可指定上述參數(shù)的值。

3.2.2 edialign用法實(shí)例

程序emma多用于全局比對(duì)，如上述9個(gè)長(zhǎng)度相差不大的血紅蛋白，而edialign采用全局比對(duì)加局部比對(duì)的方法，適用于尋找蛋白質(zhì)序列中具有局部相似性的保守結(jié)構(gòu)域或核酸序列中保守序列模體(Motif)。例如，12個(gè)人源癌胚抗原可用edialign進(jìn)行多序列比對(duì)。

edialign 12HUMAN_CEA.FASTA 12HUMAN_CEA.EDIA 12HUMAN_CEA.ALN

edialign - 多序列比對(duì)程序

12HUMAN_CEA.FASTA - FASTA格式12個(gè)人源癌胚抗原序列

12HUMAN_CEA.EDIA - edialign格式比對(duì)輸出結(jié)果文件

12HUMAN_CEA.ALN - FASTA格式比對(duì)輸出結(jié)果文件

輸出結(jié)果保存到兩個(gè)文件中，12HUMAN_CEA.EDIA是多序列比對(duì)格式，比對(duì)結(jié)果中保守區(qū)域用大寫(xiě)字母表示，每個(gè)位點(diǎn)標(biāo)有數(shù)字0-9，數(shù)字越大，保守性越高。12HUMAN_CEA.ALN為FASTA格式的比對(duì)結(jié)果。

3.3 點(diǎn)陣圖

點(diǎn)陣圖也是序列比對(duì)中常用方法，其特點(diǎn)是輸出結(jié)果直觀。EMBOSS中整合了4個(gè)點(diǎn)陣圖程序，即dottup, dotpath, dotmatcher和polydot。程序polydot用于多序列比對(duì)，而其余3個(gè)程序用于雙序列比對(duì)。運(yùn)行點(diǎn)陣圖程序時(shí)，通常需要指定滑動(dòng)窗口大小，若滑動(dòng)窗口中兩個(gè)序列片段相似性超過(guò)用戶指定的閾值，則在平面坐標(biāo)系中用點(diǎn)標(biāo)出。需要注意的是，dottup和dotpath只考慮指定大小的滑動(dòng)窗口中兩個(gè)序列片段中相同核苷酸或氨基酸，可用于核酸或蛋白質(zhì)序列比對(duì)；而dotmatcher不僅考慮相同殘基，同時(shí)根據(jù)計(jì)分矩陣考慮氨基酸殘基之間的相似性，只能用于蛋白質(zhì)序列比對(duì)。

3.3.1 dottup用法實(shí)例

下面，以河豚魚(yú)質(zhì)粒片段DNA序列為例，說(shuō)明dottup的用法。人的多藥耐藥(Multidrug Resistance, MDR)基因家族包括MDR1, MDR3等幾種不同亞型，其表達(dá)產(chǎn)物為膜通道糖蛋白，利用ATP提供能量，將藥物等細(xì)胞內(nèi)外源物質(zhì)運(yùn)送到胞外從而產(chǎn)生抗藥性。為探索MDR耐藥機(jī)制，劉勇于1998年從模式生物河豚魚(yú)(Takifugurubripes, Fugu)柯氏質(zhì)粒(Cosmid)中克隆到兩個(gè)人的MDR同源基因(補(bǔ)充材料)。這兩個(gè)河豚魚(yú)MDR基因頭尾相接串聯(lián)排列，測(cè)序拼接得到的全長(zhǎng)序列約為40 kb。該河豚魚(yú)序列片段已提交NCBI核酸序列數(shù)據(jù)庫(kù)GenBank(登錄號(hào)AF164138)。下載FASTA格式的河豚魚(yú)DNA序列片段，利用點(diǎn)陣圖程序dottup，可以輸出圖形文件，顯示這兩個(gè)基因的大體位置。

EMBOSS軟件包中dottup等程序運(yùn)行圖形文件格式可由用戶選擇，缺省為X11，若裝有圖形顯示終端(X-Terminal)，可直接在屏幕上輸出。也可保存為其它格式的圖形文件，如可縮放矢量圖形格式(Scalable Vector Graphics, SVG)和可移植網(wǎng)絡(luò)圖形格式(Portable Network Graphics, PNG)。

以下是利用EMBOSS軟件包中點(diǎn)陣圖程序dottup分別對(duì)河豚魚(yú)基因組序列片段進(jìn)行比對(duì)的命令和所用參數(shù)。

dottup AF164138.FASTA AF164138.FASTA -graph svg -goutfile AF164138 -gtitle'Cosmid'-gsubtitle'AF164138'-word 13

dottup - 繪制點(diǎn)陣圖程序

-graph svg - 輸出結(jié)果圖形格式為SVG

A164138.FASTA - 河豚魚(yú)基因組片段DNA序列

-goutfile AF164138 - 輸出結(jié)果圖形文件名

-gtitle 'Cosmid' - 輸出結(jié)果圖形標(biāo)題

-gsubtitle 'AF164138' - 輸出結(jié)果圖形副標(biāo)題

-word 13 - 滑動(dòng)窗口大小，缺省為10，此處取13，以減少背景噪聲

上述命令輸出結(jié)果顯示兩條與對(duì)角線平行的線段(見(jiàn)圖2a)，表明該基因組序列片段5’端有兩個(gè)長(zhǎng)度約為13kb相似片段，即兩個(gè)串聯(lián)重復(fù)多藥耐藥基因MDR。

用以下命令，設(shè)定比對(duì)范圍，則可進(jìn)一步確定這兩個(gè)串聯(lián)重復(fù)基因的相似性。

dottup AF164138.FASTA -send 13001 AF164138.FASTA -sbegin 13001 -send 26000 -graph svg -goutfile AF164138_GENE -gtitle'Gene'-gsubtitle'AF164138'-word 13

-send 13001 - 指定第1個(gè)序列終止位點(diǎn)為13 000，起始位點(diǎn)默認(rèn)為1

-sbegin 13001 - 指定第2個(gè)序列起始位點(diǎn)為13 001

-send 26000 - 指定第2個(gè)序列終止位點(diǎn)為26 000

從輸出結(jié)果(見(jiàn)圖2b)中可以看出，這兩個(gè)序列片段具有一定相似性，有些區(qū)域相似性較高，圖中為連接在一起的線段，而有些區(qū)域相似性較低，可能是基因中內(nèi)含子區(qū)域。

查看該基因組序列注釋信息，發(fā)現(xiàn)這兩個(gè)基因由20多個(gè)外顯子組成。利用coderet程序，提取編碼序列，運(yùn)行dottup程序，比較這兩個(gè)編碼序列的相似性。

dottup AF164138_CDS_1.FASTA AF164138_CDS_2.FASTA -graph svg -goutfile AF164138_CDS -gtitle'CDS'-gsubtitle'AF164138'-word 8

AF164138_CDS_1.FASTA - 第1個(gè)基因編碼序列

AF164138_CDS_2.FASTA - 第2個(gè)基因編碼序列

-word 8 - 序列比對(duì)時(shí)滑動(dòng)窗口大小，默認(rèn)為10，此處取8，以增加靈敏度

輸出結(jié)果(見(jiàn)圖2c)顯示，編碼區(qū)序列相似性比全長(zhǎng)基因序列相似性更高。運(yùn)行dottup程序，可進(jìn)一步比較所編碼蛋白質(zhì)序列相似性。

dottup AF164138_PRO_1.FASTA AF164138_PRO_2.FASTA -graph svg -goutfile AF164138_PRO -gtitle'Protein'-gsubtitle'AF164138'-word 6

AF164138_PRO_1.FASTA - 第1個(gè)基因編碼蛋白質(zhì)序列

AF164138_PRO_2.FASTA - 第2個(gè)基因編碼蛋白質(zhì)序列

-word 6 - 序列比對(duì)時(shí)滑動(dòng)窗口大小，缺省為10，此處取6，以增加靈敏度

輸出結(jié)果(見(jiàn)圖2d)顯示，所編碼兩個(gè)蛋白質(zhì)序列同樣具有較高相似性。

3.3.2 Dotmatcher用法實(shí)例

點(diǎn)陣圖程序dottup多用于核酸序列，而dotmatcher則可用于蛋白質(zhì)序列。下面以果蠅體節(jié)發(fā)育相關(guān)基因?yàn)槔?，說(shuō)明利用dotmatcher顯示序列中的重復(fù)片段。該基因編碼長(zhǎng)度為1 504個(gè)氨基酸的蛋白質(zhì)(UniProt序列條目名SLIT_DROME)。可用以下dotmatcher命令和參數(shù)，得到不同輸出結(jié)果。

圖2 河豚魚(yú)多藥耐藥基因點(diǎn)陣圖序列比對(duì)結(jié)果Fig.2 Dot-plot alignment output of Fugu multidrug resistance gene

dotmatcher SLIT_DROME.FASTA SLIT_DROME.FASTA -graph svg -goutfile SLIT_DROME_W10T23 -window 10 -threshold 23

dotmatcher SLIT_DROME.FASTA SLIT_DROME.FASTA -graph svg -goutfile SLIT_DROME_W24T20 -window 24 -threshold 20

dotmatcher SLIT_DROME.FASTA SLIT_DROME.FASTA -graph svg -goutfile SLIT_DROME_W38T20 -window 38 -threshold 20

dotmatcher SLIT_DROME.FASTA SLIT_DROME.FASTA -graph svg -goutfile SLIT_DROME_W38T30 -window 38 -threshold 30

dotmatcher - 繪制蛋白質(zhì)序列點(diǎn)陣圖程序

SLIT_DROME.FASTA - 果蠅體節(jié)發(fā)育相關(guān)基因蛋白質(zhì)序列

-window 10 - 滑動(dòng)窗口大小為默認(rèn)值10個(gè)氨基酸殘基

-threshold 23 - 相似性閾值為默認(rèn)值23

-window 24 - 滑動(dòng)窗口設(shè)為24個(gè)氨基酸殘基

-threshold 20 - 相似性閾值設(shè)為20

-window 38 - 滑動(dòng)窗口設(shè)為38個(gè)氨基酸殘基

-threshold 20 - 相似性閾值設(shè)為20

-window 38 - 滑動(dòng)窗口設(shè)為38個(gè)氨基酸殘基

-threshold 30 - 相似性閾值設(shè)為30

查看數(shù)據(jù)庫(kù)中SLIT_DROME序列特征表注釋信息，其N(xiāo)末端有4個(gè)區(qū)域富含亮氨酸重復(fù)片段(Leucine Rich Repeat, LRR)，每個(gè)區(qū)域由6個(gè)重復(fù)片段組成，每個(gè)重復(fù)片段約含24個(gè)氨基酸殘基，序列上有一定保守性；C末端含7個(gè)類(lèi)表皮生長(zhǎng)因子(Epidermal Growth Factor like, EGF-like)結(jié)構(gòu)域，每個(gè)結(jié)構(gòu)域約含38個(gè)氨基酸殘基。運(yùn)行程序dotmatcher對(duì)其自身進(jìn)行比對(duì)，采用不同大小的滑動(dòng)窗口和相似性閾值，可得到不同結(jié)果(見(jiàn)圖3)。當(dāng)窗口大小和相似性閾值均為默認(rèn)值時(shí)，背景噪聲較大(見(jiàn)圖3a)；當(dāng)把窗口大小改為24個(gè)殘基，相似性閾值改為20時(shí)，可清晰顯示長(zhǎng)度為24的亮氨酸重復(fù)片段(見(jiàn)圖3b)。當(dāng)窗口大小與重復(fù)單元大小相近時(shí)，所顯示的重復(fù)區(qū)域最為清晰。當(dāng)把窗口大小改為38個(gè)殘基，相似性閾值改為20時(shí)，可清晰顯示表皮生長(zhǎng)因子結(jié)構(gòu)域(見(jiàn)圖3c)。當(dāng)保持窗口大小為38而相似性閾值改為30時(shí)，可減少背景噪聲(見(jiàn)圖3d)。閾值越大，背景噪聲越小。

圖3 果蠅體節(jié)發(fā)育基因蛋白質(zhì)序列點(diǎn)陣圖分析結(jié)果Fig.3 Dot-plot results of repeat regions in protein sequence of fruitfly midline development related gene

利用點(diǎn)陣圖程序，通過(guò)設(shè)置適當(dāng)參數(shù)，可清晰顯示串聯(lián)重復(fù)基因和重復(fù)結(jié)構(gòu)域等。重復(fù)結(jié)構(gòu)域在蛋白質(zhì)分子中較為多見(jiàn)，如鈣結(jié)合蛋白(CALB1_HUMAN)含5個(gè)長(zhǎng)度為36 aa的EF-手型(EF-hand)重復(fù)單元，膜聯(lián)蛋白A1(ANXA1_HUMAN)含4個(gè)長(zhǎng)度為72-73 aa的膜聯(lián)蛋白重復(fù)單元，抗肌萎縮蛋白(DMD_HUMAN)含24個(gè)長(zhǎng)度為100-110 aa的紅膜肽(Spectrin)重復(fù)單元。而肌聯(lián)蛋白(TITIN_HUMAN)則是由多個(gè)不同重復(fù)單元組成的巨型蛋白質(zhì)，全長(zhǎng)34 350個(gè)氨基酸，含152個(gè)長(zhǎng)度為90 aa左右的免疫球蛋白結(jié)構(gòu)域(Ig-like)、132個(gè)長(zhǎng)度為95 aa左右的III型纖聯(lián)蛋白(Fibronectin type-II)、19個(gè)長(zhǎng)度為45 aa左右的Kelch重復(fù)單元、17個(gè)長(zhǎng)度為55 aa左右的RCC1重復(fù)單元、15個(gè)長(zhǎng)度為40 aa左右的WD重復(fù)單元和14個(gè)長(zhǎng)度為34 aa左右的TPR重復(fù)單元。

4 核酸序列分析

4.1 序列組分統(tǒng)計(jì)

組分分析是序列分析中最基本的方法之一。EMBOSS中用于核酸序列組分分析的程序包括geecee, freak, wordcount和compseq等，其中g(shù)eecee和freak主要用于GC含量分析，wordcount和compseq用于統(tǒng)計(jì)四種核苷酸出現(xiàn)頻率。

4.1.1 compseq用法實(shí)例

程序compseq可用于按指定長(zhǎng)度統(tǒng)計(jì)核酸序列中不同字串出現(xiàn)頻率。所謂字串，是指一定長(zhǎng)度的不同核苷酸組合。長(zhǎng)度為2的2字串共有16種，即AA, AC, AG, AT, ……, TA, TC, TG, TT；長(zhǎng)度為3的3字串有64種，即AAA, AAC, AAG, AAT, ACA, ACC, ACG, ACT, ……, TTA, TTC, TTG, TTT；4字串、5字串和6字串分別有256、1,024和4,096種。下面以三種模式微生物為例，分別統(tǒng)計(jì)6字串出現(xiàn)的次數(shù)和與期望值的比例。

compseq ECOLI_K12.FASTA -out ECOLI_K12.COMP -word 6

compseq MYCTO_H37.FASTA -out MYCTO_H37.COMP -word 6

compseq CALSU_MB4.FASTA -out CALSU_MB4.COMP -word 6

compseq - 計(jì)算指定長(zhǎng)度核苷酸組分程序

ECOLI_K12.FASTA - 大腸桿菌K12菌株基因組序列

MYCTO_H37.FASTA - 結(jié)核分枝桿菌H37菌株基因組序列

CALSU_MB4.FASTA - 泉生熱胞菌MB4菌株基因組序列

ECOLI_K12.COMP - 大腸桿菌6字串輸出結(jié)果

MYCTO_H37.COMP - 結(jié)核分枝桿菌6字串輸出結(jié)果

CALSU_MB4.COMP - 泉生熱胞菌6字串輸出結(jié)果

-word 6 - 指定字串長(zhǎng)度為6

上述程序運(yùn)行結(jié)果每個(gè)序列都生成4 096個(gè)不同的6字串，其中有的6字串頻數(shù)很低，有的6字串頻數(shù)很高(見(jiàn)表5)。這三種模式微生物在細(xì)菌基因組學(xué)研究中具有重要地位。大腸桿菌是人類(lèi)基因組計(jì)劃指定的模式微生物，結(jié)核分枝桿菌是最早完成基因組測(cè)序的致病菌，泉生熱胞菌是我國(guó)科學(xué)家于2002年完成基因組測(cè)序的第一個(gè)細(xì)菌。

表5 三種模式微生物基因組序列中的特殊6字串Table 5 Special 6-mer in three models of bacterial genome sequences

4.1.2 freak用法實(shí)例

程序compseq用于統(tǒng)計(jì)核酸序列中不同字串出現(xiàn)頻率，而程序freak則可以圖形方式輸出不同區(qū)域GC含量。以下命令可顯示小鼠alpha血紅蛋白基因DNA序列(全長(zhǎng)1 441 bp)不同區(qū)域GC含量分布。

freak V00714.FASTA -letters"GC"-plot Y -graph svg -goutfile V00714_FREAK -window 100 -step 10

freak - 統(tǒng)計(jì)DNA序列中GC含量程序

V00714.FASTA - 小鼠alpha血紅蛋白基因序列

-letters "GC" - 顯示GC含量

-plot Y - 生成圖形文件

-goutfile V00714_FREAK - 圖形格式輸出文件名

-window 100 - 滑動(dòng)窗口大小為100

-step 10 - 步長(zhǎng)為10

從程序freak輸出結(jié)果可以看出，小鼠alpha血紅蛋白基因5’端和3’端GC含量較低，而在600-1 000 bp區(qū)域GC含量較高(見(jiàn)圖4)。

圖4 小鼠alpha血紅蛋白基因不同區(qū)域GC含量Fig.4 GC content of different regions in mouse alpha hemoglobin gene

4.2 開(kāi)放讀碼框分析

EMBOSS中整合的開(kāi)放讀碼框分析程序包括plotorf, sixpack, showorf和getorf。這四個(gè)程序可以組合使用，plotorf用圖形方式輸出DNA序列中所有可能的讀碼框，即起始密碼子和終止密碼子之間、終止密碼子和終止密碼子之間的序列，包括三條正鏈和三條負(fù)鏈；sixpack給出DNA序列所有可能的編碼方式；showorf按指定編碼鏈輸出DNA序列及其編碼的氨基酸序列；getorf用于提取讀碼框序列或所編碼的氨基酸序列，也可提取編碼區(qū)上游或下游非翻譯區(qū)序列。

本文1.3.3中介紹了getorf的用法，下面介紹sixpack和showorf的用法。

4.2.1 sixpack用法實(shí)例

以豌豆開(kāi)花后特異表達(dá)基因全長(zhǎng)mRNA序列(GenBank登錄號(hào)Y12618)為例，調(diào)用EMBOSS軟件包中的sixpack程序，輸出該序列正鏈(F1-F3)和互補(bǔ)鏈(F4-F6)序列，以及6條開(kāi)放讀碼框所編碼的氨基酸。

sixpack Y12618.FASTA Y12618.SIXPACK -outseq Y12618_ORF.FASTA

sixpack - 顯示6條讀碼框程序

Y12618.FASTA - 豌豆開(kāi)花后特異表達(dá)基因FASTA格式序列

Y12618.SIXPACK - 輸出結(jié)果讀碼框和對(duì)應(yīng)的氨基酸序列

-outseq Y12618_ORF.FASTA - FASTA格式讀碼框序列文件

運(yùn)行結(jié)果顯示，正鏈第3條讀碼框(F3)第48位有起始密碼子ATG，第1 374位和1 377位有兩個(gè)連續(xù)終止密碼子TGA和TAG，表明該序列編碼區(qū)位于正鏈48-1 373位，共1 326 bp，編碼442個(gè)氨基酸。

4.2.2 showorf用法實(shí)例

上述以豌豆開(kāi)花后特異表達(dá)基因全長(zhǎng)mRNA序列為例，調(diào)用EMBOSS中showorf程序，指定第3條讀碼框，則輸出該讀碼框序列及對(duì)應(yīng)的氨基酸序列。

showorf Y12618.FASTA Y12618.SHOWORF -frames 3

showorf - 顯示指定讀碼框程序

Y12618.FASTA - 豌豆開(kāi)花后特異表達(dá)基因FASTA格式序列

Y12618.SHOWORF - 輸出結(jié)果文件

-frames 3 - 輸出第3條讀碼框(F3)

4.3 CpG島識(shí)別

CpG島是指DNA序列中富含CG雙核苷酸的區(qū)域，其順序?yàn)镃在前，G在后。為避免誤解，常用CpG表示，即胞嘧啶3’端與鳥(niǎo)嘌呤5’端通過(guò)磷酸基團(tuán)連接。CpG島通常位于基因上游啟動(dòng)子區(qū)域300-3 000 bp區(qū)域內(nèi)，該區(qū)域的特征是核苷酸G和C含量較高，且富含CpG雙核苷酸。因此，CpG島預(yù)測(cè)結(jié)果可用來(lái)推斷某個(gè)DNA序列片段中是否存在蛋白質(zhì)編碼基因。

EMBOSS中整合的CpG島分析程序包括cpgplot和cpgreport等。程序cpgplot用于預(yù)測(cè)DNA序列中的CpG島，并以圖形方式輸出結(jié)果。程序cpgreport用于計(jì)算DNA序列中CpG雙核苷酸含量，所用方法與cpgplot有所不同，靈敏度較高，但假陽(yáng)性率也較高。

4.3.1 cpgplot用法實(shí)例

人alpha血紅蛋白基因家族分布在16號(hào)染色體短臂靠近端粒處，包括5個(gè)功能基因(zeta, mu, alpah2, alpha1和theta)以及兩個(gè)假基因(HBZP和HBA1P)。該基因家族DNA序列長(zhǎng)度為43 058 bp(GenBank登錄號(hào)Z84721)。下載FASTA格式序列并運(yùn)行cpgplot程序。

cpgplot Z84721.FASTA -window 500 -minlen 500 -minoe 0.65 -minpc 0.55 -outfile Z84721.CPGPLOT -outfeat Z84721.GFF -graph svg -goutfile Z84721_CPGPLOT

cpgplot - 顯示DNA序列中CpG島程序

Z84721.FASTA - 人alpha珠蛋白基因家族DNA序列

-window 500 - 滑動(dòng)窗口大小，默認(rèn)值200

-minlen 500 - CpG島最小長(zhǎng)度(minimum length)，默認(rèn)值100

-minoe 0.65- CpG含量平均值觀察值與期望值最小比值(minimum average observed to expected

ratio)，默認(rèn)值0.6

-minpc 0.55 - GC含量平均值最小值(minimum average percentage)，默認(rèn)值0. 5

-outfile Z84721.CPGPLOT - 輸出結(jié)果文件

-outfeat Z84721.GFF - 輸出結(jié)果文件

-goutfile Z84721_CPGPLOT - 輸出圖形結(jié)果文件

上述運(yùn)行過(guò)程中，滑動(dòng)窗口大小和CpG島長(zhǎng)度均設(shè)為500 bp，雙核苷酸CpG含量觀察值與期望值比值下限設(shè)為0.65，GC含量下限設(shè)為0.55。查看該序列條目Z84721中注釋信息，預(yù)測(cè)結(jié)果與注釋信息比較吻合。若采用系統(tǒng)給定缺省參數(shù)，預(yù)測(cè)靈敏度較高，但假陽(yáng)性率也較高。結(jié)果表明，該基因組序列片段中有5個(gè)CpG島(見(jiàn)表6)。

表6 人alpha血紅蛋白基因家族序列中的CpG島Table 6 CpG island of human alpha hemoglobin gene cluster

程序cpgplot除了輸出文本文件Z84721.CPGPLOT外，還可輸出圖形文件，以波形圖方式顯示序列不同區(qū)域CpG雙核苷酸的含量和可能的CpG島位置(見(jiàn)圖5)。

圖5 程序cpgplot分析人alpha血紅蛋白基因簇序列中的CpG島輸出結(jié)果Fig.5 Cpgplot result of human alpha hemoglobin gene cluster

4.4 密碼子分析

密碼子一共有64個(gè)。密碼子具有通用性、簡(jiǎn)并性和偏好性等特點(diǎn)。除個(gè)別特殊密碼子外，通用密碼子中ATG為起始密碼子或編碼甲硫氨酸(Met, M)；TAA, TAG和TGA為終止密碼子；其余60個(gè)密碼子編碼19種氨基酸。編碼同一氨基酸的密碼子使用頻率可能不同，即密碼子使用具有偏好性(Codon Usage Bias)。分析密碼子使用頻率及偏好性，是系統(tǒng)發(fā)育和分子演化研究的常用方法。

EMBOSS中密碼子分析程序包括cusp、syco、cai和chips等。其中cusp用于統(tǒng)計(jì)核酸序列中64種密碼子使用頻率和期望值，并給出編碼同一氨基酸的不同密碼子使用比例。syco用圖形方式顯示同義密碼子使用偏好，可用于基因預(yù)測(cè)。cai用于計(jì)算密碼子適應(yīng)指數(shù)(Codon Adaptation Index)，取值范圍為0-1；cai值越大，密碼子使用偏好性越強(qiáng)。一般說(shuō)來(lái)，表達(dá)量高的基因，其密碼子使用偏好性強(qiáng)；因此，cai值可用于預(yù)測(cè)基因表達(dá)水平高低。chips用于計(jì)算有效密碼子數(shù)(Effective Number of Codon, ENC)，其范圍為20-61。ENC值越小，密碼子使用偏好性越強(qiáng)。不同物種或同一物種的不同基因，其ENC值有所不同。

4.4.1 cusp用法實(shí)例

以豌豆開(kāi)花后特異表達(dá)基因Y12618編碼區(qū)序列為例，程序cusp運(yùn)行結(jié)果為不同密碼子使用頻率。結(jié)果表明，該編碼區(qū)序列密碼子第3位偏好使用A或T。

cusp Y12618_CDS.FASTA Y12618_CDS.CUSP

cusp - 統(tǒng)計(jì)密碼子使用頻率程序

Y12618_CDS.FASTA - 豌豆開(kāi)花后特異表達(dá)基因編碼區(qū)序列

Y12618_CDS.CUSP - 密碼子頻率輸出結(jié)果文件

4.4.2 chips用法實(shí)例

以豌豆開(kāi)花后特異表達(dá)基因Y12618編碼區(qū)序列為例，運(yùn)行chips程序，可得到有效密碼子ENC值。

chips Y12618_CDS.FASTA Y12618_CDS.CHIPS

chips - 計(jì)算有效密碼子數(shù)程序

Y12618_CDS.FASTA-豌豆開(kāi)花后特異表達(dá)基因編碼區(qū)序列

Y12618_CDS.CHIPS - 輸出結(jié)果有效密碼子ENC值

4.5 重復(fù)序列尋找

重復(fù)序列在核酸序列中十分普遍。常見(jiàn)的重復(fù)序列包括串聯(lián)重復(fù)(Tandem Repeat)和倒轉(zhuǎn)重復(fù)(Inverted Repeat)。串聯(lián)重復(fù)是指某一特定序列片段連續(xù)多次重復(fù)排列，如DNA序列中微衛(wèi)星(Microsatellite)序列、短散在重復(fù)元件(Short Interspersed Nuclear Elements, SINE)和長(zhǎng)散在重復(fù)元件(Long Interspersed Nuclear Elements, LINE)等。倒轉(zhuǎn)重復(fù)是指同一條鏈上兩個(gè)序列片段通過(guò)堿基配對(duì)反向互補(bǔ)，如microRNA前體序列中的莖環(huán)結(jié)構(gòu)(Stem-loop)。

EMBOSS中重復(fù)序列分析程序包括etandam、equicktandem、einverted和palindrome等。程序etandem和equicktandem用于尋找核酸序列中串聯(lián)重復(fù)序列，前者允許空位插入，而后者不允許插入空位，運(yùn)算速度較快。程序palindrome和einverted用于尋找倒轉(zhuǎn)重復(fù)序列，前者用于尋找較短片段的回文結(jié)構(gòu)，兩個(gè)配對(duì)序列之間可以有錯(cuò)配，但不允許空位插入；而后者用于尋找較長(zhǎng)的倒轉(zhuǎn)重復(fù)，既可以有錯(cuò)配，也可以有空位插入。

4.5.1 palindrome用法實(shí)例

以擬南芥 microRNA 172c(ath-MIR172c)為例，從microRNA數(shù)據(jù)庫(kù)下載前體序列(登錄號(hào)MI0000991)，利用seqret程序?qū)⑵滢D(zhuǎn)換成FASTA格式，并將字符U替換成T。運(yùn)行程序palindrome，用交互式方式設(shè)置參數(shù)：最小反向重復(fù)序列長(zhǎng)度22，最大反向重復(fù)序列長(zhǎng)度25，反向重復(fù)序列間最大間隔100，允許錯(cuò)配核苷酸數(shù)2，輸出結(jié)果包括互相重疊的重復(fù)序列。

palindrome ATH-MIR172C.FASTA ATH-MIR172C.PAL

Finds inverted repeats in nucleotide sequence(s)

Enter minimum length of palindrome [10]:22

Enter maximum length of palindrome [100]:25

Enter maximum gap between repeated regions [100]:

Number of mismatches allowed [0]:2

Report overlapping matches [Y]:

palindrome -尋找反向重復(fù)序列程序

ATH-MIR172C.FASTA -擬南芥microRNA 172c前體FASTA格式序列

ATH-MIR172C.PAL -運(yùn)行結(jié)果輸出文件

運(yùn)行結(jié)果表明，microRNA前體序列ath-MIR 172c中17-38/96-117位為倒轉(zhuǎn)重復(fù)序列。查看microRNA數(shù)據(jù)庫(kù)miRBase中的注釋信息，該前體序列成熟miRNA位于98-118位。

4.5.2 einverted用法實(shí)例

以果蠅性別相關(guān)基因?yàn)槔瑥腉enBank下載序列(登錄號(hào)EF565211)，運(yùn)行程序einverted，參數(shù)設(shè)置采用系統(tǒng)默認(rèn)值，空位罰分12、最小分值50、匹配分值3和錯(cuò)配分值-4。輸出結(jié)果為倒轉(zhuǎn)重復(fù)文件和FASTA格式序列文件。

einverted EF565211.FASTA EF565211.EINV -outseq EF565211_EINV.FASTA

Finds inverted repeats in nucleotide sequences

Gap penalty [12]:

Minimum score threshold [50]:

Match score [3]:

Mismatch score [-4]:

einverted - 尋找倒轉(zhuǎn)重復(fù)序列程序

EF565211.FASTA - 果蠅性別相關(guān)基因FASTA格式序列

EF565211.EINV - 輸出結(jié)果倒轉(zhuǎn)重復(fù)文件

-outseq EF565211_EINV.FASTA - 輸出結(jié)果FASTA格式文件

運(yùn)行結(jié)果表明，該序列中1 617-1 966/2 355-2 699位為倒轉(zhuǎn)重復(fù)，中間有1個(gè)長(zhǎng)度為6的空位。查看該序列注釋信息，該倒轉(zhuǎn)重復(fù)序列與果蠅性別比例抑制功能有關(guān)。

5 蛋白質(zhì)序列分析

5.1 序列組分統(tǒng)計(jì)

EMBOSS中用于蛋白質(zhì)一級(jí)結(jié)構(gòu)氨基酸序列統(tǒng)計(jì)分析的程序包括pepstats, pepinfor, wordcount和 compseq等。pepstats以文本和表格方式輸出蛋白質(zhì)序列中各種氨基酸含量，并對(duì)不同類(lèi)型氨基酸進(jìn)行統(tǒng)計(jì)，如親水氨基酸和帶電氨基酸等；pepinfor則以圖形方式顯示各種類(lèi)別氨基酸在序列不同區(qū)域的分布，如疏水性氨基酸、極性氨基酸、帶電氨基酸和芳香族氨基酸等。此外，用于核酸序列組分分析的wordcount和compseq也可用于蛋白質(zhì)序列組分分析。

5.1.1 pepstats用法實(shí)例

以水稻落?？刂苹騭h4蛋白質(zhì)產(chǎn)物(UniProt序列條目Q1PIH9_ORYSI)為例，運(yùn)行程序pepstats，則可統(tǒng)計(jì)20種氨基酸出現(xiàn)頻率。

pepstats Q1PIH9_ORYSI.FASTA Q1PIH9_ORYSI.PEPSTATS

pepstats - 統(tǒng)計(jì)蛋白質(zhì)序列中不同氨基酸出現(xiàn)頻率程序

Q1PIH9_ORYSI.FASTA - 水稻落?？刂苹虻鞍踪|(zhì)FASTA格式序列

Q1PIH9_ORYSI.PEPSTATS - 水稻落粒控制基因蛋白質(zhì)氨基酸出現(xiàn)頻率輸出結(jié)果文件

運(yùn)行結(jié)果輸出20種氨基酸頻數(shù)和百分比。水稻落?？刂苹蛩幋a蛋白質(zhì)長(zhǎng)度為390個(gè)氨基酸殘基，不同氨基酸出現(xiàn)頻率很不均勻，某些氨基酸出現(xiàn)頻率較高，如脯氨酸、丙氨酸、絲氨酸和精氨酸高于平均值，而苯丙氨酸、異亮氨酸和甲硫氨酸則低于平均值。

5.1.2 Wordcount用法實(shí)例

從以上分析可以看出，水稻落?？刂苹蚓幋a蛋白質(zhì)中有些氨基酸出現(xiàn)頻率偏高。利用程序wordcount，指定不同字長(zhǎng)，可進(jìn)一步分析水稻落?？刂苹虻鞍踪|(zhì)產(chǎn)物短片段重復(fù)序列出現(xiàn)頻率。

wordcount Q1PIH9_ORYSI.FASTA Q1PIH9_ORYSI.WORD3 -word 3

wordcount Q1PIH9_ORYSI.FASTA Q1PIH9_ORYSI.WORD4 -word 4

wordcount Q1PIH9_ORYSI.FASTA Q1PIH9_ORYSI.WORD5 -word 5

wordcount - 統(tǒng)計(jì)蛋白質(zhì)序列中指定長(zhǎng)度字串出現(xiàn)頻率程序

Q1PIH9_ORYSI.FASTA - 水稻落?？刂苹蛩幋a蛋白質(zhì)FASTA格式序列

Q1PIH9_ORYSI.WORD3 - 三肽片段出現(xiàn)頻率輸出結(jié)果文件

Q1PIH9_ORYSI.WORD4 - 四肽片段出現(xiàn)頻率輸出結(jié)果文件

Q1PIH9_ORYSI.WORD5 - 五肽片段出現(xiàn)頻率輸出結(jié)果文件

結(jié)果發(fā)現(xiàn)，該蛋白質(zhì)序列中存在大量短片段重復(fù)序列(見(jiàn)表7)。

表7 水稻落粒控制基因所編碼的蛋白質(zhì)序列中短肽重復(fù)片段Table 7 Short peptide repeats in protein sequences of rice shattering control gene

5.2 序列特征位點(diǎn)識(shí)別

EMBOSS中用于蛋白質(zhì)序列特征位點(diǎn)分析的程序包括antigenic、sigcleave和digest等，其中antigenic用于抗原決定簇分析，sigcleave用于信號(hào)肽剪切位點(diǎn)分析，digest用于酶切位點(diǎn)分析。

5.2.1 antigenic用法實(shí)例

人III型癌胚抗原(CEAM3_HUMAN)胞外區(qū)35-142位為免疫球蛋白可變結(jié)構(gòu)域。利用antigenic程序，可預(yù)測(cè)該結(jié)構(gòu)域抗原決定簇，即可能的抗體結(jié)合部位。

antigenic CEAM3_HUMAN.FASTA -sbegin 35 -send 142 CEAM3_HUMAN.ANTI-minlen 10

antigenic - 抗原決定簇預(yù)測(cè)程序

CEAM3_HUMAN.FASTA - 人III型癌胚抗原FASTA格式序列

-sbegin 35 - 預(yù)測(cè)起始位點(diǎn)35

-send 142 - 預(yù)測(cè)終止位點(diǎn)142

CEAM3_HUMAN.ANTI - 輸出結(jié)果文件

-minlen 10 - 抗原決定簇最小序列長(zhǎng)度10(默認(rèn)值為6)

預(yù)測(cè)結(jié)果表明，該蛋白質(zhì)分子可能有三個(gè)抗原決定簇，分別位于第48-66, 75-86和119-129位。

5.2.2 fuzzprot用法實(shí)例

程序fuzzprot用于尋找蛋白質(zhì)序列中的序列模體。下面我們用植物特異轉(zhuǎn)錄因子家族Squamos promoter binding protein(SBP)為例，說(shuō)明fuzzprot的用法。

植物轉(zhuǎn)錄因子數(shù)據(jù)庫(kù)收錄了17個(gè)擬南芥SBP家族基因，共30種不同轉(zhuǎn)錄本，編碼17個(gè)轉(zhuǎn)錄因子。這17個(gè)轉(zhuǎn)錄因子的DNA結(jié)合結(jié)構(gòu)域長(zhǎng)度為79個(gè)氨基酸(見(jiàn)圖6)，含兩個(gè)鋅指結(jié)構(gòu)(Zinc finger)和1個(gè)核定位信號(hào)(Nuclear localization signal, NLS)。該核定位信號(hào)的序列比較保守，富含帶正電的精氨酸(Arg, R)和賴氨酸(Lys, K)。利用以下命令可以找出核定位信號(hào)在17個(gè)轉(zhuǎn)錄因子DNA結(jié)合結(jié)構(gòu)域中的位置。

fuzzpro 17ARATH_SPLD.FASTA 17ARATH_SPLD.FUZ -pattern R[RK][RK]x(6)RR[RK][KR]-pname"NLS"

fuzzpro - 序列模體尋找程序

17ARATH_SPLD.FASTA - 擬南芥17個(gè)SBP家族轉(zhuǎn)錄因子DNA結(jié)合結(jié)構(gòu)域序列

17ARATH_SPLD.FUZ - 輸出結(jié)果文件

-pattern R[RK][RK]x(6)RR[RK][KR] - 序列模體

-pname "NLS" - 序列模體名稱(chēng)

序列模體由用戶指定，保守氨基酸用大寫(xiě)字符表示，中括號(hào)內(nèi)為可選氨基酸，x為任意氨基酸，括號(hào)中為任意氨基酸個(gè)數(shù)，此處為6(輸入括號(hào)時(shí)需要加轉(zhuǎn)義符反斜杠，否則無(wú)法正常運(yùn)行)。

圖6 擬南芥17個(gè)植物特異轉(zhuǎn)錄因子SBP家族DNA結(jié)合結(jié)構(gòu)域序列圖標(biāo)Fig.6 Sequence logo of DNA binding domain in 17 SBP plant-specific transcription factors

5.3 二級(jí)結(jié)構(gòu)分析

EMBOSS中用于蛋白質(zhì)二級(jí)結(jié)構(gòu)分析的程序包括tmap, topo, pepwheel, helixturnhelix, pepcoil和garnier等，其中tmap用于跨膜螺旋預(yù)測(cè)，topo用于顯示跨膜螺旋拓?fù)浣Y(jié)構(gòu)，helixturnhelix用于預(yù)測(cè)螺旋-轉(zhuǎn)角-螺旋構(gòu)象，pepcoil用于預(yù)測(cè)無(wú)規(guī)卷曲肽段，pepwheel用于顯示alpha螺旋輪，garnier用于二級(jí)結(jié)構(gòu)預(yù)測(cè)。

5.3.1 tmap用法實(shí)例

以豌豆內(nèi)膜蛋白(UniProt序列條目PPF1_PEA)為例，用以下命令運(yùn)行程序tmap，可預(yù)測(cè)跨膜螺旋：

tmap PPF1_PEA.FASTA PPF1_PEA.TMAP -graph svg -goutfile PPF1_PEA_TMAP

map - 跨膜螺旋預(yù)測(cè)程序

PPF1_PEA.FASTA - 豌豆內(nèi)膜蛋白FASTA格式序列

PPF1_PEA.TMAP - 預(yù)測(cè)結(jié)果輸出文件

-goutfile PPF1_PEA_TMAP - 圖形輸出文件

預(yù)測(cè)結(jié)果同時(shí)以文本格式和圖形格式輸出。結(jié)果表明，豌豆內(nèi)膜蛋白序列中有4個(gè)可能的跨膜螺旋。跨膜螺旋長(zhǎng)度為20-22個(gè)氨基酸，通常疏水性氨基酸為主。提取其中的第1個(gè)跨膜螺旋序列，保存為FASTA格式序列文件PPF1_PEA_H1.FASTA，可用pepwheel程序繪制alpha螺旋輪。

5.3.2 pepwheel用法實(shí)例

對(duì)上述tmap程序預(yù)測(cè)所得第1個(gè)alpha螺旋FASTA格式序列，用以下命令運(yùn)行pepwheel可繪制alpha螺旋輪。結(jié)果表明，該跨膜螺旋輪主要由疏水氨基酸組成。

pepwheel PPF1_PEA_H1.FASTA -graph svg -goutfile PPF1_PEA_H1_PEPWHEEL

pepwheel - 繪制alpha螺旋輪程序

PPF1_PEA_H1.FASTA - 豌豆內(nèi)膜蛋白中預(yù)測(cè)到的第1個(gè)跨膜螺旋(序列如下)

>PPF1_PEA_H1 | 111-135

SVHVPYSYGFAIILLTVIVKAATLP

-goutfile PPF1_PEA_H1_PEPWHEEL - 圖形文件名

5.3.3 garnier用法實(shí)例

以抹香鯨肌紅蛋白(PDB登錄號(hào)1MBN)、人癌胚抗原N-端結(jié)構(gòu)域(PDB登錄號(hào)2QSQ)和人泛素蛋白(PDB登錄號(hào)1UBQ)為例，從PDB數(shù)據(jù)庫(kù)分別下載這3個(gè)蛋白質(zhì)分子FASTA格式序列，分別用以下命令運(yùn)行程序garnier，即可預(yù)測(cè)得到二級(jí)結(jié)構(gòu)。

garnier 1MBN.FASTA 1MBN.GARNIER

garnier 2QSQ.FASTA 2QSQ.GARNIER

garnier 1UBQ.FASTA 1UBQ.GARNIER

garnier -二級(jí)結(jié)構(gòu)預(yù)測(cè)程序

1MBN.FASTA -抹香鯨肌紅蛋白FASTA格式序列

1MBN.GARNIER - 抹香鯨肌紅蛋白預(yù)測(cè)結(jié)果文件

2QSQ.FASTA - 人癌胚抗原N-端結(jié)構(gòu)域FASTA格式序列

2QSQ.GARNIER - 人癌胚抗原N-端結(jié)構(gòu)域預(yù)測(cè)結(jié)果文件

1UBQ.FASTA - 人泛素蛋白FASTA格式序列

1UBQ.GARNIER -人泛素蛋白預(yù)測(cè)結(jié)果文件

預(yù)測(cè)結(jié)果中字母H表示alpha螺旋(Helix)，E表示beta折疊(Extended)、T表示轉(zhuǎn)角(Turn)、C表示無(wú)規(guī)卷曲(Coil)；下劃線表示預(yù)測(cè)準(zhǔn)確區(qū)域。這3種蛋白質(zhì)分子的三維空間結(jié)構(gòu)均已由實(shí)驗(yàn)測(cè)定，利用蛋白質(zhì)結(jié)構(gòu)顯示和分析軟件可觀察這3種蛋白質(zhì)分子的實(shí)際構(gòu)象。肌紅蛋白全長(zhǎng)153個(gè)氨基酸殘基，由8股alpha螺旋組成，按N-端到C-端順序編號(hào)為A-H；癌胚抗原N-端結(jié)構(gòu)域全長(zhǎng)110個(gè)氨基酸殘基，由9個(gè)beta折疊片組成，按N-端到C-端順序?yàn)锳-G(C和D之間另有C’和C”兩個(gè)beta折疊片)。泛素蛋白全長(zhǎng)76個(gè)氨基酸殘基，由5個(gè)beta折疊(A, B, D, E, G)和2個(gè)alpha螺旋(C, F)組成。與實(shí)際構(gòu)象比較表明，預(yù)測(cè)結(jié)果有一定誤差。二級(jí)結(jié)構(gòu)預(yù)測(cè)方法很多，目前預(yù)測(cè)精度為70%-80%。除EMBOSS中整合的garnier外，許多網(wǎng)站提供在線預(yù)測(cè)工具。讀者可嘗試不同程序，比較預(yù)測(cè)結(jié)果。

6 總結(jié)和討論

6.1 使用說(shuō)明

以上我們介紹了EMBOSS軟件包中一些常用程序。經(jīng)過(guò)十年多開(kāi)發(fā)，EMBOSS已成為核酸和蛋白質(zhì)序列分析常用軟件包，為廣大生物學(xué)工作者廣泛使用。EMBOSS軟件包功能齊全、用途廣泛。選修“實(shí)用生物信息技術(shù)”課程的同學(xué)，在學(xué)習(xí)EMBOSS軟件包中常用程序后，編寫(xiě)了以下順口溜。

EMBOSS軟件包，包羅萬(wàn)象真的好，

核酸蛋白都適用，功能強(qiáng)大效率高。

比對(duì)進(jìn)化設(shè)引物，翻譯酶切找重復(fù)，

程序名稱(chēng)不記得，wossname幫你找。

程序命令不會(huì)用，tfm 把你教。

參數(shù)設(shè)置技巧高，點(diǎn)點(diǎn)滴滴要記牢，

EMBOSS是法寶，活學(xué)活用不愁了。

要熟練使用EMBOSS軟件包中的程序，首先必須熟悉分子生物學(xué)基本概念，如中心法則和序列-結(jié)構(gòu)-功能關(guān)系等；掌握必要的分子生物學(xué)基礎(chǔ)知識(shí)，如基因結(jié)構(gòu)、啟動(dòng)子、外顯子、內(nèi)含子、編碼序列、RNA二級(jí)結(jié)構(gòu)、蛋白質(zhì)結(jié)構(gòu)層次、一級(jí)結(jié)構(gòu)序列特征和二級(jí)結(jié)構(gòu)單元，以及序列模體、結(jié)構(gòu)域、蛋白質(zhì)家族和蛋白質(zhì)功能等。

選擇合適的程序、設(shè)置恰當(dāng)?shù)膮?shù)，是正確、高效使用EMBOSS軟件包的關(guān)鍵。除了深入了解所研究問(wèn)題的生物學(xué)背景外，也應(yīng)搞清輸入數(shù)據(jù)的種類(lèi)、格式，掌握各種參數(shù)的含義，對(duì)所用程序的算法有所了解。同樣一個(gè)問(wèn)題，使用不同程序，運(yùn)行結(jié)果就可能不同；即使是同一個(gè)程序，參數(shù)不同，結(jié)果也可能不同。熟練使用EMBOSS軟件包提供的三個(gè)幫助程序wossname、tfm和seealso，深入理解各個(gè)程序的功能和可供設(shè)置的參數(shù)，可以在程序使用過(guò)程中起到事半功倍的效果。

需要說(shuō)明的是，EMBOSS軟件包啟動(dòng)時(shí)，人類(lèi)基因組計(jì)劃尚未完成，基因組數(shù)據(jù)分析剛剛開(kāi)始。因此，EMBOSS不是組學(xué)數(shù)據(jù)分析軟件，而是針對(duì)單個(gè)或多個(gè)蛋白質(zhì)或核酸序列分析工具，其功能相當(dāng)于基于個(gè)人計(jì)算機(jī)的共享軟件BioEdit或商業(yè)軟件DNAStar和MacVector等。從事基因組和轉(zhuǎn)錄組等高通量數(shù)據(jù)分析的讀者，可選擇Bowtie, BWA, TopHat/Cufflinks等軟件。此外，EMBOSS軟件包是單個(gè)程序的集成，各個(gè)程序之間并無(wú)聯(lián)系，而后來(lái)開(kāi)發(fā)的Galaxy平臺(tái)，則將某些工具整合而形成互相關(guān)聯(lián)的分析流程。

與所有計(jì)算機(jī)軟件均可能存在“bug”一樣，EMBOSS軟件包中某些程序在運(yùn)行時(shí)結(jié)果可能有誤。例如，點(diǎn)陣圖程序dotmatcher和dottup在比較兩個(gè)不同序列時(shí)，坐標(biāo)軸顯示有誤。此外，由于近年來(lái)UniProt數(shù)據(jù)庫(kù)格式有所調(diào)整，序列提取程序extractfeat在處理UniProt/Swiss-Prot格式輸入文件時(shí)，得不到正確結(jié)果。讀者可自行修改源代碼改正錯(cuò)誤，或與Peter Rice聯(lián)系。

6.2 程序列表

2016年發(fā)布的EMBOSS 6.6.0版包括300多個(gè)程序，本文介紹的程序只是其中一小部分。為便于讀者查詢，我們按類(lèi)別列出其中部分常用程序(見(jiàn)表8)。

表8 EMBOSS軟件包常用程序Table 8 List of commonly used programs in EMBOSS

pepstats一級(jí)結(jié)構(gòu)分析統(tǒng)計(jì)蛋白質(zhì)序列中不同種氨基酸出現(xiàn)頻率pepinfo一級(jí)結(jié)構(gòu)分析用圖形方式顯示蛋白質(zhì)序列不同氨基酸分布特征antigenic一級(jí)結(jié)構(gòu)分析預(yù)測(cè)蛋白質(zhì)序列中抗原決定簇sigcleave一級(jí)結(jié)構(gòu)分析尋找蛋白質(zhì)序列中信號(hào)肽切割位點(diǎn)digest一級(jí)結(jié)構(gòu)分析尋找蛋白質(zhì)序列中蛋白酶酶切位點(diǎn)fuzzprot一級(jí)結(jié)構(gòu)分析尋找蛋白質(zhì)序列中序列模體Tmap二級(jí)結(jié)構(gòu)分析預(yù)測(cè)蛋白質(zhì)序列中的跨膜螺旋topo二級(jí)結(jié)構(gòu)分析顯示跨膜螺旋拓?fù)浣Y(jié)構(gòu)pepwheel二級(jí)結(jié)構(gòu)分析繪制alpha螺旋輪pepcoil二級(jí)結(jié)構(gòu)分析預(yù)測(cè)無(wú)規(guī)卷曲區(qū)域helixturn-helix二級(jí)結(jié)構(gòu)分析螺旋-轉(zhuǎn)角-螺旋序列模體分析garnier二級(jí)結(jié)構(gòu)分析蛋白質(zhì)序列二級(jí)結(jié)構(gòu)預(yù)測(cè)

EMBOSS網(wǎng)站列出了所有程序的名稱(chēng)和用途，也可用wossname命令按功能分類(lèi)或字母表順序列出所有程序。

wossname -search“”

按功能分類(lèi)列出所有程序

wossname -search“”-alphabetic

按字母表順序列出所有程序

除了EMBOSS開(kāi)發(fā)團(tuán)隊(duì)自行編寫(xiě)的程序外，EMBOSS還整合了不少其它常用生物信息軟件包，如基于隱馬爾可夫模型的蛋白質(zhì)結(jié)構(gòu)域序列譜構(gòu)建和結(jié)構(gòu)域識(shí)別軟件包HEMMER、系統(tǒng)發(fā)育分析軟件包Phylip及RNA二級(jí)結(jié)構(gòu)分析和預(yù)測(cè)軟件包VIENNA等。限于篇幅，本文未介紹這些軟件包中程序的用法。

6.3 回顧和展望

EMBOSS項(xiàng)目始于上世紀(jì)九十年代，初始宗旨是取代序列分析商業(yè)軟件包GCG。上世紀(jì)八十年代，美國(guó)威斯康辛大學(xué)遺傳計(jì)算團(tuán)隊(duì)(Genetic Computing Group, GCG)開(kāi)發(fā)了基于UNIX的序列分析軟件并商業(yè)化。該軟件包整合了序列比對(duì)、酶切位點(diǎn)分析等許多工具，在美國(guó)和歐洲等西方國(guó)家十分流行[9]。早期的GCG軟件包源代碼公開(kāi)，用戶可以修改和整合自己的程序。九十年代末，GCG軟件包不再公開(kāi)源代碼。為避免GCG商業(yè)軟件的限制，歐洲分子生物學(xué)網(wǎng)絡(luò)組織EMBnet啟動(dòng)了EMBOSS項(xiàng)目，并很快取代了GCG。有關(guān)EMBOSS項(xiàng)目啟動(dòng)和實(shí)施過(guò)程，以及EMBnet的詳細(xì)情況，請(qǐng)參閱本刊擬于年內(nèi)發(fā)表的“EMBOSS和EMBnet”一文。

本世紀(jì)初，Peter Rice領(lǐng)導(dǎo)的EMBOSS研發(fā)團(tuán)隊(duì)受聘于歐洲生物信息學(xué)研究所，完成了該軟件包的主要開(kāi)發(fā)。2009年，EMBOSS項(xiàng)目曾得到英國(guó)生物技術(shù)和生命科學(xué)研究委員會(huì)(Biotechnology and Biological Science Research Council, BBSRC)資助，繼續(xù)進(jìn)行開(kāi)發(fā)。目前，該軟件包開(kāi)發(fā)項(xiàng)目已經(jīng)結(jié)束，由英國(guó)transSMART基金會(huì)Peter Rice負(fù)責(zé)維護(hù)。顯然，作為開(kāi)源軟件，EMBOSS的進(jìn)一步開(kāi)發(fā)，需要得到生物信息軟件開(kāi)發(fā)人員和廣大用戶的支持。對(duì)該軟件包開(kāi)發(fā)感興趣的讀者可與Peter Rice直接聯(lián)系，聯(lián)系方式請(qǐng)參閱補(bǔ)充材料1。

致 謝

感謝楊德昌安裝和調(diào)試EMBOSS軟件包，顏林林改正點(diǎn)陣圖程序中的錯(cuò)誤。感謝樊麗編寫(xiě)的EMBOSS順口溜。金錄佳、李宏博和趙坤認(rèn)真閱讀并校正了初稿中多處文字錯(cuò)誤。感謝匿名審稿人寶貴的修改意見(jiàn)。感謝中國(guó)科學(xué)院北京基因組研究所(國(guó)家生物信息中心)對(duì)EMBOSS網(wǎng)絡(luò)教程提供的支持。感謝北京大學(xué)生命科學(xué)學(xué)院、中國(guó)農(nóng)業(yè)科學(xué)院研究生院和中國(guó)科學(xué)院大學(xué)生命科學(xué)學(xué)院多年來(lái)對(duì)“實(shí)用生物信息技術(shù)”課程的支持。