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

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

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

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

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

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

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

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

1 SPIDER軟件包的主要函數(shù)及功能

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

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

2 BarcodingR軟件包的主要函數(shù)及功能

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)為基礎(chǔ)建立。該軟件包主要用途是進行物種識別和條形碼分析。具體函數(shù)及功能詳見表2，其中3個較為關(guān)鍵的物種識別函數(shù)依次是‘barcoding.spe.identify’，‘barcoding.spe.identify2’ 和‘bbsik’。

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

2.1 人工智能方法bpNewTraining/bpNewTrainingOnly/bpUseTrained參數(shù)

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

2.2 模糊集合理論fuzzyId參數(shù)

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

2.3 貝葉斯方法Bayesian參數(shù)

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

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

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

3 應(yīng)用實例分析

3.1 數(shù)據(jù)輸入

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

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

>data(dolomedes) #加載dolomedes數(shù)據(jù)

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

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

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

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

>ref.dna <- read.dna("ref.fas",format="fasta") #讀入?yún)⒖夹蛄?/p>

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

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

3.2 從GenBank和BOLD數(shù)據(jù)庫中下載數(shù)據(jù)

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

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

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

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

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

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

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

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

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

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

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

>data(anoteropsis) #加載anoteropsis數(shù)據(jù)

>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%閾值范圍內(nèi)進行物種識別效率計算

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

3.4 物種識別分析

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

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

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

>refdata<- TibetanMoth #獲取 TibetanMoth數(shù)據(jù)集并賦值給變量refdata

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

>ref<- output$ref.selected #定義參考數(shù)據(jù)集

>que<- out$ref.left #定義問詢數(shù)據(jù)集

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

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

#輸出結(jié)果BPidentified.txt詳見附件1

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

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

#輸出結(jié)果fuzzyIDidentified.txt詳見附件2

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

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

#輸出結(jié)果Bayidentified.txt詳見附件3

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

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

>ref<-TibetanMoth #獲取 TibetanMoth數(shù)據(jù)集并賦值給變量ref

>set.seed(10) #產(chǎn)生隨機數(shù)

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

>ref2<-out$ref.selected #定義參考數(shù)據(jù)集

>que<-out$ref.left #定義問詢數(shù)據(jù)集

>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 #將人工智能物種識別結(jié)果賦值

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

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

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

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

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

續(xù)上表

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

3.6 非編碼基因的物種識別

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

>refdata<- pineMothITS1 #獲取pineMothITS1數(shù)據(jù)集并賦值給變量refdata

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

>ref<- out$ref.selected #定義參考數(shù)據(jù)集

>que<- out$ref.left #定義問詢數(shù)據(jù)集

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

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

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

#輸出結(jié)果ITS1identified.txt詳見附件4

3.7 DNA條形碼空白區(qū)

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

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

#形碼空白區(qū)圖

圖2 DNA條形碼空白區(qū)Fig.2 DNA barcoding gap

4 結(jié)論

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

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].環(huán)境昆蟲學報，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