孫小潔，鄭方強，曾健明

(1.山東農(nóng)業(yè)大學(xué) 植物保護(hù)學(xué)院，山東 泰安 271018；2.珠海健明生物醫(yī)藥科技有限公司，廣東 珠海 519000)

基因芯片技術(shù)自20世紀(jì)80年代發(fā)展至今已產(chǎn)生了大量的基因表達(dá)數(shù)據(jù)。如何從復(fù)雜的基因大數(shù)據(jù)中進(jìn)行知識發(fā)現(xiàn)，是生物信息學(xué)研究的重要課題之一。為了滿足對高通量基因表達(dá)數(shù)據(jù)存儲不斷增長的需求，美國國家生物技術(shù)信息中心(NCBI)建立了基因表達(dá)數(shù)據(jù)庫(GEO)[1-4]，為用戶提供了可供數(shù)據(jù)提交、存儲和檢索的平臺。目前，GEO數(shù)據(jù)庫已經(jīng)收錄了累計10萬多個系列、280多萬個樣本的數(shù)據(jù)，涉及3 000多種生物[5]。

面對海量復(fù)雜的生物數(shù)據(jù)，研究者的思維方式也相應(yīng)地從數(shù)據(jù)的生成轉(zhuǎn)向?qū)?shù)據(jù)的深入挖掘和分析。數(shù)據(jù)挖掘是從大量的數(shù)據(jù)中通過算法搜索隱藏于其中信息的過程[6]。將數(shù)據(jù)挖掘方法應(yīng)用于生物信息大數(shù)據(jù)，能夠從中挖掘出有價值的信息，尋找潛在規(guī)律，進(jìn)而對相關(guān)疾病機(jī)制作出科學(xué)的詮釋，是當(dāng)前生物信息學(xué)的熱點問題之一。

基因表達(dá)芯片是采用傳統(tǒng)的基因表達(dá)量測定方法，會產(chǎn)生出大量有價值的數(shù)據(jù)，是生物信息數(shù)據(jù)挖掘工作的重要組成部分?；虮磉_(dá)芯片測序的結(jié)果是每個樣品的探針表達(dá)量，在后續(xù)分析過程中需要根據(jù)基因與探針之間的對應(yīng)關(guān)系進(jìn)行ID轉(zhuǎn)換，進(jìn)而計算基因的表達(dá)量高低。部分芯片平臺可以從Bioconductor網(wǎng)站的注釋程序包中直接獲取這種對應(yīng)關(guān)系，但只覆蓋了約90個常用的芯片，而現(xiàn)存的測序平臺有10 000多個，且日益增長；也有一些芯片平臺可以從生產(chǎn)廠家的官方網(wǎng)站或GEO數(shù)據(jù)庫的通用公共許可證(GPL)平臺信息表格中查找；更多芯片平臺則是僅提供了探針I(yè)D與序列信息，而未提供現(xiàn)成的探針與基因的對應(yīng)關(guān)系[7-8]。

準(zhǔn)確的探針注釋是芯片數(shù)據(jù)下游分析的前提，確保能對分析結(jié)果進(jìn)行正確的生物學(xué)解釋。目前的注釋存在兩個主要問題：其一是基因ID沒有一個統(tǒng)一的標(biāo)準(zhǔn)，每個數(shù)據(jù)庫都使用其特定的基因ID，主流的有Official_Gene_ID、NCBI的Entrez_Gene_ID、Genebank GI號、Gene Accession、RefSeq_ accession、Ensembl_Gene_ID等；此外還有Vaga gene ID、havana_gene_ID、ena等[9-10]?；騃D的復(fù)雜多樣，導(dǎo)致已有的芯片注釋依據(jù)的基因ID也不統(tǒng)一；另外，芯片注釋是根據(jù)以往的參考基因組設(shè)計和比對的，而參考基因組的版本多樣，且時常更新。參考基因組存儲于Ensembl[11]、UCSC Genome Browser[12-13]以及NCBI 3個數(shù)據(jù)庫，每個數(shù)據(jù)庫中都存放了多個參考基因組版本。不同的基因芯片注釋依據(jù)的參考基因組版本不統(tǒng)一，更新速度較慢，有些甚至不更新。

基因芯片注釋過時，ID不統(tǒng)一的混亂現(xiàn)狀，使存放在GEO數(shù)據(jù)庫中大量有價值的數(shù)據(jù)無法利用起來，給芯片數(shù)據(jù)挖掘工作帶來了較大的困難，如果直接使用過時的注釋文件，勢必導(dǎo)致后續(xù)分析結(jié)果與最新的基因注釋大相徑庭。因此，以最新的基因組為參考，對探針序列進(jìn)行重新注釋，是芯片數(shù)據(jù)分析過程中至關(guān)重要的工作。Yin等[14]整合了多個數(shù)據(jù)庫中的斑馬魚基因注釋，將Affymetrix公司的斑馬魚基因表達(dá)芯片探針序列映射到整合的轉(zhuǎn)錄本中，大幅增加了檢測到的基因數(shù)量、差異基因和可變剪切數(shù)量。同年，Barbosa-Morais等[15]發(fā)現(xiàn)Illumina公司提供的許多芯片原始注釋并不可靠，并針對BeadArrays系列芯片開發(fā)了基于Perl語言的寡核苷酸芯片技術(shù)的重新注釋工具(ReMOAT)；Arloth等[16]也開發(fā)了Illumina芯片重注釋的Perl工具，使用該工具注釋的Human-HT12 v4 芯片有約25%的探針注釋與公司提供的原始注釋不同，并與ReMOAT比較發(fā)現(xiàn)能注釋到更多的探針。近年來，多項長鏈非編碼RNA(lncRNA)的差異分析研究都用到了重注釋，例如非小細(xì)胞肺癌亞型的特異性lncRNA及潛在功能分析[17]。

本文搭建了一套簡便靈活的表達(dá)芯片通用自主注釋流程，以期可以對已有注釋的經(jīng)典芯片平臺進(jìn)行重注釋，并致力于應(yīng)用在無注釋但提供探針序列信息的任一表達(dá)芯片平臺上。

1 系統(tǒng)與方法

1.1 開發(fā)環(huán)境

硬件環(huán)境：云服務(wù)器，16核心，32G內(nèi)存，硬盤1T；操作系統(tǒng)：Ubuntu 16.04.5。

1.2 R軟件及主要程序包

R軟件版本為3.5.2，可從https://mirrors.tuna.tsinghua.edu.cn/CRAN/bin/獲取。

R程序包Rsubread、Rsamtools[18]、refGenome和GenomicRanges，可從http://www.bioconductor.org/獲取，也可在R語言界面使用BiocManager::install()命令安裝。

1.3 數(shù)據(jù)準(zhǔn)備

流程的輸入文件是芯片探針序列文件，通?？梢栽贕EO數(shù)據(jù)庫或芯片廠家官方網(wǎng)站下載探針平臺信息表格，刪除掉多余信息，只留下2列。第一列是探針id(Probe_id)，第二列是探針序列(Sequence)，數(shù)據(jù)結(jié)構(gòu)見表1。

表1 探針序列文件格式Table 1 File formats of probe sequence

推薦以逗號為分隔符，存為csv格式，命名為“GPLxxx.id2sequence.csv”，存放于工作目錄下。

1.4 參考基因組及注釋文件下載

從Ensembl數(shù)據(jù)庫下載最新的人類參考基因組(Reference Genome)Homo_sapiens.GRCh38.dna.primary_assembly.fa和對應(yīng)版本的基因組注釋(Genome Annotation)文件Homo_sapiens.GRCh38.94.gtf，小鼠參考基因組Mus_musculus.GRCm38.dna.primary_assembly.fa和對應(yīng)版本的基因組注釋Mus_musculus.GRCm38.95.gtf，存放于同一目錄下。使用本流程需輸入?yún)⒖蓟蚪M和注釋文件的存放路徑。

1.5 表達(dá)芯片探針自主注釋流程

表達(dá)芯片探針自主注釋流程(圖1)基于R語言，整合了多個R程序包。先讀取芯片和探針的對應(yīng)關(guān)系文件，并將其轉(zhuǎn)換為fasta格式(一種序列存儲格式，是本流程使用的參考基因組序列格式。每條序列的第一行以“>”開頭，跟隨“>”的是序列的ID號及描述信息；第二行開始是序列內(nèi)容；第二條序列另起一行，仍然由“>”開始，以此類推)。將探針序列比對到參考基因組(也稱參考序列，是一個數(shù)字化核酸序列數(shù)據(jù)庫，由科學(xué)家組裝，作為一個物種的一組基因的代表性例子[19-20])，生成BAM格式的比對結(jié)果文件，獲得探針序列在基因組中的位置信息；讀取最新參考基因組的注釋文件，獲得基因序列在基因組中的位置信息。將探針序列與基因序列的位置信息分別轉(zhuǎn)換成Grange對象(即存儲一組基因位置信息的容器，每個基因位置信息由染色體名稱、開始位置、結(jié)束位置和正點鏈來描述)，尋找二者在基因組上的位置重疊區(qū)域，就獲得了基因與探針的對應(yīng)關(guān)系，將其組合為一個數(shù)據(jù)框，導(dǎo)出為csv格式的表格。

圖1 基于R語言的基因表達(dá)芯片注釋流程Fig.1 An R workflow for annotation of geneexpression microarray

根據(jù)參考基因組構(gòu)建索引是序列比對的重要前提，索引僅取決于參考基因組，與需注釋的芯片平臺數(shù)據(jù)無關(guān)，但構(gòu)建索引耗時長、需要較大的內(nèi)存，且會生成約15G的大文件，是限速步驟。流程中對該步驟進(jìn)行了邏輯判斷，同一物種的芯片平臺注釋僅在首次運行時構(gòu)建索引，不會重復(fù)構(gòu)建，后續(xù)進(jìn)行其他芯片平臺注釋時，整個流程可在3 min以內(nèi)迅速完成。其中，基因組注釋為利用生物信息學(xué)方法和工具，對基因組所有基因的生物學(xué)功能進(jìn)行高通量注釋，包括基因識別和基因功能注釋兩個方面，常存為gtf和gff格式[20]；SAM(Sequence Alignment/Map)格式為一種通用的比對格式，用來存儲reads到參考序列的比對信息；BAM(Binary Alignment Map)是SAM的二進(jìn)制格式[21]。

1.6 流程運行

準(zhǔn)備好R軟件R程序包、參考基因組、注釋文件和探針序列文件后，用戶需要提供：

1)參考基因組名稱，如“Homo_sapiens.GRCh38.dna.primary_assembly.fa”；

2)注釋文件名稱，如“Homo_sapiens.GRCh38.94.gtf”；

3)參考基因組和注釋文件的存放路徑，如“/home/u1239/xijieprobeid/ref”；

4)GEO數(shù)據(jù)庫中的芯片平臺登錄號，如“GPL570”；

5)探針序列文件名稱，如“GPL570.id2sequence.csv”。

在對不同平臺進(jìn)行自主注釋時，用戶僅需在附件的Rmd格式文件開頭修改以上內(nèi)容，使用render()命令運行。

1.7 流程輸出文件解讀

輸出文件是探針與基因的位置信息和對應(yīng)關(guān)系，格式為csv。探針與基因的位置各用6列信息描述，列名解釋如下。

seqnames：原指序列名稱，這里指的是染色體或scaffold序號；

start：序列比對的起始位置；

end：序列比對的終止位置；

width：比對覆蓋的堿基數(shù)；

strand：染色體或scaffold的正負(fù)鏈信息；

id：基因或探針id。

2 流程測試

本文以目前應(yīng)用最廣泛、樣本量最大的兩個人類全基因組范圍表達(dá)量芯片GPL570、GPL10558和曾使用的小鼠的全基因組表達(dá)量芯片GPL21163為例，進(jìn)行重注釋；以無注釋的人類長鏈非編碼RNA表達(dá)量芯片GPL16956為例，進(jìn)行自主注釋，以測試流程的有效性。

2.1 GPL570重注釋

Human Genome U133 Plus 2.0 Array(GPL570)是Affymetrix公司的經(jīng)典產(chǎn)品，用于測定整個基因組范圍的基因表達(dá)量。自2008年問世以來廣受歡迎，且沿用至今，已有5 000多個系列、總計將近150 000個樣品的測序結(jié)果被提交到GEO數(shù)據(jù)庫，是目前樣品數(shù)最多、應(yīng)用最廣泛的基因芯片。該芯片有兩個版本的注釋文件，分別來自Affymetrix公司官網(wǎng)的注釋表格和Biocductor中的專用注釋程序包hgu133plus2.db。

該芯片設(shè)計有54 675個探針集，但每個探針集對應(yīng)的序列則有8～69條不等，總計604 258條，具體序列數(shù)統(tǒng)計結(jié)果見表2。

表2 GPL570探針集對應(yīng)的序列數(shù)統(tǒng)計Table 2 The number of sequences corresponding to the probe sets

由表2可知：絕大多數(shù)的探針集包含11條序列。在數(shù)據(jù)分析過程中發(fā)現(xiàn)，同一探針集的不同序列對應(yīng)的基因基本一致，因此完成序列比對后，探針集與基因的重復(fù)對應(yīng)關(guān)系需要去除。

使用自主注釋流程，計算得出：比對到基因組的序列數(shù)為581 910，占全部序列的比例為96.30%。最終552 760條序列成功映射到基因組，注釋表格去除重復(fù)的探針-基因映射關(guān)系后，剩余62 350條，其中有的探針對應(yīng)多個基因，有的基因?qū)?yīng)多個探針，因此分別對映射成功的探針數(shù)、映射到的基因個數(shù)進(jìn)行統(tǒng)計，并與Affymetrix公司和Biocductor中該芯片的注釋程序包hgu133plus2.db做比較，結(jié)果以韋恩圖表示(圖2)。由圖2可知：3種不同注釋共有的探針數(shù)為38 158，共有的基因數(shù)為19 234，3種注釋兩兩之間各有交集，說明3種注釋間絕大多數(shù)探針和基因的對應(yīng)關(guān)系是一致的。由于算法和依賴的參考基因組注釋版本的不同，3種注釋又各自單獨匹配到了一些不同的對應(yīng)關(guān)系，Affymetrix官網(wǎng)注釋和hgu133plus2.db程序包分別覆蓋到了41 597個(占全部探針總數(shù)的76.08%)、40 964個(占全部探針總數(shù)的74.92%)探針，并分別匹配到了22 268、21 869個基因。

值得注意的是，自主注釋流程總共注釋到了48 978個探針(占全部探針總數(shù)的89.58%)、26 963個基因，其中單獨匹配到的基因數(shù)為7 107，在原有的兩種注釋中都沒有發(fā)現(xiàn)。因此，根據(jù)基因本體論(GO)對新注釋到的編碼蛋白的基因(protein-coding gene)進(jìn)行富集分析，以驗證其正確性。

mapped_probe為比對到的探針數(shù)，mapped_gene為比對到的基因數(shù)；Bio為hguplus2.db程序包，Aff為Affymetrix官網(wǎng)注釋，Mine為自主注釋圖2 自主注釋與Affymetrix官網(wǎng)注釋及hgu133plus2.db程序包的對比Fig.2 Comparison of new annotations with Affymetrix annotations and hgu133plus2.db package

結(jié)果顯示：有411個基因成功富集到了4 275個GO條目，其中有3 178個GO條目屬于生物學(xué)過程，418個GO條目屬于細(xì)胞組分，679個GO條目屬于分子功能。這些能夠富集到GO條目的基因具有已知的生物學(xué)功能，可能會影響到表達(dá)芯片數(shù)據(jù)分析的GO富集分析結(jié)果，這也從側(cè)面說明了自主注釋的必要性。

人類基因組(HGNC)數(shù)據(jù)庫分別根據(jù)基因家族(gene family)和生物學(xué)分類(biotype)對部分基因進(jìn)行了分類。根據(jù)這兩種分類方式，分別對3種注釋匹配到的基因數(shù)量的差異進(jìn)行了比較。

選取全部的生物學(xué)分類和基因數(shù)量排名前20的基因家族統(tǒng)計結(jié)果繪制韋恩圖，結(jié)果顯示：自主注釋均比原有注釋匹配到了數(shù)量相同或者更多的基因，說明自主注釋流程幾乎在每種分類的基因中都尋找到了新的基因，這與3種注釋匹配到的基因總數(shù)比較情況相對應(yīng)。

2.2 GPL10558重注釋

HumanHT-12 V4.0 expression beadchip(GPL10558)是Illumina公司表達(dá)芯片的典型代表，可測定全基因組范圍的基因表達(dá)量，已有2 000多個系列，總計80 000多個樣品的測序結(jié)果被提交到GEO數(shù)據(jù)庫。該芯片共設(shè)計了48 107個探針，經(jīng)自主注釋，比對到參考基因組的探針數(shù)為44 302，占全部探針總數(shù)的92.10%。注釋成功的有39 226個，占全部探針總數(shù)的81.54%。注釋到的基因數(shù)為25 610個。

2.3 GPL21163重注釋

Agilent-074809 SurePrint G3 Mouse GE v2 8x60K Microarray(GPL21163)是Agilent公司生產(chǎn)的小鼠全基因組范圍的基因表達(dá)量芯片。該芯片共設(shè)計了56 745個探針，其中有153個未提供探針序列，因此有效探針數(shù)為56 592個，目前可用的探針注釋表格文件存放在GEO數(shù)據(jù)庫中，能夠注釋到46 289個探針。經(jīng)自主注釋，比對到參考基因組的探針數(shù)為52 451，占全部探針的92.68%，注釋成功的有45 692個，占探針總數(shù)的84.68%，注釋到的基因數(shù)為27 682個。

Gu等[22]使用了該芯片平臺，其排名前20的差異基因中的Ighg1基因(探針I(yè)D為A_55_P2066173，ENSAMBEL ID為ENSMUST00000103420)，是現(xiàn)有的注釋文件并未比對到的，如果直接使用現(xiàn)有注釋信息，將會影響分析結(jié)果。使用本文的自主注釋流程，能夠比對到45 692個探針，其結(jié)果文件中包含了Ighg1基因，這從側(cè)面驗證了本流程的有效性。

2.4 GPL16956自主注釋

Agilent-062918 OE Human lncRNA Microarray V4.0 028004(GPL16956)是Agilent公司于2015年生產(chǎn)的lncRNA表達(dá)芯片。目前沒有可用的探針注釋。該芯片共設(shè)計了58 944個探針，經(jīng)自主注釋，比對到參考基因組的探針數(shù)為51 869，占全部探針的88.00%。注釋成功的有31 146個，占探針總數(shù)的76.15%。注釋到的基因數(shù)為44 883個,4個測試數(shù)據(jù)統(tǒng)計匯總見表3。

由表3可知:本流程能夠較好地對多種表達(dá)芯片產(chǎn)品的探針序列進(jìn)行注釋，從比對成功率、注釋成功率、注釋到的探針數(shù)及基因數(shù)等指標(biāo)來看，均表現(xiàn)良好、穩(wěn)定，具有很好的普適性。

表3 4個測試數(shù)據(jù)統(tǒng)計匯總Table 3 Statistical summary of four test datasets

3 結(jié)論

1)本文成功開發(fā)了基于R語言的芯片探針自主注釋流程，能夠以最新的參考基因組和基因組注釋文件為基礎(chǔ)獲得最新的探針注釋，為數(shù)據(jù)挖掘工作提供了新的有力工具。

2)以GPL570為例，對本流程的自主注釋與原有注釋進(jìn)行了比較。GPL570芯片的重注釋覆蓋到的探針，基因總數(shù)及各分類基因數(shù)與原有注釋相比有所增加；經(jīng)GO富集分析驗證，本流程單獨匹配到的編碼蛋白的基因有411個能夠富集到GO條目，這部分基因具有已知的生物學(xué)功能，將會影響到表達(dá)芯片數(shù)據(jù)分析的GO富集分析結(jié)果，而原有的兩種注釋沒有比對到這些基因。一方面表明參考基因組和注釋文件較幾年前進(jìn)行了更新，比對到了原本無法比對到基因組的探針序列，另一方面證明了本流程的可靠性和時效性。

3)測試數(shù)據(jù)GPL10558、GPL21163的重注釋，分別覆蓋到了81.54%和84.68%的探針，GPL16956的自主注釋覆蓋到了76.15%的探針，能夠滿足后續(xù)分析的需求，證明了本流程的實用性和有效性。

本流程考慮到對不同芯片平臺的普適性，設(shè)置了最嚴(yán)格的比對參數(shù)，未將存在堿基錯配的探針計算在內(nèi)。在今后的研究中，筆者所在團(tuán)隊將進(jìn)一步細(xì)化流程，根據(jù)不同芯片的設(shè)計特點設(shè)置允許錯配的堿基個數(shù)，再根據(jù)實際情況設(shè)定參數(shù)，并依據(jù)本流程開發(fā)專用網(wǎng)頁工具，允許科研工作者上傳探針序列數(shù)據(jù)，獲得注釋結(jié)果文件。