EMBOSS軟件包序列分析程序應用實例

羅 靜 初

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

1 概述

1.1 簡介

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

2001年起，筆者在北京大學開設“實用生物信息技術”研究生課程[4]，曾介紹EMBOSS軟件包中部分常用程序。2020年新冠病毒疫情期間，為北京大學生物信息學專業本科生開設“Linux基礎及其在生物信息領域中的應用”線上課程，較為系統地介紹了EMBOSS軟件包主要程序。

為推廣EMBOSS軟件包在生物信息學研究中的應用，本文基于教學中的一些實例，介紹EMBOSS軟件包中主要程序的用途、用法，及其運行結果的生物學意義。文中所舉實例的蛋白質序列源自國際蛋白質數據庫UniProt[5]，核酸序列源自美國國家生物技術信息中心(National Center for Biotechnology Information, NCBI)的核酸序列數據庫GenBank和參考序列數據庫RefSeq。文中對所舉實例的生物學背景作了簡要說明，并給出相關文獻，具體細節可參閱“實用生物信息技術”課程(Applied Bioinformatics Course, ABC)教學網站(https://bigd.big.ac.cn/education/)。獲中國科學院北京基因組研究所大數據中心支持，本文補充材料已上載到該所國家生物信息中心網站(https://bigd.big.ac.cn/education/)。文中所涉及的網站請參閱補充材料1，所舉實例中的輸入數據可通過補充材料2中的鏈接下載，運行結果可通過補充材料3中的鏈接查看。計劃將文中主要程序實例改編成網絡教程，并不斷進行擴充和更新。

1.2 用戶界面

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

1.2.1 命令行方式

EMBOSS基于UNIX操作系統開發，所有程序均可在Linux系統上用命令行方式執行。國內外不少生物信息學相關科研機構、高等學校和公司企業的Linux服務器上裝有EMBOSS軟件包。讀者可向服務器管理員申請Linux系統賬號，登錄后即可通過命令行方式運行該軟件包中的程序。MacOS系統用戶可下載EMBOSS源代碼并編譯、安裝該軟件包。基于Windows系統版本mEMBOSS可在磁盤操作系統(Disk Operation System, DOS)環境下以命令行方式運行。

1.2.2 圖形界面

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

1.2.3 瀏覽器界面

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

1.3 運行模式

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

1.3.1 交互式

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

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

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

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

上述運行過程說明如下。

1)用戶輸入程序名needle，按回車鍵后開始運行。

2)系統給出所運行程序的簡單說明Needleman-Wunsch global alignment of two sequences，并提示用戶輸入用于比對的序列文件名Input sequence。

3)輸入第一個序列文件名HBA_HUMAN.FASTA后系統提示輸入第二個序列文件名。

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

5)系統提示輸入延伸空位罰分值Gap extension penalty，并給出默認值0.5，按回車鍵采用默認值。

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

1.3.2 參數式

上述交互式運行模式適用于初學者，用戶可根據提示信息確定需要輸入的文件名和參數，許多情況下可先使用默認值，在分析運行結果后再次運行時則可適當調節參數，以得到更好的結果。對于熟練用戶，則可采用參數模式運行程序，即在命令行中直接給出輸入文件名、輸出文件名和程序參數。例如，用needle程序對人、小鼠和大鼠三種哺乳動物alpha血紅蛋白進行兩兩比對，則可分別采用以下命令。

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程序運行過程說明如下。

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

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

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

上述命令的運行結果包括比對分值、相同位點數及百分比、相同及相似位點數及百分比(見表1)。

表1 人、小鼠和大鼠alpha血紅蛋白兩兩比對結果

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

1.3.3 菜單式

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

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

上述豌豆開花后特異表達基因序列分析的第一步，需要提取mRNA序列中自起始密碼子ATG到終止密碼子之間的編碼區序列。該序列已遞交NCBI核酸序列數據庫GenBank(登錄號Y12618)，并作了初步注釋，序列全長1 523個核苷酸，包括第48-1 373位編碼區、5’端和3’端非翻譯區(Untranslated region, UTR)。

從核酸序列數據庫GenBank中下載FASTA格式序列文件Y12618.FASTA，采用菜單式運行編碼區提取程序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)

…… (選項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

…… (選項2-6省略)

Type of output [0]: 1

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

上述運行過程說明如下。

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

2)輸入豌豆開花后特異表達基因mRNA序列Y12618.FASTA。

3)系統顯示0-23種可選遺傳密碼表，按回車鍵選擇默認通用密碼表。

4)系統顯示最小讀碼框長度，默認為30，輸入1 000，僅獲取長度大于1 000 bp的編碼區序列。

5)系統顯示最大讀碼框長度，按回車鍵選擇默認值1 000 000。

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

1.4 三個幫助程序

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

1.4.1 wossname

第一個程序為wossname，其含義為What’s the name，可通過輸入關鍵詞查找特定用途的程序名稱。例如，可用以下命令找到所有點陣圖(Dot-plot)序列比對程序。

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

第二個幫助程序為tfm，其含義為The file manual，可用來顯示某個程序的使用方法和可選參數，內容十分詳盡，命令參數部分列出可供用戶選擇的所有參數及其數據類型，并給出默認值和可選范圍(見表2)。

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

命令tfm的功能與Linux系統man命令的功能類似，詳細列出某程序所有幫助信息。此外，也可以在程序運行時用-help參數列出該程序簡略信息，包括EMBOSS軟件包的版本。

1.4.3 seealso

第三個幫助程序為seealso，即英文See also，其含義為列出EMBOSS軟件包中與某程序相關的其它程序。例如，以下命令列出與點陣圖程序dotmatcher相關的其它點陣圖程序。

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軟件包整合的程序幾乎涵蓋了序列分析的所有方面。本文按功能分類，簡要介紹部分常用程序，包括序列處理程序、序列比對程序、核酸序列分析程序和蛋白質序列分析程序，對其中一些具有代表性的程序結合實例給出具體操作步驟和運行結果，并對運行結果的生物學意義略加說明。

序列處理是序列分析的基礎，下面我們分別介紹最為常用的格式轉換、序列提取和序列變換三類序列處理程序。

2.1 格式轉換

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

EMBOSS包括多個序列格式轉換程序，此處介紹最為常用的seqret。該程序是EMBOSS軟件包開發的第一個程序，除了格式轉換外，還有其它許多功能。

2.1.1 seqret用法實例

從NCBI核酸序列數據庫下載GenBank格式豌豆開花后特異表達基因mRNA序列(登錄號Y12618)，以Y12618.GB為文件名保存。可用seqret程序將其轉換成FASTA格式。采用參數式方法，直接在命令行指定輸入文件Y12618.GB和輸出文件Y12618.FASTA。

seqret Y12618.GB Y12618.FASTA

seqret - 格式轉換程序

Y12618.GB - GenBank格式豌豆開花后特異表達基因mRNA序列文件

Y12618.FASTA - FASTA格式輸出結果文件

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

seqret PPF1_PEA.SW PPF1_PEA.FASTA

seqret - 格式轉換程序

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

PPF1_PEA.FASTA - FASTA格式輸出結果文件

2.2 序列提取

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

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

這類程序中最為常用的有coderet和extractfeat。coderet用法比較簡單，可用于提取mRNA、編碼序列和所編碼的蛋白質序列。extractfeat功能十分強大，用法也比較靈活，可提取更多種類的子序列，包括外顯子、內含子、重復序列和多聚腺苷酸信號等。此外，通過設定特征表類型(Type)、標簽(Tag)和標簽值(Value)等參數，也可提取用戶指定的一個或幾個特定子序列。

2.2.1 coderet用法實例

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

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]:

調用coderet程序，輸入小鼠alpha血紅蛋白GenBank格式文件V00714.GB，程序以GenBank序列條目中基因座(LOCUS)名稱v00714為默認輸出文件名，將輸出結果分別保存到5個文件中，分別為列表文件(v00714.coderet)、編碼區序列文件(v00714.cds)、mRNA序列文件(v00714.mrna)、所編碼的蛋白質序列文件(v00714.prot)和非編碼區序列文件(v00714.noncoding)。按EMBOSS軟件包習慣，默認輸出文件名用小寫字母表示。

2.2.2 extractfeat用法實例

以上述小鼠alpha血紅蛋白編碼基因為例，根據序列注釋信息，該基因包括兩個內含子，分別位于501-622和827-960，用以下命令可提取這兩個內含子的序列：

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

extractfeat - 子序列提取程序

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

V00714.INTRON - 輸出結果內含子序列

-type "intron" - 指定提取注釋信息中的內含子intron

2.3 序列變換

EMBOSS軟件包中的序列變換程序包括reveseq, msbar和 shuffleseq等。程序revseq將已知序列轉換成反向互補序列，msbar對已知序列進行突變，shuffleseq則用于產生隨機序列。

程序msbar可用于對已知序列進行單點或多點隨機突變，突變方式可以是替換、插入或刪除，突變位點可以是單個核苷酸點突變(Point mutation)，也可插入或刪除一個序列片段(Block mutation)，還可插入或刪除一個密碼子(Codon mutation)。程序運行默認方式為交互式，即屏幕顯示交互菜單，用戶可以自行選擇突變次數(Number of times)，也可按上述三種不同突變時，選擇插入(Insertion)、刪除(Deletion)、替換(Substitution)、復制(Duplication)和移動(Move)等不同突變方式。以豌豆開花后特異表達基因編碼區序列為例，以下命令對該序列進行1次單點插入突變、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個選項與點突變相同，略)

Types of block mutations to perform [0]: 3

Codon mutation operations

(6個選項與點突變相同，略)

Types of codon mutations to perform [0]: 4

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

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

msbar - 序列變換程序

Y12618_CDS.FASTA - FASTA格式豌豆開花后特異表達基因編碼區序列

Y12618_CDS_NEW.FASTA - 隨機突變輸出結果

-count 1 - 突變次數1次

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

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

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

突變結果可用雙序列比對程序needle驗證。由于程序msbar基于隨機突變，即使運行時設定的參數完全一致，兩次運行結果并不相同。

2.4 序列顯示

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

2.4.1 infoseq用法實例

程序infoseq簡單實用，可用于顯示序列名稱、格式及登錄號等基本信息，并可統計序列長度。對于核酸序列，還能統計GC含量。可用以下命令顯示豌豆開花后特異表達基因的基本信息。

infoseq Y12618.FASTA -outfile Y12618.INFO

infoseq - 序列信息顯示程序

Y12618.FASTA - FASTA格式豌豆開花后特異表達基因mRNA序列

Y12618.INFO - 輸出結果文件

程序infoseq既可用于GenBank和Swiss-Prot等格式，也可用于FASTA格式；既可用于單個序列，也可用于多序列文件。可用以下命令顯示12個人源癌胚抗原(Carcinoembryonic Antigen, CEA)蛋白質分子的基本信息[8]。

infoseq 12HUMAN_CEA.FASTA -outfile 12HUMAN_CEA.INFO

infoseq - 序列信息顯示程序

12HUMAN_CEA.FASTA - 12個FASTA格式人癌胚抗原蛋白質序列

12HUMAN_CEA.INFO - 輸出結果文件

2.4.2 showseq用法實例

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

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

showseq - 核酸子序列顯示程序

Y12618.GB - GenBank格式豌豆開花后特異表達基因序列

Y12618.SHOWSEQ - 輸出結果文件

-sbegin 1 - 指定子序列開始位點為1

-send 120 - 指定子序列終止位點為120

-format 3 - 指定輸出格式

-upper "48-120" - 指定48-120位用大寫字母表示

2.4.3 showfeat用法實例

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

showfeat V00714.GB V00714.SHOWFEAT

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

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

V00714.SHOWFEAT - 輸出結果文件

3 序列比對

序列比對在生物信息學中占有重要地位，是核酸和蛋白質序列分析的基礎。EMBOSS軟件包整合了十多個序列比對程序，包括雙序列比對、多序列比對、數據庫搜索，以及基于點陣圖的可視化序列比對等。

3.1 雙序列比對

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

3.1.1 needle用法實例

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

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

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

UniProt/Swiss-Prot蛋白質數據庫中收錄了12個人源癌胚抗原蛋白質CEA(見圖1，根據德國Ludwig-Maximilians大學Zimmermann教授CEA網站改編)。圖中標有N的結構域為免疫球蛋白可變結構域IgV；標有A和B的結構域為免疫球蛋白恒定結構域IgC，各分3種亞型(A1, A2, A3和B1, B2, B3)。嵌入磷脂雙層膜的箭頭表示糖基磷脂酰肌醇(Glycosylphosphatidylinositol, GPI)膜結合位點，穿過磷脂雙層膜的螺旋表示跨膜螺旋(Transmembrane helix, TMH)。

人的III型(CEAM3_HUMAN)和IV型(CEAM4_HUMAN)CEA分子長度接近，各含1個可變結構域、1個跨膜螺旋區和1個膜內區。用全局比對程序needle對其進行序列比對，命令如下。

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

needle-EMBOSS程序，用于全局序列比對

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

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

CEAM3-CEAM4.NEEDLE-輸出結果文件

-gapo 20- 起始空位罰分

-gape 2-延伸空位罰分

分析比對結果可以發現，這兩個亞型的CEA分子具有較高的相似性，僅在C-端有1個長度為8個氨基酸殘基的插入。上述命令中，起始空位罰分設為20，而不用默認值10，延伸空位罰分設為2，而不用默認值0.5，可用來避免不必要的插入或刪除。

圖1 UniProt/Swiss-Prot數據庫中12個人源癌胚抗原蛋白質分子Fig.1 Twelve human CEA proteins in UniProt/Swiss-Prot

3.1.2 water用法實例

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

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

water-雙序列局部比對程序

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

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

CEAM3-CEAM5.WATER-輸出結果文件

-sbegin 35-指定比對序列起始位點35

-send 142-指定比對序列終止位點142

-gapo 20-起始空位罰分

-gape 2-延伸空位罰分

比對結果表明，這兩個CEA分子N-端可變結構域序列具有較高相似性。

3.2 多序列比對

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

3.2.1 emma用法實例

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

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

可用以下命令對上述9個血紅蛋白進行多序列比對，除輸出FASTA格式序列比對文件外，同時輸出Newick格式分支圖文件，可用MEGA軟件顯示其樹形分支結構。

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

emma - 多序列比對程序

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

9HUMAN_HB.ALN - FASTA格式輸出結果文件

9HUMAN_HB.DND - Newick格式輸出結果文件

程序emma有許多可調參數，包括計分矩陣、空位罰分、比對方式和輸出格式等，可用菜單運行模式，即運行時加-options參數，即可指定上述參數的值。

3.2.2 edialign用法實例

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

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

edialign - 多序列比對程序

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

12HUMAN_CEA.EDIA - edialign格式比對輸出結果文件

12HUMAN_CEA.ALN - FASTA格式比對輸出結果文件

輸出結果保存到兩個文件中，12HUMAN_CEA.EDIA是多序列比對格式，比對結果中保守區域用大寫字母表示，每個位點標有數字0-9，數字越大，保守性越高。12HUMAN_CEA.ALN為FASTA格式的比對結果。

3.3 點陣圖

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

3.3.1 dottup用法實例

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

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

以下是利用EMBOSS軟件包中點陣圖程序dottup分別對河豚魚基因組序列片段進行比對的命令和所用參數。

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

dottup - 繪制點陣圖程序

-graph svg - 輸出結果圖形格式為SVG

A164138.FASTA - 河豚魚基因組片段DNA序列

-goutfile AF164138 - 輸出結果圖形文件名

-gtitle 'Cosmid' - 輸出結果圖形標題

-gsubtitle 'AF164138' - 輸出結果圖形副標題

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

上述命令輸出結果顯示兩條與對角線平行的線段(見圖2a)，表明該基因組序列片段5’端有兩個長度約為13kb相似片段，即兩個串聯重復多藥耐藥基因MDR。

用以下命令，設定比對范圍，則可進一步確定這兩個串聯重復基因的相似性。

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個序列終止位點為13 000，起始位點默認為1

-sbegin 13001 - 指定第2個序列起始位點為13 001

-send 26000 - 指定第2個序列終止位點為26 000

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

查看該基因組序列注釋信息，發現這兩個基因由20多個外顯子組成。利用coderet程序，提取編碼序列，運行dottup程序，比較這兩個編碼序列的相似性。

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個基因編碼序列

AF164138_CDS_2.FASTA - 第2個基因編碼序列

-word 8 - 序列比對時滑動窗口大小，默認為10，此處取8，以增加靈敏度

輸出結果(見圖2c)顯示，編碼區序列相似性比全長基因序列相似性更高。運行dottup程序，可進一步比較所編碼蛋白質序列相似性。

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個基因編碼蛋白質序列

AF164138_PRO_2.FASTA - 第2個基因編碼蛋白質序列

-word 6 - 序列比對時滑動窗口大小，缺省為10，此處取6，以增加靈敏度

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

3.3.2 Dotmatcher用法實例

點陣圖程序dottup多用于核酸序列，而dotmatcher則可用于蛋白質序列。下面以果蠅體節發育相關基因為例，說明利用dotmatcher顯示序列中的重復片段。該基因編碼長度為1 504個氨基酸的蛋白質(UniProt序列條目名SLIT_DROME)。可用以下dotmatcher命令和參數，得到不同輸出結果。

圖2 河豚魚多藥耐藥基因點陣圖序列比對結果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 - 繪制蛋白質序列點陣圖程序

SLIT_DROME.FASTA - 果蠅體節發育相關基因蛋白質序列

-window 10 - 滑動窗口大小為默認值10個氨基酸殘基

-threshold 23 - 相似性閾值為默認值23

-window 24 - 滑動窗口設為24個氨基酸殘基

-threshold 20 - 相似性閾值設為20

-window 38 - 滑動窗口設為38個氨基酸殘基

-threshold 20 - 相似性閾值設為20

-window 38 - 滑動窗口設為38個氨基酸殘基

-threshold 30 - 相似性閾值設為30

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

圖3 果蠅體節發育基因蛋白質序列點陣圖分析結果Fig.3 Dot-plot results of repeat regions in protein sequence of fruitfly midline development related gene

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

4 核酸序列分析

4.1 序列組分統計

組分分析是序列分析中最基本的方法之一。EMBOSS中用于核酸序列組分分析的程序包括geecee, freak, wordcount和compseq等，其中geecee和freak主要用于GC含量分析，wordcount和compseq用于統計四種核苷酸出現頻率。

4.1.1 compseq用法實例

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

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 - 計算指定長度核苷酸組分程序

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

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

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

ECOLI_K12.COMP - 大腸桿菌6字串輸出結果

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

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

-word 6 - 指定字串長度為6

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

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

4.1.2 freak用法實例

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

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

freak - 統計DNA序列中GC含量程序

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

-letters "GC" - 顯示GC含量

-plot Y - 生成圖形文件

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

-window 100 - 滑動窗口大小為100

-step 10 - 步長為10

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

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

4.2 開放讀碼框分析

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

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

4.2.1 sixpack用法實例

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

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

sixpack - 顯示6條讀碼框程序

Y12618.FASTA - 豌豆開花后特異表達基因FASTA格式序列

Y12618.SIXPACK - 輸出結果讀碼框和對應的氨基酸序列

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

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

4.2.2 showorf用法實例

上述以豌豆開花后特異表達基因全長mRNA序列為例，調用EMBOSS中showorf程序，指定第3條讀碼框，則輸出該讀碼框序列及對應的氨基酸序列。

showorf Y12618.FASTA Y12618.SHOWORF -frames 3

showorf - 顯示指定讀碼框程序

Y12618.FASTA - 豌豆開花后特異表達基因FASTA格式序列

Y12618.SHOWORF - 輸出結果文件

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

4.3 CpG島識別

CpG島是指DNA序列中富含CG雙核苷酸的區域，其順序為C在前，G在后。為避免誤解，常用CpG表示，即胞嘧啶3’端與鳥嘌呤5’端通過磷酸基團連接。CpG島通常位于基因上游啟動子區域300-3 000 bp區域內，該區域的特征是核苷酸G和C含量較高，且富含CpG雙核苷酸。因此，CpG島預測結果可用來推斷某個DNA序列片段中是否存在蛋白質編碼基因。

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

4.3.1 cpgplot用法實例

人alpha血紅蛋白基因家族分布在16號染色體短臂靠近端粒處，包括5個功能基因(zeta, mu, alpah2, alpha1和theta)以及兩個假基因(HBZP和HBA1P)。該基因家族DNA序列長度為43 058 bp(GenBank登錄號Z84721)。下載FASTA格式序列并運行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 - 滑動窗口大小，默認值200

-minlen 500 - CpG島最小長度(minimum length)，默認值100

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

ratio)，默認值0.6

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

-outfile Z84721.CPGPLOT - 輸出結果文件

-outfeat Z84721.GFF - 輸出結果文件

-goutfile Z84721_CPGPLOT - 輸出圖形結果文件

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

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

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

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

4.4 密碼子分析

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

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

4.4.1 cusp用法實例

以豌豆開花后特異表達基因Y12618編碼區序列為例，程序cusp運行結果為不同密碼子使用頻率。結果表明，該編碼區序列密碼子第3位偏好使用A或T。

cusp Y12618_CDS.FASTA Y12618_CDS.CUSP

cusp - 統計密碼子使用頻率程序

Y12618_CDS.FASTA - 豌豆開花后特異表達基因編碼區序列

Y12618_CDS.CUSP - 密碼子頻率輸出結果文件

4.4.2 chips用法實例

以豌豆開花后特異表達基因Y12618編碼區序列為例，運行chips程序，可得到有效密碼子ENC值。

chips Y12618_CDS.FASTA Y12618_CDS.CHIPS

chips - 計算有效密碼子數程序

Y12618_CDS.FASTA-豌豆開花后特異表達基因編碼區序列

Y12618_CDS.CHIPS - 輸出結果有效密碼子ENC值

4.5 重復序列尋找

重復序列在核酸序列中十分普遍。常見的重復序列包括串聯重復(Tandem Repeat)和倒轉重復(Inverted Repeat)。串聯重復是指某一特定序列片段連續多次重復排列，如DNA序列中微衛星(Microsatellite)序列、短散在重復元件(Short Interspersed Nuclear Elements, SINE)和長散在重復元件(Long Interspersed Nuclear Elements, LINE)等。倒轉重復是指同一條鏈上兩個序列片段通過堿基配對反向互補，如microRNA前體序列中的莖環結構(Stem-loop)。

EMBOSS中重復序列分析程序包括etandam、equicktandem、einverted和palindrome等。程序etandem和equicktandem用于尋找核酸序列中串聯重復序列，前者允許空位插入，而后者不允許插入空位，運算速度較快。程序palindrome和einverted用于尋找倒轉重復序列，前者用于尋找較短片段的回文結構，兩個配對序列之間可以有錯配，但不允許空位插入；而后者用于尋找較長的倒轉重復，既可以有錯配，也可以有空位插入。

4.5.1 palindrome用法實例

以擬南芥 microRNA 172c(ath-MIR172c)為例，從microRNA數據庫下載前體序列(登錄號MI0000991)，利用seqret程序將其轉換成FASTA格式，并將字符U替換成T。運行程序palindrome，用交互式方式設置參數：最小反向重復序列長度22，最大反向重復序列長度25，反向重復序列間最大間隔100，允許錯配核苷酸數2，輸出結果包括互相重疊的重復序列。

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 -尋找反向重復序列程序

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

ATH-MIR172C.PAL -運行結果輸出文件

運行結果表明，microRNA前體序列ath-MIR 172c中17-38/96-117位為倒轉重復序列。查看microRNA數據庫miRBase中的注釋信息，該前體序列成熟miRNA位于98-118位。

4.5.2 einverted用法實例

以果蠅性別相關基因為例，從GenBank下載序列(登錄號EF565211)，運行程序einverted，參數設置采用系統默認值，空位罰分12、最小分值50、匹配分值3和錯配分值-4。輸出結果為倒轉重復文件和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 - 尋找倒轉重復序列程序

EF565211.FASTA - 果蠅性別相關基因FASTA格式序列

EF565211.EINV - 輸出結果倒轉重復文件

-outseq EF565211_EINV.FASTA - 輸出結果FASTA格式文件

運行結果表明，該序列中1 617-1 966/2 355-2 699位為倒轉重復，中間有1個長度為6的空位。查看該序列注釋信息，該倒轉重復序列與果蠅性別比例抑制功能有關。

5 蛋白質序列分析

5.1 序列組分統計

EMBOSS中用于蛋白質一級結構氨基酸序列統計分析的程序包括pepstats, pepinfor, wordcount和 compseq等。pepstats以文本和表格方式輸出蛋白質序列中各種氨基酸含量，并對不同類型氨基酸進行統計，如親水氨基酸和帶電氨基酸等；pepinfor則以圖形方式顯示各種類別氨基酸在序列不同區域的分布，如疏水性氨基酸、極性氨基酸、帶電氨基酸和芳香族氨基酸等。此外，用于核酸序列組分分析的wordcount和compseq也可用于蛋白質序列組分分析。

5.1.1 pepstats用法實例

以水稻落粒控制基因sh4蛋白質產物(UniProt序列條目Q1PIH9_ORYSI)為例，運行程序pepstats，則可統計20種氨基酸出現頻率。

pepstats Q1PIH9_ORYSI.FASTA Q1PIH9_ORYSI.PEPSTATS

pepstats - 統計蛋白質序列中不同氨基酸出現頻率程序

Q1PIH9_ORYSI.FASTA - 水稻落粒控制基因蛋白質FASTA格式序列

Q1PIH9_ORYSI.PEPSTATS - 水稻落粒控制基因蛋白質氨基酸出現頻率輸出結果文件

運行結果輸出20種氨基酸頻數和百分比。水稻落粒控制基因所編碼蛋白質長度為390個氨基酸殘基，不同氨基酸出現頻率很不均勻，某些氨基酸出現頻率較高，如脯氨酸、丙氨酸、絲氨酸和精氨酸高于平均值，而苯丙氨酸、異亮氨酸和甲硫氨酸則低于平均值。

5.1.2 Wordcount用法實例

從以上分析可以看出，水稻落粒控制基因編碼蛋白質中有些氨基酸出現頻率偏高。利用程序wordcount，指定不同字長，可進一步分析水稻落粒控制基因蛋白質產物短片段重復序列出現頻率。

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 - 統計蛋白質序列中指定長度字串出現頻率程序

Q1PIH9_ORYSI.FASTA - 水稻落粒控制基因所編碼蛋白質FASTA格式序列

Q1PIH9_ORYSI.WORD3 - 三肽片段出現頻率輸出結果文件

Q1PIH9_ORYSI.WORD4 - 四肽片段出現頻率輸出結果文件

Q1PIH9_ORYSI.WORD5 - 五肽片段出現頻率輸出結果文件

結果發現，該蛋白質序列中存在大量短片段重復序列(見表7)。

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

5.2 序列特征位點識別

EMBOSS中用于蛋白質序列特征位點分析的程序包括antigenic、sigcleave和digest等，其中antigenic用于抗原決定簇分析，sigcleave用于信號肽剪切位點分析，digest用于酶切位點分析。

5.2.1 antigenic用法實例

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

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

antigenic - 抗原決定簇預測程序

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

-sbegin 35 - 預測起始位點35

-send 142 - 預測終止位點142

CEAM3_HUMAN.ANTI - 輸出結果文件

-minlen 10 - 抗原決定簇最小序列長度10(默認值為6)

預測結果表明，該蛋白質分子可能有三個抗原決定簇，分別位于第48-66, 75-86和119-129位。

5.2.2 fuzzprot用法實例

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

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

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

fuzzpro - 序列模體尋找程序

17ARATH_SPLD.FASTA - 擬南芥17個SBP家族轉錄因子DNA結合結構域序列

17ARATH_SPLD.FUZ - 輸出結果文件

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

-pname "NLS" - 序列模體名稱

序列模體由用戶指定，保守氨基酸用大寫字符表示，中括號內為可選氨基酸，x為任意氨基酸，括號中為任意氨基酸個數，此處為6(輸入括號時需要加轉義符反斜杠，否則無法正常運行)。

圖6 擬南芥17個植物特異轉錄因子SBP家族DNA結合結構域序列圖標Fig.6 Sequence logo of DNA binding domain in 17 SBP plant-specific transcription factors

5.3 二級結構分析

EMBOSS中用于蛋白質二級結構分析的程序包括tmap, topo, pepwheel, helixturnhelix, pepcoil和garnier等，其中tmap用于跨膜螺旋預測，topo用于顯示跨膜螺旋拓撲結構，helixturnhelix用于預測螺旋-轉角-螺旋構象，pepcoil用于預測無規卷曲肽段，pepwheel用于顯示alpha螺旋輪，garnier用于二級結構預測。

5.3.1 tmap用法實例

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

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

map - 跨膜螺旋預測程序

PPF1_PEA.FASTA - 豌豆內膜蛋白FASTA格式序列

PPF1_PEA.TMAP - 預測結果輸出文件

-goutfile PPF1_PEA_TMAP - 圖形輸出文件

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

5.3.2 pepwheel用法實例

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

pepwheel PPF1_PEA_H1.FASTA -graph svg -goutfile PPF1_PEA_H1_PEPWHEEL

pepwheel - 繪制alpha螺旋輪程序

PPF1_PEA_H1.FASTA - 豌豆內膜蛋白中預測到的第1個跨膜螺旋(序列如下)

>PPF1_PEA_H1 | 111-135

SVHVPYSYGFAIILLTVIVKAATLP

-goutfile PPF1_PEA_H1_PEPWHEEL - 圖形文件名

5.3.3 garnier用法實例

以抹香鯨肌紅蛋白(PDB登錄號1MBN)、人癌胚抗原N-端結構域(PDB登錄號2QSQ)和人泛素蛋白(PDB登錄號1UBQ)為例，從PDB數據庫分別下載這3個蛋白質分子FASTA格式序列，分別用以下命令運行程序garnier，即可預測得到二級結構。

garnier 1MBN.FASTA 1MBN.GARNIER

garnier 2QSQ.FASTA 2QSQ.GARNIER

garnier 1UBQ.FASTA 1UBQ.GARNIER

garnier -二級結構預測程序

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

1MBN.GARNIER - 抹香鯨肌紅蛋白預測結果文件

2QSQ.FASTA - 人癌胚抗原N-端結構域FASTA格式序列

2QSQ.GARNIER - 人癌胚抗原N-端結構域預測結果文件

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

1UBQ.GARNIER -人泛素蛋白預測結果文件

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

6 總結和討論

6.1 使用說明

以上我們介紹了EMBOSS軟件包中一些常用程序。經過十年多開發，EMBOSS已成為核酸和蛋白質序列分析常用軟件包，為廣大生物學工作者廣泛使用。EMBOSS軟件包功能齊全、用途廣泛。選修“實用生物信息技術”課程的同學，在學習EMBOSS軟件包中常用程序后，編寫了以下順口溜。

EMBOSS軟件包，包羅萬象真的好，

核酸蛋白都適用，功能強大效率高。

比對進化設引物，翻譯酶切找重復，

程序名稱不記得，wossname幫你找。

程序命令不會用，tfm 把你教。

參數設置技巧高，點點滴滴要記牢，

EMBOSS是法寶，活學活用不愁了。

要熟練使用EMBOSS軟件包中的程序，首先必須熟悉分子生物學基本概念，如中心法則和序列-結構-功能關系等；掌握必要的分子生物學基礎知識，如基因結構、啟動子、外顯子、內含子、編碼序列、RNA二級結構、蛋白質結構層次、一級結構序列特征和二級結構單元，以及序列模體、結構域、蛋白質家族和蛋白質功能等。

選擇合適的程序、設置恰當的參數，是正確、高效使用EMBOSS軟件包的關鍵。除了深入了解所研究問題的生物學背景外，也應搞清輸入數據的種類、格式，掌握各種參數的含義，對所用程序的算法有所了解。同樣一個問題，使用不同程序，運行結果就可能不同；即使是同一個程序，參數不同，結果也可能不同。熟練使用EMBOSS軟件包提供的三個幫助程序wossname、tfm和seealso，深入理解各個程序的功能和可供設置的參數，可以在程序使用過程中起到事半功倍的效果。

需要說明的是，EMBOSS軟件包啟動時，人類基因組計劃尚未完成，基因組數據分析剛剛開始。因此，EMBOSS不是組學數據分析軟件，而是針對單個或多個蛋白質或核酸序列分析工具，其功能相當于基于個人計算機的共享軟件BioEdit或商業軟件DNAStar和MacVector等。從事基因組和轉錄組等高通量數據分析的讀者，可選擇Bowtie, BWA, TopHat/Cufflinks等軟件。此外，EMBOSS軟件包是單個程序的集成，各個程序之間并無聯系，而后來開發的Galaxy平臺，則將某些工具整合而形成互相關聯的分析流程。

與所有計算機軟件均可能存在“bug”一樣，EMBOSS軟件包中某些程序在運行時結果可能有誤。例如，點陣圖程序dotmatcher和dottup在比較兩個不同序列時，坐標軸顯示有誤。此外，由于近年來UniProt數據庫格式有所調整，序列提取程序extractfeat在處理UniProt/Swiss-Prot格式輸入文件時，得不到正確結果。讀者可自行修改源代碼改正錯誤，或與Peter Rice聯系。

6.2 程序列表

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

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

pepstats一級結構分析統計蛋白質序列中不同種氨基酸出現頻率pepinfo一級結構分析用圖形方式顯示蛋白質序列不同氨基酸分布特征antigenic一級結構分析預測蛋白質序列中抗原決定簇sigcleave一級結構分析尋找蛋白質序列中信號肽切割位點digest一級結構分析尋找蛋白質序列中蛋白酶酶切位點fuzzprot一級結構分析尋找蛋白質序列中序列模體Tmap二級結構分析預測蛋白質序列中的跨膜螺旋topo二級結構分析顯示跨膜螺旋拓撲結構pepwheel二級結構分析繪制alpha螺旋輪pepcoil二級結構分析預測無規卷曲區域helixturn-helix二級結構分析螺旋-轉角-螺旋序列模體分析garnier二級結構分析蛋白質序列二級結構預測

EMBOSS網站列出了所有程序的名稱和用途，也可用wossname命令按功能分類或字母表順序列出所有程序。

wossname -search“”

按功能分類列出所有程序

wossname -search“”-alphabetic

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

除了EMBOSS開發團隊自行編寫的程序外，EMBOSS還整合了不少其它常用生物信息軟件包，如基于隱馬爾可夫模型的蛋白質結構域序列譜構建和結構域識別軟件包HEMMER、系統發育分析軟件包Phylip及RNA二級結構分析和預測軟件包VIENNA等。限于篇幅，本文未介紹這些軟件包中程序的用法。

6.3 回顧和展望

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

本世紀初，Peter Rice領導的EMBOSS研發團隊受聘于歐洲生物信息學研究所，完成了該軟件包的主要開發。2009年，EMBOSS項目曾得到英國生物技術和生命科學研究委員會(Biotechnology and Biological Science Research Council, BBSRC)資助，繼續進行開發。目前，該軟件包開發項目已經結束，由英國transSMART基金會Peter Rice負責維護。顯然，作為開源軟件，EMBOSS的進一步開發，需要得到生物信息軟件開發人員和廣大用戶的支持。對該軟件包開發感興趣的讀者可與Peter Rice直接聯系，聯系方式請參閱補充材料1。

致 謝

感謝楊德昌安裝和調試EMBOSS軟件包，顏林林改正點陣圖程序中的錯誤。感謝樊麗編寫的EMBOSS順口溜。金錄佳、李宏博和趙坤認真閱讀并校正了初稿中多處文字錯誤。感謝匿名審稿人寶貴的修改意見。感謝中國科學院北京基因組研究所(國家生物信息中心)對EMBOSS網絡教程提供的支持。感謝北京大學生命科學學院、中國農業科學院研究生院和中國科學院大學生命科學學院多年來對“實用生物信息技術”課程的支持。