如何使用GEOquery和limma完成芯片数据的差异表达分析-阿里云开发者社区

开发者社区> 徐洲更> 正文

如何使用GEOquery和limma完成芯片数据的差异表达分析

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

如何分析芯片数据

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

数据的获取

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

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

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

img_7a09e1acaeee919dfc151ee56fe6e78b.jpe
数据下载

先说第一种,可以直接点击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分为两个部分。第一个部分是以"!"开头的样本的所有信息,如实验平台、处理、以及分组等信息。第二个部分则是后面的表达量信息,

img_5d6bfce7c146b3059634522ec3e06051.jpe
Series Matrix Files
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)
img_e82bea7cf46a5037a0c8c78b446ef338.jpe
GEOquery解析的信息

其中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)

更多分析

找到的差异表达基因后续要做GO/KEGG分析,可以在生信技能树公众号中搜索,要是基础太差,就付费购买GEO数据库表达芯片处理之R语言流程吧。

版权声明:本文内容由阿里云实名注册用户自发贡献,版权归原作者所有,阿里云开发者社区不拥有其著作权,亦不承担相应法律责任。具体规则请查看《阿里云开发者社区用户服务协议》和《阿里云开发者社区知识产权保护指引》。如果您发现本社区中有涉嫌抄袭的内容,填写侵权投诉表单进行举报,一经查实,本社区将立刻删除涉嫌侵权内容。

相关文章
使用NAT网关轻松为单台云服务器设置多个公网IP
在应用中,有时会遇到用户询问如何使单台云服务器具备多个公网IP的问题。 具体如何操作呢,有了NAT网关这个也不是难题。
26064 0
阿里云服务器ECS远程登录用户名密码查询方法
阿里云服务器ECS远程连接登录输入用户名和密码,阿里云没有默认密码,如果购买时没设置需要先重置实例密码,Windows用户名是administrator,Linux账号是root,阿小云来详细说下阿里云服务器远程登录连接用户名和密码查询方法
10388 0
阿里云服务器端口号设置
阿里云服务器初级使用者可能面临的问题之一. 使用tomcat或者其他服务器软件设置端口号后,比如 一些不是默认的, mysql的 3306, mssql的1433,有时候打不开网页, 原因是没有在ecs安全组去设置这个端口号. 解决: 点击ecs下网络和安全下的安全组 在弹出的安全组中,如果没有就新建安全组,然后点击配置规则 最后如上图点击添加...或快速创建.   have fun!  将编程看作是一门艺术,而不单单是个技术。
9757 0
使用OpenApi弹性释放和设置云服务器ECS释放
云服务器ECS的一个重要特性就是按需创建资源。您可以在业务高峰期按需弹性的自定义规则进行资源创建,在完成业务计算的时候释放资源。本篇将提供几个Tips帮助您更加容易和自动化的完成云服务器的释放和弹性设置。
11738 0
windows server 2008阿里云ECS服务器安全设置
最近我们Sinesafe安全公司在为客户使用阿里云ecs服务器做安全的过程中,发现服务器基础安全性都没有做。为了为站长们提供更加有效的安全基础解决方案,我们Sinesafe将对阿里云服务器win2008 系统进行基础安全部署实战过程! 比较重要的几部分 1.
8285 0
阿里云服务器如何登录?阿里云服务器的三种登录方法
购买阿里云ECS云服务器后如何登录?场景不同,阿里云优惠总结大概有三种登录方式: 登录到ECS云服务器控制台 在ECS云服务器控制台用户可以更改密码、更换系.
11480 0
阿里云服务器ECS登录用户名是什么?系统不同默认账号也不同
阿里云服务器Windows系统默认用户名administrator,Linux镜像服务器用户名root
3133 0
阿里云服务器如何登录?阿里云服务器的三种登录方法
购买阿里云ECS云服务器后如何登录?场景不同,云吞铺子总结大概有三种登录方式: 登录到ECS云服务器控制台 在ECS云服务器控制台用户可以更改密码、更换系统盘、创建快照、配置安全组等操作如何登录ECS云服务器控制台? 1、先登录到阿里云ECS服务器控制台 2、点击顶部的“控制台” 3、通过左侧栏,切换到“云服务器ECS”即可,如下图所示 通过ECS控制台的远程连接来登录到云服务器 阿里云ECS云服务器自带远程连接功能,使用该功能可以登录到云服务器,简单且方便,如下图:点击“远程连接”,第一次连接会自动生成6位数字密码,输入密码即可登录到云服务器上。
21131 0
+关注
徐洲更
生信媛公众号编辑、生信必修课之软件安装课程作者
284
文章
0
问答
文章排行榜
最热
最新
相关电子书
更多
《2021云上架构与运维峰会演讲合集》
立即下载
《零基础CSS入门教程》
立即下载
《零基础HTML入门教程》
立即下载