根据分组信息做差异分析- 这个一文不够的

通过前面的讲解,我们顺利的了解了GEO数据库以及如何下载其数据,得到我们想要的表达矩阵,也成功的使用了GSEA这个分析套路。

历史目录: 解读GEO数据存放规律及下载,一文就够 解读SRA数据库规律一文就够 从GEO数据库下载得到表达矩阵 一文就够 GSEA分析一文就够(单机版+R语言版)

但最常用的其实是差异分析,下面我们来细细讲解:

提到表达量数据分析,不管是通过芯片技术还是高通量测序技术得到的表达量矩阵,我们都需要根据样本的分组信息来对所检测到的所有基因或者蛋白分子来做差异分析,想找到显著性变化的生物大分子。而我们通常会用封装好的包来替我们完成这一过程,比如limma,samr,genefilter,BioNet,deseq,edgeR等等.

它们的本质就是对表达量矩阵做一个归一化,让不同组样本的表达量具有可比性,然后利用理想的统计分布检验函数来计算差异的显著性。

最简单的T检验

img

上面是一个简单的例子,抽取两个基因的表达量来做差异分析,选取最简单的T检验来做,在R里面完成如下:

img

对这两个基因的表达量检验结果如下;

img

可以看到,虽然都有差异,但是统计学意义不显著,原因可能是case组里面的样本表达量内部差异太大。

其实如果已经拿到了表达矩阵,直接在excel里面也可以进行T检验,但是芯片数据,现在比较流行的limma这个R包,封装好的差异分析函数来做。上面代码的文字版如下:

exprSet <-matrix(c(10,7,1,9,1,3,5,5,
                   3,6,8,3,2,6,5,5),
                 nrow=2,byrow = T
                 )
## 构造表达矩阵,一般是读取表达矩阵文件
# ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE77nnn/GSE77705/suppl/GSE77705_RPKM.xlsx 
## 比如下载上面的表格,另存为csv格式,用read.csv读取到R里面
group_list <- c(rep('case',4),rep('control',4))
## 自己看看表达矩阵是如何分组的,这里是前面4列是case,后面4列是control
t.test(exprSet[1,]~group_list)
t.test(exprSet[2,]~group_list)
## 依次提取表达矩阵的每一行基因的表达值,根据分组信息做T检验

对所有基因都依次做了T检验之后,就要根据检验结果来挑选差异基因了。这个挑选的过程不确定性很大,但是一般要综合考虑表达量的变化值(fold change)和统计显著性量度(p值、q值等)。

如果芯片测到了两万个基因的表达量,一般要挑选500~1000个左右的差异基因,我一般选择p值在0.05一下,表达量的变化值(fold change)在2倍的sd之外的基因。关于这个最基础的差异分析,我在我的生信菜鸟团博客写过无数阅读超10000的教程,大家可以自行检索慢慢品味:

  • 用excel表格做差异分析
  • 用samr包对芯片数据做差异分析
  • 用limma包对芯片数据做差异分析
  • 用DESeq进行差异分析的源代码
  • 用R语言的DESeq2包来对RNA-seq数据做差异分析
  • 用limma包的voom函数来对RNA-seq数据做差异分析

使用成熟的R包来做差异分析-limma

用基因芯片的手段来探针基因表达量的技术虽然已经在逐步被RNA-seq技术取代,但毕竟经历了十多年的发展了,在GEOarrayexpress数据库里面存储的全球研究者数据都已经超过了50PB了!实在是很可观,里面还是有非常多等待挖掘的地方! 现在我们要讲的就是基因表达芯片数据的一种分析方式,差异分析,主角就是limma这个package,虽然一直都有各种差异分析统计方法提出来,但是limma包绝对是其中的佼佼者,我们也不多说废话了,直接上教程!

使用这个包需要三个数据:

  • 表达矩阵
  • 分组矩阵
  • 差异比较矩阵

而且总共也只有三个步骤:

  • lmFit
  • eBayes
  • topTable

首先是表达矩阵,这个非常简单了,你可以自己从芯片原始的cel文件开始用affy包读取,然后用rma或者mas5函数做归一化处理,最后得到表达矩阵。也可以直接下载表达矩阵然后用read.table等函数读取到R里面。或者直接下载一个含有表达矩阵的数据对象!

suppressPackageStartupMessages(library(CLL))
data(sCLLex)
exprSet=exprs(sCLLex)   ##sCLLex是依赖于CLL这个package的一个对象
samples=sampleNames(sCLLex)
pdata=pData(sCLLex)
group_list=as.character(pdata[,2])
dim(exprSet)
## [1] 12625    22
exprSet[1:5,1:5]
##           CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL CLL15.CEL
## 1000_at    5.743132  6.219412  5.523328  5.340477  5.229904
## 1001_at    2.285143  2.291229  2.287986  2.295313  2.662170
## 1002_f_at  3.309294  3.318466  3.354423  3.327130  3.365113
## 1003_s_at  1.085264  1.117288  1.084010  1.103217  1.074243
## 1004_at    7.544884  7.671801  7.474025  7.152482  6.902932
group_list  ##这个也是存储在sCLLex对象里面的信息,描述了样本的状态,需要根据它来进行样本分类
##  [1] "progres." "stable"   "progres." "progres." "progres." "progres."
##  [7] "stable"   "stable"   "progres." "stable"   "progres." "stable"  
## [13] "progres." "stable"   "stable"   "progres." "progres." "progres."
## [19] "progres." "progres." "progres." "stable"

可以看到我们得到了表达矩阵exprSet,它列是各个样本,行是每个探针ID,一个纯粹的表达矩阵,必须是数字型的! 我们可以简单做该表达矩阵做一个QC检测

par(cex = 0.7)
n.sample=ncol(exprSet)
if(n.sample>40) par(cex = 0.5)
cols <- rainbow(n.sample*1.2)
boxplot(exprSet, col = cols,main="expression value",las=2)

img

通过这些boxplot可以看到各个芯片直接数据还算整齐,可以进行差异比较

然后我们制作分组矩阵

suppressMessages(library(limma))
## Warning: package 'limma' was built under R version 3.1.1
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
design
##           progres. stable
## CLL11.CEL        1      0
## CLL12.CEL        0      1
## CLL13.CEL        1      0
## CLL14.CEL        1      0
## CLL15.CEL        1      0
## CLL16.CEL        1      0
## CLL17.CEL        0      1
## CLL18.CEL        0      1
## CLL19.CEL        1      0
## CLL20.CEL        0      1
## CLL21.CEL        1      0
## CLL22.CEL        0      1
## CLL23.CEL        1      0
## CLL24.CEL        0      1
## CLL2.CEL         0      1
## CLL3.CEL         1      0
## CLL4.CEL         1      0
## CLL5.CEL         1      0
## CLL6.CEL         1      0
## CLL7.CEL         1      0
## CLL8.CEL         1      0
## CLL9.CEL         0      1
## attr(,"assign")
## [1] 1 1
## attr(,"contrasts")
## attr(,"contrasts")$`factor(group_list)`
## [1] "contr.treatment"

这个design就是我们制作好的分组矩阵,需要根据我们下载的芯片数据的实验设计方案来,此处例子是CLL疾病的探究,22个样本分成了两组,你们自己的数据只需要按照同样的方法制作即可! 最后一个差异比较矩阵,此处不能省略,即使这里的数据只分成了两组,所以肯定是它们两组直接进行比较啦!如果分成多于两个组,那么就需要制作更复杂的差异比较矩阵。

contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
contrast.matrix ##这个矩阵声明,我们要把progres.组跟stable进行差异分析比较
##           Contrasts
## Levels     progres.-stable
##   progres.               1
##   stable                -1

那么我们已经制作好了必要的输入数据,下面开始讲如何使用limma这个包来进行差异分析了!

##step1
fit <- lmFit(exprSet,design)
##step2
fit2 <- contrasts.fit(fit, contrast.matrix) ##这一步很重要,大家可以自行看看效果
fit2 <- eBayes(fit2)  ## default no trend !!!
##eBayes() with trend=TRUE
##step3
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput) 
#write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
head(nrDEG)
##               logFC  AveExpr         t      P.Value  adj.P.Val        B
## 39400_at -1.0284628 5.620700 -5.835799 8.340576e-06 0.03344118 3.233915
## 36131_at  0.9888221 9.954273  5.771526 9.667514e-06 0.03344118 3.116707
## 33791_at  1.8301554 6.950685  5.736161 1.048765e-05 0.03344118 3.051940
## 1303_at  -1.3835699 4.463438 -5.731733 1.059523e-05 0.03344118 3.043816
## 36122_at  0.7801404 7.259612  5.141064 4.205709e-05 0.10619415 1.934581
## 36939_at  2.5471980 6.915045  5.038301 5.362353e-05 0.11283285 1.736846

可以看到就三个步骤!最后输出的nrDEG就是我们的差异分析结果! 至于结果该如何解释,大家可以仔细阅读说明书啦!

下面给大家一些我以前收费写的教程,大放送咯:

  • http://www.bio-info-trainee.com/bioconductor_China/software/limma.html
  • https://github.com/bioconductor-china/basic/blob/master/makeContrasts.md
  • http://www.bio-info-trainee.com/bioconductor_China/software/DESeq2.html
  • http://www.bio-info-trainee.com/bioconductor_China/software/limma_voom.html

写在最后

为什么说一文不够呢?

其实很简单,还有很多复杂的分组以及比较我没办法详细讲解,字数受限,包括配对样本的差异分析,时间序列的分析,还有不同差异分析结果的比较。

其实,只要入了门,后面的分析,无非就是时间问题,当然,需要用心的钻研。

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

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

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

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏新智元

谷歌大规模机器学习:模型训练、特征工程和算法选择 (32PPT下载)

【新智元导读】在 ThingsExpo 会议上,谷歌软件工程师 Natalia Ponomareva 作了有关如何在大规模机器学习中取得成功的讲座。Natali...

39010
来自专栏编程

2017年11月R新包推荐

一. 文档概述 11月份,在R官方(CRAN)共计发布了237个新包,本文选摘了40个R包,包含以下几个类别:计算方法、数据、数据科学、科学、社会科学、工具及可...

2028
来自专栏语言、知识与人工智能

【腾讯知文】任务型对话机器人简介

2K11
来自专栏PPV课数据科学社区

32页ppt干货|谷歌大规模机器学习:模型训练、特征工程和算法选择

谷歌机器学习:实际应用技巧 ? ? 什么是机器学习(ML)? 从概念上讲:给定(训练)数据,发现一些潜在的模式并将这个模式应用于新数据。 ML 的类型:监督学习...

49010
来自专栏AI研习社

机器学习者必知的 5 种深度学习框架

本文为雷锋字幕组编译的技术博客,原标题 The 5 Deep Learning Frameworks Every Serious Machine Learner...

1543
来自专栏腾讯技术工程官方号的专栏

FPGA异构计算在图片处理上的应用以及HEVC算法原理介绍

作者介绍:chaningwang,2008年毕业于中国科学院研究生院,主攻FPGA高性能计算、FPGA图像处理等方向。 先后在华为、怡化公司从事FPGA开发...

3786
来自专栏企鹅号快讯

强化学习从入门到放弃

重要概念 强化学习(REinforcement Learning)(个人理解):在设定的规则下,通过训练让机器学习完成特定的任务。 强化学习的目的是学习一个策略...

3325
来自专栏AI科技评论

郑文琛:基于网络功能模块的图特征学习 | AI 研习社79期大讲堂

AI研习社按:图是一种常见的数据结构,可以被用于许多不同的预测任务。如何从图数据学习有效特征是个重要的问题。我们的新概念是从点和边出发,拓展到更高阶的子图结构(...

1094
来自专栏追不上乌龟的兔子

Neo4j中的图形算法:15种不同的图形算法及其功能

只有你拥有使用图形分析的技巧,并且图形分析能快速提供你需要的见解时,它才具有价值。因而最好的图形算法易于使用,快速执行,并且产生有权威的结果。

2.8K3
来自专栏PPV课数据科学社区

只需七步就能掌握Python数据准备

摘要: 本文主要讲述了如何在python中用七步就能完成中数据准备。 上图为CRISP-DM模型中的数据准备   下面七个步骤涵盖了数据准备的概念,个别任务...

3257

扫码关注云+社区