金 倩，阮倩倩，陳 芬，羅桂杰，張愛兵

DNA條形碼專用R包及其主要功能簡介

金 倩1,2，阮倩倩1，陳 芬1，羅桂杰1，張愛兵2*

(1.江蘇省農業科學院宿遷農科所，江蘇宿遷 223800；2.首都師范大學生命科學學院，北京 100048)

SPIDER和BarcodingR軟件包作為DNA條形碼研究的專用R軟件包，整合了基于DNA條形碼的物種識別算法和數據分析策略，避免了各種非商業用途軟件的局限性，為DNA條形碼的研究提供了便捷性和可操作性。本文對SPIDER和BarcodingR軟件包的主要功能進行整體介紹，并以西藏夜蛾科COI、松毛蟲屬ITS等數據為例，進行主要函數的演示分析，希望為昆蟲DNA條碼的應用分析提供幫助。

DNA條形碼；R語言；SPIDER；BarcodingR

近年來，伴隨著DNA條形碼的廣泛應用和后基因組時代的到來，生物多樣性的研究變得更加深入和便捷。截至2016年12月30日，BOLD系統中共記錄物種261170個，條形碼序列5281492條(http://www.barcodinglife.org)。在DNA條形碼數據快速增長的背景下，物種識別和構建類群全面的參考數據庫成為生物多樣性研究的核心任務，進而推動了物種識別算法等理論研究的不斷發展。基于DNA條形碼的識別方法及有效性，包括經典的基于距離、基于進化樹和基于特征的方法(Sarkaretal., 2002；Blaxteretal., 2005；Toffolietal., 2008)。近年來研究者也試圖提出各種新的方法，如，人工智能方法(Zhangetal., 2012a)、模糊集合理論(Zhangetal., 2012b)、構建識別效率與距離閾值之間模型的方法(Virgilioetal., 2012)、基于貝葉斯理論和形態有限介入的方法(Jinetal., 2013)等。

上述理論方法的實現對用戶的生物信息學背景要求較高，有些軟件因缺乏友好的用戶界面或涉及非商業用途而被限制使用，對試圖通過DNA條形碼技術研究生物類群的實例研究者而言，復雜的算法和程序往往使他們望而生畏。SPIDER軟件包(Brownetal., 2012)和BarcodingR軟件包(Zhangetal., 2016)的相繼提出避免了以上非商業用途軟件的局限性，發揮R語言本身在數據處理、數據統計和數據可視化方面較強的優勢，為DNA條形碼的研究提供了便捷性和可操作性。本文對SPIDER和BarcodingR軟件包的主要功能進行系統的介紹，并以西藏夜蛾科COI(Jinetal., 2013)、松毛蟲屬ITS(Daietal., 2013)等數據為例，進行主要函數的演示分析。

1 SPIDER軟件包的主要函數及功能

SPIDER軟件包以軟件包“APE”(Paradisetal., 2004)，“PHYCLUST”(Chen and Dorman, 2010)，“PEGAS”(Paradis, 2010)和“ADEGENET”(Jombart, 2008)為基礎建立，涉及到條形碼分析的主要函數及功能詳見表1。該軟件包首先可以對序列的屬性進行統計，例如可以根據dataStat函數計算數據的屬、種及樣本數量；其次，能夠實現經典的基于距離樹的物種識別效率計算，需要注意的是數據集中的所有序列必須事先進行識別鑒定，并建立參考條形碼數據庫。在計算識別效率時，每一條序列被作為未知序列問詢數據庫中參考序列，如果識別的結果與已知結果一致，表明識別成功，否則為失敗；最后，該軟件包還可以檢驗系統發育樹的單系性，能夠在Genbank和BOLD數據庫中直接獲取數據，實現滑動窗口分析等功能。

表1 SPIDER軟件包主要函數及功能Table 1 List of main functions and their descriptions of the ‘SPIDER’ package

2 BarcodingR軟件包的主要函數及功能

BarcodingR軟件包以軟件包“APE”(Paradisetal., 2004)，“ADEGENET”(Jombart, 2008)，“MCLUST”(Fraley and Raftery, 1999)，“PHYCLUST”(Chen and Dorman, 2010)，“CCLUST” (https://cran.r-project.org/)，“PEGAS”(Paradis, 2010)，“SEQINR”(Charifetal., 2005)，“NNET” (Ripley, 1996)，“E1071”(https://cran.r-project.org/)和“SPIDER”(Brownetal., 2012)為基礎建立。該軟件包主要用途是進行物種識別和條形碼分析。具體函數及功能詳見表2，其中3個較為關鍵的物種識別函數依次是‘barcoding.spe.identify’，‘barcoding.spe.identify2’ 和‘bbsik’。

barcoding.spe.identify函數進行物種識別主要通過人工智能方法、模糊集合理論和貝葉斯方法，通過以下參數進行控制：

2.1 人工智能方法bpNewTraining/bpNewTrainingOnly/bpUseTrained參數

當數據集較小時(<500條序列)，可以選擇bpNewTraining參數，訓練集和問詢集同時計算。當數據集較大時(>500條序列)需要通過bpNewTrainingOnly參數先進行人工智能訓練，然后再通過bpUseTrained進行問詢集的物種識別。訓練的模型參數被存放于一個臨時文件中，當運行bpUseTrained時加載訓練結果并進行物種識別。

2.2 模糊集合理論fuzzyId參數

每一條問詢序列的識別結果都被賦予一個FMF值(范圍從0～1)，代表序列分配的可能性(Zhangetal., 2012a)。為了提高識別效率，在搜索參考數據過程中采用KNN搜索算法替代MD算法(http://www.stats.ox.ac.uk/pub/MASS4/)。

2.3 貝葉斯方法Bayesian參數

按照金倩等(Jinetal., 2013)提出的基于貝葉斯理論和形態有限介入的方法對蛾類進行物種識別，主要通過Bayesian函數實現。

表2 BarcodingR軟件包主要函數及功能Table 2 List of main functions and their descriptions of the ‘BarcodingR’ package

除了COI序列之外，非編碼基因ITS1和ITS2也被廣泛應用于動植物類群中。但因其間隔區的存在影響了多序列比對的準確性(Zhangetal., 2008；Zhangetal., 2012b)，故張愛兵等(Zhangetal., 2012b)提出在非編碼基因中可以免除序列比對過程直接進行物種識別，但該方法受限于MATLAB環境下使用。BarcodingR軟件包規避了語言環境的局限性，通過kmer方法實現非編碼基因的序列比對和識別。

3 應用實例分析

3.1 數據輸入

SPIDER軟件包可以通過ape軟件包中read.dna函數read.nexus.data函數讀入分子數據，通過R基礎包中的read.table讀入形態數據，也可以加載軟件包本身的數據。以讀入軟件包中數據和用戶自己的數據為例，R語言的命令代碼如下，#后面的文字表示對R語言命令的注釋。

>library(spider) #加載spider軟件包

>data(dolomedes) #加載dolomedes數據

>setwd("c:\My documents\R") #設定數據集所在路徑

>dataset <- read.dna("mySequences.fas",format="fasta") #讀入"mySequences.fas"數據

BarcodingR軟件包的輸入較復雜，需提前準備好參考序列和問詢序列，格式為常用的.fasta或.fas格式文件。參考序列的.fasta文件大于號之后必須含有分類信息“>seqID, species_names”，例如：“>seq1, Noctuidae_Blepharosis_paspa”，“>seq2,Apamea_devastator”。問詢序列為一般的.fasta格式，例如“>seqID”。通過ape軟件包中read.dna函數或adegent軟件包中fasta2DNAbin函數讀入數據(Paradisetal., 2004)。

>library(ape) #加載ape軟件包

>ref.dna <- read.dna("ref.fas",format="fasta") #讀入參考序列

>library(adegenet) #加載adegenet軟件包

>que.dna <- fasta2DNAbin("que.fas") #讀入問詢序列

3.2 從GenBank和BOLD數據庫中下載數據

APE軟件包中的read.GenBank函數能夠直接下載GenBank中的序列，并存儲為DNAbin格式的對象中，但該函數不能獲得物種名字或者基因區域等信息，因此SPIDER軟件包中的read.GB函數對其進行了改進，以GenBank的登錄號為向量，并檢索記錄對應的序列。以西藏夜蛾數據為例(Jinetal., 2013)，在GenBank中下載序列步驟如下：

>seq <- 392408:392727 #將392408至392727編號賦值給seq

>seq <- paste("JX", seq, sep="") #在每個編號前添加"JX"

>TibetData <- read.GB(seq) #下載數據并賦值給TibetData

SPIDER軟件包中的search.BOLD函數在BOLD數據庫中搜尋某類群的條形碼數據，并返回搜尋對象的樣本編號，read.BOLD函數通過這些樣本編號下載對應的序列。以查找夜蛾科物種Himalaeaunica為例(Jinetal., 2013)，R語言的命令代碼如下：

>nums <- search.BOLD("Himalaeaunica") #在BOLD數據庫中查找夜蛾科Himalaea

#unica類群并將樣本編號返回給nums

>NoctuidaeSp <- read.BOLD(nums) #下載mums對應的序列并賦值給NoctuidaeSp

3.3 基于距離的物種識別效率分析

以SPIDER軟件包中Anoteropsis數據為例(Vink and Paterson, 2003)，進行threshID函數計算，R語言的命令如下，結果如圖1所示，“correct”代表在1%的閾值范圍內問詢序列識別正確的是11個，“incorrect”代表在1%的閾值范圍內問詢序列識別錯誤的是2個，“no id”代表在1%的閾值范圍內問詢序列沒有匹配到任何一個個體是20個。

>library(spider) #加載spider軟件包

>data(anoteropsis) #加載anoteropsis數據

>aa <- strsplit(dimnames(anoteropsis)[[1]], split="_") #將序列名字按照“_”分開

>anoSpp <- sapply(aa, function(x) paste(x[1], x[2], sep="_")) #將屬名種名提取出來

>inputDist<-dist.dna(anoteropsis,pairwise.deletion=TRUE) #計算K2P距離

>table(threshID(inputDist,anoSpp)) #1%閾值范圍內進行物種識別效率計算

圖1 基于1%閾值的物種識別結果Fig.1 Species identification results with 1% threshold

3.4 物種識別分析

以BarcodingR軟件包中的夜蛾科數據為例(Jinetal., 2013)。首先需要加載軟件包和數據，然后通過控制barcoding.spe.identify函數中的method參數控制不同的物種識別算法，R語言的命令如下：

>install.packages("BarcodingR") #安裝BarcodingR軟件包

>library(BarcodingR) #加載BarcodingR軟件包

>refdata<- TibetanMoth #獲取 TibetanMoth數據集并賦值給變量refdata

>output<- sample.ref(refdata, sample.porp = 0.9,sample.level = "species") #隨機抽樣

>ref<- output$ref.selected #定義參考數據集

>que<- out$ref.left #定義問詢數據集

>bp<- barcoding.spe.identify(ref,que,method ="bpNewTraining") #人工智能物種識別

>save.ids(outfile="BPidentified.txt",bp) #輸出結果

#輸出結果BPidentified.txt詳見附件1

>fuzzyID<- barcoding.spe.identify(ref,que,method ="fuzzyId") #模糊集合物種識別

>save.ids(outfile="fuzzyIDidentified.txt",fuzzyID) #輸出結果

#輸出結果fuzzyIDidentified.txt詳見附件2

>Bay<- barcoding.spe.identify(ref,que,method ="Bayesian") #貝葉斯物種識別

>save.ids(outfile="Bayidentified.txt",Bay) #輸出結果

#輸出結果Bayidentified.txt詳見附件3

3.5 不同物種識別方法結果的一致性檢驗及匯總

通過BarcodingR軟件包中的函數consensus.identify，實現人工智能方法、模糊集合理論和貝葉斯方法的一致性檢驗，以西藏夜蛾數據為例(Jinetal., 2013)，不同方法的物種一致性檢驗R命令如下，結果詳見表3。

>ref<-TibetanMoth #獲取 TibetanMoth數據集并賦值給變量ref

>set.seed(10) #產生隨機數

>out<-sample.ref(ref,sample.porp=0.95,sample.level="species") #隨機抽樣

>ref2<-out$ref.selected #定義參考數據集

>que<-out$ref.left #定義問詢數據集

>bsi0<-barcoding.spe.identify(ref2, que, method ="bpNewTrainingOnly") #人工智能物種

#識別

>bsi1<-barcoding.spe.identify(ref2, que, method ="bpUseTrained") #人工智能物種識別

>bsi2<-barcoding.spe.identify(ref2, que, method ="fuzzyId") #模糊集合物種識別

>bsi3<-barcoding.spe.identify(ref2, que, method ="Bayesian") #貝葉斯物種識別

>que.IDs<-gsub(",.+","",rownames(que)) #提取序列ID

>bpid<-bsi1$output_identified$spe.Identified #將人工智能物種識別結果賦值

>fuzzyid<-bsi2$output_identified$spe.Identified #將模糊集合物種識別結果賦值

>Bayesianid<-bsi3$output_identified$output_identified #將貝葉斯物種識別結果賦值

>identifications<-data.frame(queIDs = que.IDs, pid = bpid, fuzzyid = fuzzyid, Bayesianid = Bayesianid) #將三種方法的物種識別結果輸出放入數據框identifications中

>ccs<-consensus.identify(identifications) #一致性檢驗

表3 基于不同的物種識別方法的物種一致性Table 3 List of main functions and their descriptions of the ‘BarcodingR’ package

續上表

編號Code問詢序列IDQueIDs物種識別一致性結果Consensus.id投票數Votes4ML0829139TNoctuidaePoliamortua35RKZ0908014TNoctuidaeAnartatrifolii36ML0829120MNoctuidaePhlogophorasubpurpurea37RKZ0908008TNoctuidaeAtracheaspproblematic38ML0829004MNoctuidaeDiarsiasp3problematic39LS0909017MNoctuidaePerissandriasheljuzhkoiAcuta310JZ0907051MNoctuidaePerissandriaficta311ML0829149TNoctuidaePerissandriasikkima312LS0909022MNoctuidaeXestiacnigrum313LS0909041RNoctuidaeXestiacnigrum314BM0830061TNoctuidaeXestiacnigrum315LS0909044MNoctuidaeAthethisspn316LS0909045TNoctuidaeAthethisspn317BM0901038RNoctuidaeHermonassasp318BSC0902032MNoctuidaeHermonassaspn319ML0829144MNoctuidaeAgrotisjusta320JZ0907026MNoctuidaeAgrotisscotacra321BM0830066TNoctuidaeAgrotismacrobscura322RKZ0908002MNoctuidaeEuxoakotzschi323JZ0907038TNoctuidaeEuxoasp1problematic324SN0906006TNoctuidaeDichagyrisastigmata325SN0906015TNoctuidaeDichagyrisastigmata326ML0829147MNoctuidaeunknown1sp1problematic327JZ0907044MNoctuidaeAtracheaparvispina328JZ0907047TNoctuidaeApamealateritiaObfuscate329JZ0907027TNoctuidaeBlepharosisspproblematic330JZ0907032TNoctuidaeBlepharosispaspa331LS0909019TNoctuidaeMacdunnoughiacrassisigna332JZ0907007TNoctuidaeMacdunnoughiacrassisigna333BM0830013TNoctuidaeEphesiabutleri334BM0830016TNoctuidaeEphesiabutleri335SN0906001TNoctuidaeEphesiafulminea336BM0901016TNoctuidaeCatocalahyperconnexa337JZ0907031BNoctuidaeBlepharosispaspa238ML0829131TNoctuidaeApameadevastator3

3.6 非編碼基因的物種識別

以BarcodingR軟件包中的非編碼基因數據ITS1為例(Daietal., 2013)，通過函數barcoding.spe.identify2實現，R語言的命令如下：

>refdata<- pineMothITS1 #獲取pineMothITS1數據集并賦值給變量refdata

>out<- sample.ref(refdata,sample.porp = 0.9,sample.level="species") #隨機抽樣

>ref<- out$ref.selected #定義參考數據集

>que<- out$ref.left #定義問詢數據集

>ITS1identified<- barcoding.spe.identify2(ref,que,kmer=10,optimization="T") #基于kmer

#的模糊集合理論物種識別

>save.ids(outfile="ITS1identified.txt",ITS1identified) #輸出結果

#輸出結果ITS1identified.txt詳見附件4

3.7 DNA條形碼空白區

通過函數barcoding.gap實現，該函數輸出結果為DNA條形碼空白區(如圖2)，同時還能顯示種內和種間距離的最大值、最小值、中位數和平均值。R語言的命令如下:

>b.gap <- barcoding.gap(ref = TibetanMoth, dist = “K80”) #計算種內種間距離并繪制條

#形碼空白區圖

圖2 DNA條形碼空白區Fig.2 DNA barcoding gap

4 結論

SPIDER軟件包(Brownetal., 2012)可以進行常規的條形碼統計分析，能夠實現經典的基于距離的物種識別效率計算。BarcodingR軟件包(Zhangetal., 2016)是SPIDER的后續補充，將近年來新提出的物種識別方法進行了整合，包括人工智能方法、模糊集合理論、基于貝葉斯理論和形態有限介入的方法。兩個軟件包開發利用R語言平臺的免費性、開源性和可擴展性為DNA條形碼的研究提供了便捷。SPIDER軟件包和BarcodingR軟件包都可以從網站中直接下載并使用，SPIDER軟件包的網址為http://spider.r-forge.r-project.org/，BarcodingR軟件包的網址為：https://cran.r-project.org/web/packages/BarcodingR/index.html或者https://github.com/zhangab2008/BarcodingR.git。二者都適用于Microsoft Windows，MAC OSX 和LINUX系統。

References)

Brown SDJ, Collins RA, Boyer S,etal.SPIDER: An R package for the analysis of species identity and evolution, with particular reference to DNA barcoding [J].MolecularEcologyResources, 2012, 12 (3): 562-565.

Charif D, Thioulouse J, Lobry JR,etal.Online synonymous codon usage analyses with the ade4 and seqinR packages [J].Bioinformatics, 2005, 21 (4): 545-547.

Dai QY, Gao Q, Wu CS,etal.Phylogenetic reconstruction and DNA barcoding for closely related pine moth species (Dendrolimus) in China with multiple gene markers [J].PLoSONE, 2012, 7 (4): e32544.

Fraley C, Raftery AE.MCLUST: Software for model-based cluster analysis [J].JournalofClassification, 1999, 16 (2): 297-306.

Jin Q, Han HL, Hu XM,etal.Quantifying species diversity with a DNA barcoding-based method: Tibetan moth species (Noctuidae) on the Qinghai-Tibetan Plateau [J].PLoSONE, 2013, 8 (5): e64428.

Jombart T.adegenet: A R package for the multivariate analysis of genetic markers [J].Bioinformatics, 2008, 24 (11): 1403-1405.

Paradis E, Claude J, Strimmer K.APE: Analyses of phylogenetics and evolution in R language [J].Bioinformatics, 2014, 20 (2): 289-290.

Paradis E.pegas: An R package for population genetics with an integrated-modular approach [J].Bioinformatics, 2010, 26 (3): 419-420.

Ripley BD.Pattern Recognition and Neural Networks [M].New York: Cambridge University Press, 1996.

Vink CJ, Paterson AM.Combined molecular and morphological phylogenetic analyses of the New Zealand wolf spider genusAnoteropsis(Araneae: Lycosidae) [J].MolecularPhylogeneticsandEvolution, 2003, 28 (3): 576-587.

Zhang AB, Feng J, Ward RD,etal.A new method for species identification via protein-coding and noncoding DNA barcodes by combining machine learning with bioinformatics methods [J].PLoSONE, 2012b, 7 (2): e30986.

Zhang AB, Hao MD, Yang CQ,etal.BarcodingR: An integrated R package for species identification using DNA barcodes [J].MethodsinEcologyandEvolution, 2016: 1-8.

Zhang AB, Muster C, Liang HB,etal.A fuzzy-set-theory-based approach to analyse species membership in DNA barcoding [J].MolecularEcology, 2012a, 21 (8): 1848-1863.

Zhang AB, Sikes DS, Muster C,etal.Inferring species membership using DNA sequences with back-propagation neural networks [J].SystematicBiology, 2008, 57 (2): 202-215.

Main functions and descriptions of R packages used for DNA barcoding

JIN Qian1,2, RUAN Qian-Qian1, CHEN Fen1, LUO Gui-Jie1, ZHANG Ai-Bing2*

(1.Suqian Institute of Agricultural Sciences, Jiangsu Academy of Agricultural Sciences, Suqian 223800, Jiangsu Province, China; 2.College of Life Sciences, Capital Normal University, Beijing 100048, China)

SPIDERand BarcodingR, are new R packages implementing a number of useful functions for DNA barcoding analyses and associated research into species speciation.Both of them avoid limitations of noncommercial computer programs and provide convenience and operability for user.Here we demonstrated the main functions of the two packages with Tibet moth, pine moth datasets.

DNA barcoding; R language; SPIDER; BarcodingR

特邀稿件InvitedReview

金倩，阮倩倩，陳芬,等.DNA條形碼專用R包及其主要功能簡介[J].環境昆蟲學報，2017,39(3):485-492.

國家自然科學基金青年科學基金項目(31601877)

金倩，女，1987年生，博士研究生，研究方向為物種界定及昆蟲多樣性，E-mail: jinhongyu2001@163.com

*通訊作者Author for correspondence，E-mail: zhangab2008@mail.cnu.edu.cn

Received: 2017-03-01；接受日期Accepted: 2017-04-27

Q963

A

1674-0858(2017)03-0485-08