基于R語言的基因表達芯片注釋流程

孫小潔，鄭方強，曾健明

(1.山東農業大學 植物保護學院，山東 泰安 271018；2.珠海健明生物醫藥科技有限公司，廣東 珠海 519000)

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

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

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

準確的探針注釋是芯片數據下游分析的前提，確保能對分析結果進行正確的生物學解釋。目前的注釋存在兩個主要問題：其一是基因ID沒有一個統一的標準，每個數據庫都使用其特定的基因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]。基因ID的復雜多樣，導致已有的芯片注釋依據的基因ID也不統一；另外，芯片注釋是根據以往的參考基因組設計和比對的，而參考基因組的版本多樣，且時常更新。參考基因組存儲于Ensembl[11]、UCSC Genome Browser[12-13]以及NCBI 3個數據庫，每個數據庫中都存放了多個參考基因組版本。不同的基因芯片注釋依據的參考基因組版本不統一，更新速度較慢，有些甚至不更新。

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

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

1 系統與方法

1.1 開發環境

硬件環境：云服務器，16核心，32G內存，硬盤1T；操作系統：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 數據準備

流程的輸入文件是芯片探針序列文件，通常可以在GEO數據庫或芯片廠家官方網站下載探針平臺信息表格，刪除掉多余信息，只留下2列。第一列是探針id(Probe_id)，第二列是探針序列(Sequence)，數據結構見表1。

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

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

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

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

1.5 表達芯片探針自主注釋流程

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

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

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

1.6 流程運行

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

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

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

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

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

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

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

1.7 流程輸出文件解讀

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

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

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

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

width：比對覆蓋的堿基數；

strand：染色體或scaffold的正負鏈信息；

id：基因或探針id。

2 流程測試

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

2.1 GPL570重注釋

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

該芯片設計有54 675個探針集，但每個探針集對應的序列則有8～69條不等，總計604 258條，具體序列數統計結果見表2。

表2 GPL570探針集對應的序列數統計Table 2 The number of sequences corresponding to the probe sets

由表2可知：絕大多數的探針集包含11條序列。在數據分析過程中發現，同一探針集的不同序列對應的基因基本一致，因此完成序列比對后，探針集與基因的重復對應關系需要去除。

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

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

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

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

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

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

2.2 GPL10558重注釋

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

2.3 GPL21163重注釋

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

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

2.4 GPL16956自主注釋

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

由表3可知:本流程能夠較好地對多種表達芯片產品的探針序列進行注釋，從比對成功率、注釋成功率、注釋到的探針數及基因數等指標來看，均表現良好、穩定，具有很好的普適性。

表3 4個測試數據統計匯總Table 3 Statistical summary of four test datasets

3 結論

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

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

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

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