专栏首页生物信息云基因芯片数据分析(二):读取芯片数据

基因芯片数据分析(二):读取芯片数据

上一篇文章(基因芯片数据分析(一):

在microarray的处理中,第一步就是读取数据。无论是自己的保存在本地的数据,还是在线保存的数据,对于不同公司的芯片可以使用不同的软件包读取。在这里,我们说的在线数据,主要是指保存在GEO (Gene Expression Omnibus) 数据库中的数据,当然GEO的数据可先下载后再读入。

读取本地数据

我们可以使用affy包中的ReadAffy函数读取cel文件。

library(affy) ##加载库文件
Data <- ReadAffy("/path/of/CELs") ##读取工作目录下的CEL文件

ReadAffy()格式如下:

ReadAffy(..., filenames=character(0),
              widget=getOption("BioC")$affy$use.widgets,
              compress=getOption("BioC")$affy$compress.cel,
              celfile.path=NULL,
              sampleNames=NULL,
              phenoData=NULL,
              description=NULL,
              notes="",
              rm.mask=FALSE, rm.outliers=FALSE, rm.extra=FALSE,
              verbose=FALSE,sd=FALSE, cdfname = NULL)

参数

说明

需要读入的文件名,可以用逗号间隔,输入多个CEL文件

filenames

文件名列表构成的字符向量(character vector)

phenoData

an AnnotatedDataFrame object, a character of length one, or a data.frame.

description

MIAME对象,它是对microarray实验的完整描述,包括名称,实验室,通信方式,实验摘要,网址,样品,杂交信息,对照,预处理,PubMedId,等等。除继承的方法外,还提供了提取相应信息的方法。

notes

注释

compress

CEL文件是否为压缩文件,支持zip和gzip

rm.mask

should the spots marked as ‘MASKS’ set to NA?

rm.outliers

should the spots marked as ‘OUTLIERS’ set to NA?

rm.extra

if TRUE, then overrides what is in rm.mask and rm.oultiers.

verbose

verbosity flag.

widget

a logical specifying if widgets should be used.

celfile.path

文件所在的目录,缺省时为R的工作目录

sampleNames

样品名列表构成的字符向量(character vector)

sd

是否读入CEL文件中的标准差?默认不读入,可以节省大量的内存。

cdfname

指定CDF库的文件名。如果设置为NULL,程序会自动从标准库中下载。

对于Affymetrix Exon/Gene ST Arrays,我们不能使用affy包来读取,我们需要使用oligo或者xps来进行分析。这里介绍oligo包。

library(oligo)
geneCELs <- list.celfiles("path/to/cels", full.names=TRUE)
affyGeneFS <- read.celfiles(geneCELs)

对于NimbleGen数据(XYS文件),可以使用oligo读取。

library(oligo)
xysFiles <- list.celfiles('myXYSs', full.names=TRUE)
rawData <- read.xysfiles(xysFiles)

对于Agilent的gpr文件或者txt文件,可以使用limma的read.maimages来读取。需要注意的是,对于不同的文件source参数有多种选择:"generic", "agilent", "agilent.median", "agilent.mean", "arrayvision", "arrayvision.ARM", "arrayvision.MTM", "bluefuse", "genepix", "genepix.custom", "genepix.median", "imagene", "imagene9", "quantarray", "scanarrayexpress", "smd.old", "smd", "spot" or "spot.close.open"。对于单色芯片,注意将green.only设置成TRUE.

library(limma)
data <- read.maimages(files=filelist, source = "agilent")

读取在线数据

在GEO数据库中保存有大量的microarray的原始数据。许多文章在发表之前,作者为了提高文章的可重复性,都会将高通量的数据提交至GEO数据库当中,以方便审稿人以及公从读者调验。

本文以GSE46106数据为例,讲述如何从GEO上下载数据。分为两种方式,第一,直接从GEO上下载表达数据,第二,直接从GEO上下载CEL文件,然后以读取本地数据的方式读入。

首先我们示例如何下载表达数据。

library(GEOquery)
gset <- getGEO("GSE46106", GSEMatrix =TRUE)
length(gset)

gset只是从geo从抓回的数据,它可能是多个数据,所以返回结果保存在了一个list当中。为了以后操作的方便,对于返回长度为1的list,我将其结果从list中抽取出来。就是类似:a <- list(a=c(1,2,3)) 以后每次访问c(1,2,3),我都要写成a[[1]]这样,感觉不方便,于是 a <- a[[1]] 这样以后访问c(1,2,3)就只需要写成a就可以了。

gset <- gset[[1]]
head(pData(gset)[,1:5])
# load NCBI platform annotation
gpl <- annotation(gset)
platf <- getGEO(gpl, AnnotGPL=TRUE)
ncbifd <- data.frame(attr(dataTable(platf), "table"))
eset <- exprs(gset)
head(eset[,1:5])
head(ncbifd[,1:5])

其次我们示例如何下载原始的CEL文件。

getGEOSuppFiles("GSE46106")
setwd("GSE46106/")
dir()
untar("GSE46106_RAW.tar")
files <- dir(pattern="gz$")
sapply(files, gunzip)
filelist <- dir(pattern="CEL$")
library(affy)
library(annotate)
data <- ReadAffy(filenames=filelist)
affydb<-annPkgName(data@annotation,type="db")
require(affydb, character.only=TRUE)
eset<-rma(data,verbose=FALSE)
eset.e <- exprs(eset)
library(annaffy)
symbols<-as.character(aafSymbol(as.character(rownames(eset)),affydb))
genes<-as.character(aafUniGene(as.character(rownames(eset)),affydb))

本文分享自微信公众号 - MedBioInfoCloud(MedBioInfoCloud),作者:DoubleHelix

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2019-11-20

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 分子对接教程 | (8) PyMOL可视化对接结果

    或者通过命令cd F:\AutoDock来实现,这与window的dos命令行和Linux系统的cd(Change Directory)命令一样。

    DoubleHelix
  • 基因芯片数据挖掘分析表达差异基因

    基因芯片(genechip)(又称DNA芯片、生物芯片)的原型是80年代中期提出的。基因芯片的测序原理是杂交测序方法,即通过与一组已知序列的核酸探针杂交进行核酸...

    DoubleHelix
  • 免疫组织化学实验

    下面是2个视频的中的操作步骤,关于免疫组化,我们之前也发过文章:免疫组化实验,可以进行参考。

    DoubleHelix
  • 资源 | 一名英文不好的程序员的救赎

    起初我并不在意这点,和大多数人一样,以为能写得一手代码,足够应付工作就行,英文好不好并不重要。

    程序员鱼皮
  • 实时音视频(TRTC)如何打印日志?

    开发者在集成 TRTCSDK 过程中,需要调试查看SDK内部打印的日志,可以参考以下代码实现日志输出。

    腾讯云-yyuanchen
  • 简述 synchronized 和 java.util.concurrent.locks.Lock 的异同 ?

    主要不同点:Lock 有比 synchronized 更精确的线程语义和更好的性能。

    MickyInvQ
  • DeepMind医疗领域大突破:AI眼病诊断工具堪比专家,准确率达94%!

    【新智元导读】Deepmind与英国NHS旗下的医院合作开发了一款AI眼部诊断工具,通过对眼部OCT图像的扫描,可识别出50多种威胁到视力的眼科疾病,准确率高达...

    新智元
  • synchronized与Lock 擂台之战

    面试官:说说synchronized和Lock(或ReentrantLock)的区别 Java 1.5之后,对共享变量访问的协调机制除了之前的synchron...

    用户1260737
  • HLS lesson1-软件工程师眼里的FPGA架构

    Vivado HLS 国内目前也是正在兴起,就我所知目前比较好的两家是华为和展讯两家ESL部门了,这是一门加速硬件设计的神器! 1.Vivado HLS的设计...

    瓜大三哥
  • 人工智能---开启未来

    转自:人民邮电报社 近年来,人工智能浪潮席卷全球。不久前,谷歌AlphaGo大胜世界围棋冠军李世石,更是引发激烈讨论。阿姆斯特朗曾经说过:“我们对下一年作出太多...

    昱良

扫码关注云+社区

领取腾讯云代金券