基于R语言的基因表达芯片注释流程

1、基于R语言的基因表达芯片注释流程摘 要：基于R语言，将R程序包Rsubread、Rsamtools、refGenome和GenomicRanges整合为一个完整的流程，实现了 基因表达芯片探针序列的自主注释。以应用范围最广的GPL570,GPL10558和曾使用的GPL21163芯片平台为测 试数据进行重注释，并将GPL570的新注释与现存的注释做比较；对较新的长链非编码RNA表达芯片GPL16956进 行自主注释，以测试流程的实用性。结果表明：GPL570的自主注释覆盖到了 89. 58%的探针，GPL1055B、GPL21163 和GPL16956的自主注释分别覆盖到了 81. 54%、8

2、4. 68%和76. 15%的探针。在GPL570新注释单独比对到的7 107 个基因中，有411个编码蛋白的基因能够富集到GO条目，而另外两种注释未能比对到这些基因，证明了本流程的 可靠性和先进性。因此，本流程实用、有效，为数据挖掘工作提供了新的有力工具。关键词：基因表达芯片（GEO）;数据挖掘；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年代发展至今已产 生了大量的基因表达数据。如何从复杂的基因大 数据中进行知识发现，是生物信息学研究的重要课 题之一。为了满足对高通量基因表达数据存储不 断增长的需求，美国国家生物技术信息中心（NCBI） 建立了基因表达数据库（GEO），为用户提供了 可供数据提交、存储和检索的平台。目前,GEO数 据库已经收录了累计10万多个系列、280多万个样 本的数据，涉及3 000多种生物。面对海量复杂的生物数据，研究者的思维方式 也相应地从数据的生成转向对数据的深入挖掘和 分析。数据挖掘是从大量的数据中通过算法搜索 隐藏于其中信息的过程面。将数据挖掘方法应

7、用 于生物信息大数据，能够从中挖掘出有价值的信 息,寻找潜在规律，进而对相关疾病机制作出科学 的诠释，是当前生物信息学的热点问题之一。基因表达芯片是采用传统的基因表达量测定 方法，会产生出大量有价值的数据,是生物信息数 据挖掘工作的重要组成部分。基因表达芯片测序 的结果是每个样品的探针表达量，在后续分析过程 中需要根据基因与探针之间的对应关系进行ID转 换，进而计算基因的表达量高低。部分芯片平台可 以从Bioconductor网站的注释程序包中直接获取这 种对应关系,但只覆盖了约90个常用的芯片，而现 存的测序平台有10 000多个,且日益增长;也有一 些芯片平台可以从生产厂家的官方网站或GE

8、O数 据库的通用公共许可证（GPL）平台信息表格中查 找;更多芯片平台则是仅提供了探针ID与序列信 息，而未提供现成的探针与基因的对应关系项。准确的探针注释是芯片数据下游分析的前提，确 保能对分析结果进行正确的生物学解释。目前的注 释存在两个主要问题:其一是基因ID没有一个统一 的标准，每个数据库都使用其特定的基因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等。基因ID的复杂多样，导致已有 的芯片注释依据的基因ID也不统一;另外，芯片注释 是根据以往的参考基因组设计和比对的，而参考基因 组的版本多样，且时常更新。参考基因组存储于 Ensembl1、UCSC Genome Browser213以及 NCBI 3 个数据库,每个数据库中都存放了多个参考基因组版 本。不同的基因芯片注释依据的参考基因组版本不 统一,更新速度较慢，有些甚至不更新。基因芯片注释过时,ID不统一的混乱现状，使存 放在GEO数据库中大量有价值的数据无法利用起 来，给芯片数据挖掘工作带来了较大的困难，如果直 接使用过时的注释文件,势必导致后续

10、分析结果与最 新的基因注释大相径庭。因此,以最新的基因组为参 考，对探针序列进行重新注释，是芯片数据分析过程 中至关重要的工作。Yin等皿整合了多个数据库中 的斑马鱼基因注释，将Affymetrix公司的斑马鱼基因 表达芯片探针序列映射到整合的转录本中，大幅增加 了检测到的基因数量、差异基因和可变剪切数量。同 年,BarbosalNorais等发现Illumina公司提供的许 多芯片原始注释并不可靠,并针对BeadArrays系列芯 片开发了基于Perl语言的寡核苷酸芯片技术的重新 注释工具（ReMOAT） ; Arloth等句也开发了 Illumina 芯片重注释的Perl工具，使用该工具注

11、释的Human- HT12 v4芯片有约25%的探针注释与公司提供的原 始注释不同，并与ReMOAT比较发现能注释到更多 的探针。近年来，多项长链非编码RNA（lncRNA）的 差异分析研究都用到了重注释,例如非小细胞肺癌亚 型的特异性lncRNA及潜在功能分析3。本文搭建了一套简便灵活的表达芯片通用自 主注释流程，以期可以对已有注释的经典芯片平台 进行重注释,并致力于应用在无注释但提供探针序 列信息的任一表达芯片平台上。1系统与方法1.1开发环境硬件环境：云服务器,16核心，32G内存，硬盘 1T;操作系统: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数据准备流程的输入文件是芯片探针序列文件，通常可以 在GEO数据库或芯片厂家官方网站下载探针平台信 息表格，删除掉多余信息,只留下2列。第一列是探 针id（ Probe_id），第二列是探针序列（Sequence），数 据结构见表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数据库下载最新的人类参考基因组 （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，存放于同一目录下。使用 本流程需输入参考基因组和注释文件的存放路径。 L5

15、表达芯片探针自主注释流程表达芯片探针自主注释流程（图1）基于R语 言，整合了多个R程序包。先读取芯片和探针的对 应关系文件,并将其转换为fasta格式（一种序列存 储格式，是本流程使用的参考基因组序列格式。每 条序列的第一行以$”开头，跟随$%的是序列的 ID号及描述信息;第二行开始是序列内容;第二条 序列另起一行，仍然由$”开始，以此类推）。将探 针序列比对到参考基因组（也称参考序列，是一个 数字化核酸序列数据库，由科学家组装，作为一个 物种的一组基因的代表性例子：1920：），生成BAM格 式的比对结果文件，获得探针序列在基因组中的位图1基于R语言的基因表达芯片注释流程Fig. 1 An

16、R workflow for annotation of geneexpression microarray置信息;读取最新参考基因组的注释文件，获得基 因序列在基因组中的位置信息。将探针序列与基 因序列的位置信息分别转换成Grange对象（即存储 一组基因位置信息的容器，每个基因位置信息由染 色体名称、开始位置、结束位置和正点链来描述）， 寻找二者在基因组上的位置重叠区域，就获得了基 因与探针的对应关系，将其组合为一个数据框，导 出为csv格式的表格。根据参考基因组构建索引是序列比对的重要 前提，索引仅取决于参考基因组，与需注释的芯片 平台数据无关，但构建索引耗时长、需要较大的内 存，且会生

17、成约15G的大文件,是限速步骤。流程 中对该步骤进行了逻辑判断，同一物种的芯片平台 注释仅在首次运行时构建索引，不会重复构建，后 续进行其他芯片平台注释时，整个流程可在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数据库中的芯片平台登录号，如 “GPL570”；5）探针序列文件名称，如$GPL570.id2sequence.csv在对不同平台进行自主注释时,用户仅需在附 件的Rmd格式文件开头修改以上内容，使用render （）命令运行。1.7流程输出

19、文件解读输出文件是探针与基因的位置信息和对应关 系,格式为CSV。探针与基因的位置各用6列信息描 述，列名解释如下。seqnames:原指序列名称，这里指的是染色体或 scaffold 序号；start:序列比对的起始位置；end:序列比对的终止位置；width:比对覆盖的碱基数；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数据库，是目前样品 数最多、应用最广泛的基因芯片。该芯片有两个版本 的注释文件,分别来自Affymetrix公司官网的注释表格 和Biocductor中的专用注释程序包hgu133plus2. db。该芯片设计有54 675个探针集,但每个探针集 对应的

21、序列则有869条不等，总计604 258条，具 体序列数统计结果见表2。表2 GPL570探针集对应的序列数统计Table 2 The number of sequences corresponding to the probe sets序列数891011131415162069探针集数51654 130442482401由表2可知：绝大多数的探针集包含11条序 列。在数据分析过程中发现，同一探针集的不同序 列对应的基因基本一致，因此完成序列比对后，探 针集与基因的重复对应关系需要去除。使用自主注释流程，计算得出：比对到基因组 的序列数为581 910,占全部序列的比例为96. 30%。 最终

22、552 760条序列成功映射到基因组，注释表格 去除重复的探针-基因映射关系后，剩余62 350条， 其中有的探针对应多个基因，有的基因对应多个探 针,因此分别对映射成功的探针数、映射到的基因 个数进行统计，并与Affymetrix公司和Biocductor中 该芯片的注释程序包hgu133plus2. db做比较，结果 以韦恩图表示（图2）。由图2可知：3种不同注释 共有的探针数为38 158,共有的基因数为19 234, 3种注释两两之间各有交集，说明3种注释间绝大 多数探针和基因的对应关系是一致的。由于算法 和依赖的参考基因组注释版本的不同，3种注释又 各自单独匹配到了一些不同的对应关系

23、，Affymetrix 官网注释和hgu133plus2. db程序包分别覆盖到了 41 597个（占全部探针总数的76. 08%）、40 964个 （占全部探针总数的74. 92%）探针,并分别匹配到 了22 26821 869 个基因。值得注意的是，自主注释流程总共注释到了 48 978个探针（占全部探针总数的89. 58%）、26 963 个基因，其中单独匹配到的基因数为7 107,在原有 的两种注释中都没有发现。因此，根据基因本体论 （GO）对新注释到的编码蛋白的基因（protein-coding gene）进行富集分析，以验证其正确性。结果显示：有411个基因成功富集到了 4 275

24、 个GO条目，其中有3 178个GO条目属于生物学过 程,418个GO条目属于细胞组分,679个GO条目属 于分子功能。这些能够富集到GO条目的基因具有 已知的生物学功能，可能会影响到表达芯片数据分 析的GO富集分析结果，这也从侧面说明了自主注 释的必要性。人类基因组（HGNC）数据库分别根据基因家族 （gene family）和生物学分类（biotype）对部分基因进 行了分类。根据这两种分类方式，分别对3种注释 匹配到的基因数量的差异进行了比较。选取全部的生物学分类和基因数量排名前20mapped_probe为比对到的探针数,mapped_gene为比对到的基因数;Biomapped_pr

25、obe为比对到的探针数,mapped_gene为比对到的基因数;Bio为hguplus2. db程序包,Aff为Affymetrix官网注释,Mine为自主注释 图2自主注释与Affymetrix官网注释及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数据库。 该芯片共设计了 48 107个探针，经自主注释，比对 到参考基因组的探针数为44 302,占全部探针总数 的92.10%。注释成功的有39 226个，占全部探针 总数的81.54%。注释到的基因数为25 610个。2. 3 GPL21163 重注释Agilent-074809 SurePrint G3 Mouse GE 2 8x60K Microarray( GPL21163)是Agilent公司生产的小鼠全基 因组范围的基因表达量芯片。该芯片共设计了 56 745个探针，

27、其中有153个未提供探针序列,因此有 效探针数为56 592个，目前可用的探针注释表格文件 存放在GEO数据库中，能够注释到46 289个探针。 经自主注释，比对到参考基因组的探针数为52 451, 占全部探针的92. 68%，注释成功的有45 692个，占探 针总数的84. 68%，注释到的基因数为27 682个。Gu等药使用了该芯片平台，其排名前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个测试 数