芯片数据分析,so easy?

我最早接触的高通量数据就是RNA-seq,后来接触的也基本是高通量测序结果而不是芯片数据,因此我从来没有分析过一次芯片数据,而最近有一个学员在看生信技能树在腾讯课堂发布的课程GEO数据库表达芯片处理之R语言流程(阅读原文购买)遇到了问题问我请教,为了解决这个问题,我花了一个晚上时间学习这方面的分析。

:这篇文章不会介绍R语言的安装和使用,也不会介绍GEO数据库的结构,默认你都知道,不知道可以搜索Jimmy的教程,或者购买视频学习。

数据的获取

数据获取有两种方式,R包GEOquery解析和手动下载。其中前面一种最方便,完成了手动数据下载和Bioconductor常见数据结构ExpressionSet的构造,关于这个数据结构的具体介绍看Bioconductor的介绍或者视频,简言之,就是用于存放 实验信息, 分组信息表达信息, 方便后续调用。

library(GEOquery)
gset <- getGEO("GSE13535", GSEMatrix =TRUE, AnnotGPL=TRUE )
show(gset)

一般而言GEOquery解析都是首选,除非你提供的GSE号还没被GEOquery记录或者说网络速度感人,以及你不觉得别人提供的矩阵是你所需要的,你才会决定去手工下载。分为两种情况,一种是下载赛默飞的下机原始数据格式CEL,一种是下载单个样本表达量向量或者含有所有样本的表达量矩阵。

先说第一种,可以直接点击http下载到tar打包的数据, 然后解压缩得到所有的CEL文件

setwd("F:/Project/GEO_project/")library(affy)
affy.data <- ReadAffy()
length(affy.data)
# 13
eset.rma <- rma(affy.data)
exprSet <- exprs(eset.rma)
write.table(exprSet, "expr_rma_matrix.txt", quote=F, sep="\t")
  • ReadAffy: 读取当前文件下的CEL格式文件,同时第一次还会从bioconductor上下载hugene10stv1用来注释cel文件。
  • rma: 基于robust multi-arrary average(RMA)算法衡量表达量,从而将AffyBatch对象转换成ExpressionSet
  • exprs: 获取ExpressionSet中的表达量矩阵
  • write.table: 将表达量矩阵信息保存到本地

然后是第二种,以所有样本的表达矩阵为例,可以用浏览器到ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE42nnn/GSE42589/matrix/下载,如果你会用Linux的话,可以用wget -4 ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE42nnn/GSE42589/matrix/GSE42589_series_matrix.txt.gz, 才1.7M。解压缩这个文件后,有一个txt文件, 这个txt分为两个部分。第一个部分是以”!”开头的样本的所有信息,如实验平台、处理、以及分组等信息。第二个部分则是后面的表达量信息,

expr.df <- read.table(file = "GSE42589_series_matrix.txt", header =TRUE,
comment.char = "!", row.names=1)

可以从这个角度理解这三种方法: 最开始得到的都是CEL文件,CEL文件需要一系列的步骤才能转换成表达矩阵,例如去除批次效应、质控和过滤等,得到的表达矩阵在上传时会增加元数据信息(处理方法、分组信息),就成为我们下载的GSEXXXX_series_matrix.txt.gz. 通过手工解析加R语言简单操作得到了R语言中的数据框(data.frame), 而GEOquery能够帮助我们完成下载和解析这两个步骤。

三者的优先级为:GEOquery > 手工下载表达量矩阵文件 > 手工下载原始的CEL文件。

使用limma进行差异表达分析

limma的核心函数是lmFit和eBayes, 前者是用于线性拟合,后者根据前者的拟合结果进行统计推断。

lmFit至少需要两个输入,一个是表达矩阵,一个是分组对象。

表达矩阵必须是matrix类数据结构,每一列都是存放一个样本,每一行是一个探针信息或者是注释后的基因名。这里就是向我提问的人出错的原因,他在读入数据时,read.table少了参数,row.names= 1,导致第一列是探针信息。

# 使用GEOquery
exprSet <- exprs(gset[[1]])
# 基于matrix
expr.df <- read.table(file = "GSE42589_series_matrix.txt", header =TRUE,
comment.char = "!", row.names=1)
# 从cel文件开始
exprSet <- exprs(eset.rma)

试验设计矩阵: 没有试验设计矩阵对象,limma就不知道如何比较。分组数据可以手工从之前的matrix.gz整理,整理到一个excel,然后用R读取,或者就是直接从Geoquery的结果中解析。

pData <- pData(gset[[1]])
view(pData)

其中title部分告诉了我们分组信息,2小时和18小时,每个时间段又有vehicle control, PE1.3 embolized, PE2.0 embolized,也就是2x2的双因素试验设计, 我们可以现在R语言里构建实验设计的数据框。

sample <- pData$geo_accession
treat_time <- rep(c("2h","18h"),each=11)
treat_type <- rep(rep(c("vehicle_control","PE1.3_embolized","PE2.0_embolized"), c(3,4,4)),
times=2)
design_df <- data.frame(sample, treat_time, treat_type)

根据Limma的使用手册的”9.5 Interaction Models: 2 X 2 Factorial Design”进行手续的分析。这里仅仅展示单个因素的分析过程,多个因素看文档依葫芦画瓢就行。

构建单因素试验设计矩阵,进行线性拟合

TS <- paste(design_df$treat_time, design_df$treat_type, sep=".")
TS
TS <- factor(TS, levels = unique(TS))
design <- model.matrix(~0+TS)
fit <- lmFit(exprSet, design)

然后根据我们要回答的问题,来设置比较对象。比如不同时间段下控制组哪些基因发生了差异报答,处理18小时后,处理组相对于对照组有哪些基因发生差异表达,也就是做多次t检验。

cont.matrix <- makeContrasts(
vs1  = TS18.vehicle_control-TS2.vehicle_control, # 对照组在前后的差异表达基因
vs2  = TS18.PE2.0_embolized-TS2.PE2.0_embolized, # PE2.0处理前后的差异基因
vs3  = TS18.PE1.3_embolized-TS2.PE1.3_embolized, # PE1.3在处理前后差异基因# 处理18小时候,PE2.0相对于对照变化的基因再与PE1.3与对照的差异比较
diff = (TS18.PE2.0_embolized-TS18.vehicle_control)-(TS18.PE1.3_embolized-TS18.vehicle_control),
levels = design
)

fit2 <- contrasts.fit(fit, cont.matrix)
results <- decideTests(fit2)

最后的结果可以用韦恩图展示vennDiagram(results)

更多分析

你要的R语言处理GEO芯片数据的视频出炉啦! (请务必先查看课程介绍哦

教程列表: 解读GEO数据存放规律及下载,一文就够 解读SRA数据库规律一文就够 从GEO数据库下载得到表达矩阵 一文就够 GSEA分析一文就够(单机版+R语言版) 根据分组信息做差异分析- 这个一文不够的 差异分析得到的结果注释一文就够

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2018-05-18

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

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏林浩威的专栏

使用机器学习算法打造一个简单的“微博指数”

写这篇文章的契机,是我在某天看完腾讯指数的推送后,突发奇想,想自己实现类似这样的一个东西,感觉蛮好玩的。然后就在上周末,利用了一些空余时间,写了一个简单的舆情监...

7544
来自专栏数据科学与人工智能

【机器学习】从零实现来理解机器学习算法

从零开始实现机器学习算法的好处 我推广了从零开始实现机器学习算法的观念。 我认为你可以学到很多关于算法是如何工作的。我也认为,作为一名开发者,它提供了一个学习用...

2269
来自专栏AI研习社

呵,我复现一篇深度强化学习论文容易吗

去年,OpenAI和DeepMind联手做了当时最酷的实验,不用经典的奖励信号来训练智能体,而是根据人类反馈进行强化学习的新方法。有篇博客专门讲了这个实验 Le...

1612
来自专栏流柯技术学院

性能测试学习之三—— PV->TPS转换模型&TPS波动模型

TPS = ( (80%*总PV)/(24*60*60*(T/24)))/服务器数量

2492
来自专栏互联网技术栈

七周七并发模型

以上每种模型都有各自的甜区。  这个模型适用于解决并发问题、并行问题,还是两者皆可?  这个模型适用于哪种并行架构?  这个模型是否有利于我们写出容...

852
来自专栏媒矿工厂

TICO : UHD生产系统中的Mezzanine Codec

本帖参考M. Miyazaki等人发表在SMPTE Motion Imaging Journal的文章UHDTV-2 Mezzanine Compression...

2802
来自专栏人工智能头条

聚类分析算法在Netflix服务器异常自动侦测中的应用

1652
来自专栏AI研习社

博客 | Github开源人体姿态识别项目OpenPose中文文档

人体姿态识别与估计的应用场景:抖音尬舞机、体育动作教学、3D健身教练、3D试衣、绘画辅助、游戏人物动作采集。

1.3K4

数据“厨师”ETL竞赛:今天的数据能做些什么?

它们是一个烹饪比赛的电视系列节目,享有盛名的厨师们撸起袖子,争相做出完美的菜肴。基于一个设定的主题,结合厨师们的经验,创造力和想象力,将可能有问题的食材转化为最...

2415
来自专栏Albert陈凯

R语言的优劣

https://www.zhihu.com/question/19611094 作者:艾华丰 链接:https://www.zhihu.com/questio...

2613

扫码关注云+社区

领取腾讯云代金券