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

1、基于R語言的基因表達芯片注釋流程摘 要：基于R語言，將R程序包Rsubread、Rsamtools、refGenome和GenomicRanges整合為一個完整的流程，實現(xiàn)了 基因表達芯片探針序列的自主注釋。以應用范圍最廣的GPL570,GPL10558和曾使用的GPL21163芯片平臺為測 試數(shù)據(jù)進行重注釋，并將GPL570的新注釋與現(xiàn)存的注釋做比較；對較新的長鏈非編碼RNA表達芯片GPL16956進 行自主注釋，以測試流程的實用性。結果表明：GPL570的自主注釋覆蓋到了 89. 58%的探針，GPL1055B、GPL21163 和GPL16956的自主注釋分別覆蓋到了 81. 54%、8

2、4. 68%和76. 15%的探針。在GPL570新注釋單獨比對到的7 107 個基因中，有411個編碼蛋白的基因能夠富集到GO條目，而另外兩種注釋未能比對到這些基因，證明了本流程的 可靠性和先進性。因此，本流程實用、有效，為數(shù)據(jù)挖掘工作提供了新的有力工具。關鍵詞：基因表達芯片（GEO）;數(shù)據(jù)挖掘；R語言An R workflow for annotation of gene expression microarrayAbstract: Based on the R language，the packages Rsubread，Rsamtools，refGenome，and GenomicRa

3、nges are integrated into a complete workflow to realize the self-annotation of the microarray gene expression.The most widely applied chip platform GPL570，GPL10558 and GPL21163 used as re-annotating datasets and the new annotation of GPL570 is compared with existing one. Self-annotation of the relat

4、ively new lincRNA expression chip GPL16956 is accomplished to test the practicality of the workflow.The annotation coverage rate of GPL570 was 89. 58% whereas the rate of GPL10558，GPL21163 and GPL16956 were 81. 54%， 84.68% and 76. 15%. Among the unique 7 107 genes in this workflow，411 protein-coding

5、 gene were enriched to GO terms whereas the other two existing annotations could not，indicating the reliability and advancement of this study.Therefore，this workflow is practical and effective，and provides a new powerful tool for data mining.Keywords： gene expression microarray ( GEO) ; data mining;

6、 R langrage基因芯片技術自20世紀80年代發(fā)展至今已產 生了大量的基因表達數(shù)據(jù)。如何從復雜的基因大 數(shù)據(jù)中進行知識發(fā)現(xiàn)，是生物信息學研究的重要課 題之一。為了滿足對高通量基因表達數(shù)據(jù)存儲不 斷增長的需求，美國國家生物技術信息中心（NCBI） 建立了基因表達數(shù)據(jù)庫（GEO），為用戶提供了 可供數(shù)據(jù)提交、存儲和檢索的平臺。目前,GEO數(shù) 據(jù)庫已經收錄了累計10萬多個系列、280多萬個樣 本的數(shù)據(jù)，涉及3 000多種生物。面對海量復雜的生物數(shù)據(jù)，研究者的思維方式 也相應地從數(shù)據(jù)的生成轉向對數(shù)據(jù)的深入挖掘和 分析。數(shù)據(jù)挖掘是從大量的數(shù)據(jù)中通過算法搜索 隱藏于其中信息的過程面。將數(shù)據(jù)挖掘方法應

7、用 于生物信息大數(shù)據(jù)，能夠從中挖掘出有價值的信 息,尋找潛在規(guī)律，進而對相關疾病機制作出科學 的詮釋，是當前生物信息學的熱點問題之一。基因表達芯片是采用傳統(tǒng)的基因表達量測定 方法，會產生出大量有價值的數(shù)據(jù),是生物信息數(shù) 據(jù)挖掘工作的重要組成部分?；虮磉_芯片測序 的結果是每個樣品的探針表達量，在后續(xù)分析過程 中需要根據(jù)基因與探針之間的對應關系進行ID轉 換，進而計算基因的表達量高低。部分芯片平臺可 以從Bioconductor網(wǎng)站的注釋程序包中直接獲取這 種對應關系,但只覆蓋了約90個常用的芯片，而現(xiàn) 存的測序平臺有10 000多個,且日益增長;也有一 些芯片平臺可以從生產廠家的官方網(wǎng)站或GE

8、O數(shù) 據(jù)庫的通用公共許可證（GPL）平臺信息表格中查 找;更多芯片平臺則是僅提供了探針I(yè)D與序列信 息，而未提供現(xiàn)成的探針與基因的對應關系項。準確的探針注釋是芯片數(shù)據(jù)下游分析的前提，確 保能對分析結果進行正確的生物學解釋。目前的注 釋存在兩個主要問題:其一是基因ID沒有一個統(tǒng)一 的標準，每個數(shù)據(jù)庫都使用其特定的基因ID,主流的 有 Official _ Gene _ ID、NCBI 的 Entrez _ Gene _ ID、 Genebank GI 號、Gene Accession、RefSeq _ accession、 Ensembl_Gene_ID 等;此外還有 Vaga gene ID、

9、havana_ gene_ID、ena等?；騃D的復雜多樣，導致已有 的芯片注釋依據(jù)的基因ID也不統(tǒng)一;另外，芯片注釋 是根據(jù)以往的參考基因組設計和比對的，而參考基因 組的版本多樣，且時常更新。參考基因組存儲于 Ensembl1、UCSC Genome Browser213以及 NCBI 3 個數(shù)據(jù)庫,每個數(shù)據(jù)庫中都存放了多個參考基因組版 本。不同的基因芯片注釋依據(jù)的參考基因組版本不 統(tǒng)一,更新速度較慢，有些甚至不更新。基因芯片注釋過時,ID不統(tǒng)一的混亂現(xiàn)狀，使存 放在GEO數(shù)據(jù)庫中大量有價值的數(shù)據(jù)無法利用起 來，給芯片數(shù)據(jù)挖掘工作帶來了較大的困難，如果直 接使用過時的注釋文件,勢必導致后續(xù)

10、分析結果與最 新的基因注釋大相徑庭。因此,以最新的基因組為參 考，對探針序列進行重新注釋，是芯片數(shù)據(jù)分析過程 中至關重要的工作。Yin等皿整合了多個數(shù)據(jù)庫中 的斑馬魚基因注釋，將Affymetrix公司的斑馬魚基因 表達芯片探針序列映射到整合的轉錄本中，大幅增加 了檢測到的基因數(shù)量、差異基因和可變剪切數(shù)量。同 年,BarbosalNorais等發(fā)現(xiàn)Illumina公司提供的許 多芯片原始注釋并不可靠,并針對BeadArrays系列芯 片開發(fā)了基于Perl語言的寡核苷酸芯片技術的重新 注釋工具（ReMOAT） ; Arloth等句也開發(fā)了 Illumina 芯片重注釋的Perl工具，使用該工具注

11、釋的Human- HT12 v4芯片有約25%的探針注釋與公司提供的原 始注釋不同，并與ReMOAT比較發(fā)現(xiàn)能注釋到更多 的探針。近年來，多項長鏈非編碼RNA（lncRNA）的 差異分析研究都用到了重注釋,例如非小細胞肺癌亞 型的特異性lncRNA及潛在功能分析3。本文搭建了一套簡便靈活的表達芯片通用自 主注釋流程，以期可以對已有注釋的經典芯片平臺 進行重注釋,并致力于應用在無注釋但提供探針序 列信息的任一表達芯片平臺上。1系統(tǒng)與方法1.1開發(fā)環(huán)境硬件環(huán)境：云服務器,16核心，32G內存，硬盤 1T;操作系統(tǒng):Ubuntu 16. 04. 5。1.2 R軟件及主要程序包R 軟件版本為 3.5.

12、2,可從 https: / /mirrors.tuna. /CRAN /bin/ 獲取。R 程序包 Rsubread、Rsamtools8、refGenome 和 GenomicRanges，可從 http: / HYPERLINK file:/獲 /獲 取,也可在R語言界面使用BiocManager: : install（）命 令安裝。1.3數(shù)據(jù)準備流程的輸入文件是芯片探針序列文件，通?？梢?在GEO數(shù)據(jù)庫或芯片廠家官方網(wǎng)站下載探針平臺信 息表格，刪除掉多余信息,只留下2列。第一列是探 針id（ Probe_id），第二列是探針序列（Sequence），數(shù) 據(jù)結構見表1。表1探針序列文件格式

13、Table 1 File formats of probe sequenceProbe_id探針序列（Sequence）Probe_id_1Probe_id_2Probe_id_3GAATAAAGAACAATCTGCTGATGATCCCTCCGTGGATCTGATTCGTGTAACCATGTGATACGAGGGCGCGTAGTTTGCATTATCGTTTTTATCGTTTCAACCGACAGATGTATGTAAGGCCAACGTGCTCAAATCTTCATACAGAAAGAT推薦以逗號為分隔符，存為CSV格式,命名為 GPLxxx. id2sequence. csv”，存放于工作目錄下。1.4

14、參考基因組及注釋文件下載從Ensembl數(shù)據(jù)庫下載最新的人類參考基因組 （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，存放于同一目錄下。使用 本流程需輸入?yún)⒖蓟蚪M和注釋文件的存放路徑。 L5

15、表達芯片探針自主注釋流程表達芯片探針自主注釋流程（圖1）基于R語 言，整合了多個R程序包。先讀取芯片和探針的對 應關系文件,并將其轉換為fasta格式（一種序列存 儲格式，是本流程使用的參考基因組序列格式。每 條序列的第一行以$”開頭，跟隨$%的是序列的 ID號及描述信息;第二行開始是序列內容;第二條 序列另起一行，仍然由$”開始，以此類推）。將探 針序列比對到參考基因組（也稱參考序列，是一個 數(shù)字化核酸序列數(shù)據(jù)庫，由科學家組裝，作為一個 物種的一組基因的代表性例子：1920：），生成BAM格 式的比對結果文件，獲得探針序列在基因組中的位圖1基于R語言的基因表達芯片注釋流程Fig. 1 An

16、R workflow for annotation of geneexpression microarray置信息;讀取最新參考基因組的注釋文件，獲得基 因序列在基因組中的位置信息。將探針序列與基 因序列的位置信息分別轉換成Grange對象（即存儲 一組基因位置信息的容器，每個基因位置信息由染 色體名稱、開始位置、結束位置和正點鏈來描述）， 尋找二者在基因組上的位置重疊區(qū)域，就獲得了基 因與探針的對應關系，將其組合為一個數(shù)據(jù)框，導 出為csv格式的表格。根據(jù)參考基因組構建索引是序列比對的重要 前提，索引僅取決于參考基因組，與需注釋的芯片 平臺數(shù)據(jù)無關，但構建索引耗時長、需要較大的內 存，且會生

17、成約15G的大文件,是限速步驟。流程 中對該步驟進行了邏輯判斷，同一物種的芯片平臺 注釋僅在首次運行時構建索引，不會重復構建，后 續(xù)進行其他芯片平臺注釋時，整個流程可在3 min 以內迅速完成。其中,基因組注釋為利用生物信息 學方法和工具，對基因組所有基因的生物學功能進 行高通量注釋，包括基因識別和基因功能注釋兩個 方面，常存為gtf和gff格式如；SAM （ Sequence Alignment/Map）格式為一種通用的比對格式，用來 存儲reads到參考序列的比對信息；BAM （ Binary Alignment Map）是SAM的二進制格式。16流程運行準備好R軟件R程序包、參考基因組、

18、注釋文 件和探針序列文件后，用戶需要提供：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在對不同平臺進行自主注釋時,用戶僅需在附 件的Rmd格式文件開頭修改以上內容，使用render （）命令運行。1.7流程輸出

19、文件解讀輸出文件是探針與基因的位置信息和對應關 系,格式為CSV。探針與基因的位置各用6列信息描 述，列名解釋如下。seqnames:原指序列名稱，這里指的是染色體或 scaffold 序號；start:序列比對的起始位置；end:序列比對的終止位置；width:比對覆蓋的堿基數(shù)；strand:染色體或scaffold的正負鏈信息； id:基因或探針id。2流程測試本文以目前應用最廣泛、樣本量最大的兩個人 類全基因組范圍表達量芯片GPL570、GPL10558和 曾使用的小鼠的全基因組表達量芯片GPL21163為 例,進行重注釋；以無注釋的人類長鏈非編碼RNA 表達量芯片GPL16956為例，

20、進行自主注釋，以測試 流程的有效性。2.1 GPL570重注釋Human Genome U133 Plus 2. 0 Array（ GPL570）是 Affymetrix公司的經典產品，用于測定整個基因組范 圍的基因表達量。自2008年問世以來廣受歡迎，且 沿用至今，已有5 000多個系列、總計將近150 000個 樣品的測序結果被提交到GEO數(shù)據(jù)庫，是目前樣品 數(shù)最多、應用最廣泛的基因芯片。該芯片有兩個版本 的注釋文件,分別來自Affymetrix公司官網(wǎng)的注釋表格 和Biocductor中的專用注釋程序包hgu133plus2. db。該芯片設計有54 675個探針集,但每個探針集 對應的

21、序列則有869條不等，總計604 258條，具 體序列數(shù)統(tǒng)計結果見表2。表2 GPL570探針集對應的序列數(shù)統(tǒng)計Table 2 The number of sequences corresponding to the probe sets序列數(shù)891011131415162069探針集數(shù)51654 130442482401由表2可知：絕大多數(shù)的探針集包含11條序 列。在數(shù)據(jù)分析過程中發(fā)現(xiàn)，同一探針集的不同序 列對應的基因基本一致，因此完成序列比對后，探 針集與基因的重復對應關系需要去除。使用自主注釋流程，計算得出：比對到基因組 的序列數(shù)為581 910,占全部序列的比例為96. 30%。 最終

22、552 760條序列成功映射到基因組，注釋表格 去除重復的探針-基因映射關系后，剩余62 350條， 其中有的探針對應多個基因，有的基因對應多個探 針,因此分別對映射成功的探針數(shù)、映射到的基因 個數(shù)進行統(tǒng)計，并與Affymetrix公司和Biocductor中 該芯片的注釋程序包hgu133plus2. db做比較，結果 以韋恩圖表示（圖2）。由圖2可知：3種不同注釋 共有的探針數(shù)為38 158,共有的基因數(shù)為19 234, 3種注釋兩兩之間各有交集，說明3種注釋間絕大 多數(shù)探針和基因的對應關系是一致的。由于算法 和依賴的參考基因組注釋版本的不同，3種注釋又 各自單獨匹配到了一些不同的對應關系

23、，Affymetrix 官網(wǎng)注釋和hgu133plus2. db程序包分別覆蓋到了 41 597個（占全部探針總數(shù)的76. 08%）、40 964個 （占全部探針總數(shù)的74. 92%）探針,并分別匹配到 了22 26821 869 個基因。值得注意的是，自主注釋流程總共注釋到了 48 978個探針（占全部探針總數(shù)的89. 58%）、26 963 個基因，其中單獨匹配到的基因數(shù)為7 107,在原有 的兩種注釋中都沒有發(fā)現(xiàn)。因此，根據(jù)基因本體論 （GO）對新注釋到的編碼蛋白的基因（protein-coding gene）進行富集分析，以驗證其正確性。結果顯示：有411個基因成功富集到了 4 275

24、 個GO條目，其中有3 178個GO條目屬于生物學過 程,418個GO條目屬于細胞組分,679個GO條目屬 于分子功能。這些能夠富集到GO條目的基因具有 已知的生物學功能，可能會影響到表達芯片數(shù)據(jù)分 析的GO富集分析結果，這也從側面說明了自主注 釋的必要性。人類基因組（HGNC）數(shù)據(jù)庫分別根據(jù)基因家族 （gene family）和生物學分類（biotype）對部分基因進 行了分類。根據(jù)這兩種分類方式，分別對3種注釋 匹配到的基因數(shù)量的差異進行了比較。選取全部的生物學分類和基因數(shù)量排名前20mapped_probe為比對到的探針數(shù),mapped_gene為比對到的基因數(shù);Biomapped_pr

25、obe為比對到的探針數(shù),mapped_gene為比對到的基因數(shù);Bio為hguplus2. db程序包,Aff為Affymetrix官網(wǎng)注釋,Mine為自主注釋 圖2自主注釋與Affymetrix官網(wǎng)注釋及hgul33plus2. db程序包的對比Fig. 2 Comparison of new annotations with Affymetrix annotations and hgul33plus2. db package2. 2 GPL10558 重注釋HumanHT-12 V4.0 expression beadchip ( GPL10558) 是Illumina公司表達芯片的典型代

26、表，可測定全基 因組范圍的基因表達量，已有2 000多個系列，總計 80 000多個樣品的測序結果被提交到GEO數(shù)據(jù)庫。 該芯片共設計了 48 107個探針，經自主注釋，比對 到參考基因組的探針數(shù)為44 302,占全部探針總數(shù) 的92.10%。注釋成功的有39 226個，占全部探針 總數(shù)的81.54%。注釋到的基因數(shù)為25 610個。2. 3 GPL21163 重注釋Agilent-074809 SurePrint G3 Mouse GE 2 8x60K Microarray( GPL21163)是Agilent公司生產的小鼠全基 因組范圍的基因表達量芯片。該芯片共設計了 56 745個探針，

27、其中有153個未提供探針序列,因此有 效探針數(shù)為56 592個，目前可用的探針注釋表格文件 存放在GEO數(shù)據(jù)庫中，能夠注釋到46 289個探針。 經自主注釋，比對到參考基因組的探針數(shù)為52 451, 占全部探針的92. 68%，注釋成功的有45 692個，占探 針總數(shù)的84. 68%，注釋到的基因數(shù)為27 682個。Gu等藥使用了該芯片平臺，其排名前20的差 異基因中的Ighg1基因(探針I(yè)D為A_55_P2066173, ENSAMBEL ID 為 ENSMUST00000103420)，是現(xiàn)有的 注釋文件并未比對到的,如果直接使用現(xiàn)有注釋信 息,將會影響分析結果。使用本文的自主注釋流程, 能夠比對到45 692個探針，其結果文件中包含了 Ighg1 基因，這從側面驗證了本流程的有效性2.4 GPL16956自主注釋Agilent - 062918 OE Human lncRNA Microarray V4. 0 028004( GPL16956)是 Agilent 公司于 2015 年 生產的lncRNA表達芯片。目前沒有可用的探針注 釋。該芯片共設計了 58 944個探針，經自主注釋, 比對到參考基因組的探針數(shù)為51 869，占全部探針 的88. 00%。注釋成功的有31 146個，占探針總數(shù) 的76. 15%。注釋到的基因數(shù)為44 883個,4個測試 數(shù)