搜索
查看: 4446|回复: 6

芯片数据分析【第四讲:根据分组信息做差异分析】

[复制链接]

13

主题

33

帖子

241

积分

版主

Rank: 7Rank: 7Rank: 7

积分
241
发表于 2016-12-14 14:32:48 | 显示全部楼层 |阅读模式
第四讲:根据分组信息做差异分析

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

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

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


可以看到,虽然都有差异,但是统计学意义不显著,原因可能是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/s ... /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之外的基因


用excel表格做差异分析
用samr包对芯片数据做差异分析 或者(http://blog.qiubio.com:8080/archives/2083/3)
用limma包对芯片数据做差异分析

用DESeq进行差异分析的源代码
用R语言的DESeq2包来对RNA-seq数据做差异分析
用limma包的voom函数来对RNA-seq数据做差异分析
http://www.bio-info-trainee.com/bioconductor_China/software/limma.html



后面的代码,估计不会R的人,简直就是看天书的了。


还有:DEG_resamplingPvalules_BioNet  和 DEG_rowttest_GENEFILTER ,当然还有自定义的批量t-test咯:https://github.com/jmzeng1314/humanid


library(ALL)
data(ALL)
library(limma)
group_list <- factor(c(substr(unlist(ALL$BT), 0, 1)))
design <- model.matrix(~ -1+ group_list )
colnames(design)<- c("B", "T")
contrast.matrix <- makeContrasts(B-T, levels=design)
contrast.matrix
fit <- lmFit(ALL, design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
tmp1 <- topTable(fit2,coef = 1,number = Inf)
tmp1$gene <- rownames(tmp1)
tmp2 <- humanid::batch_t_test(exprs(ALL),(1:length(group_list))[group_list == 'B'],(1:length(group_list))[group_list == 'T'])
tmp2$gene <- rownames(tmp2)

library(BioNet)
library(ALL)
data(ALL)
mat <- exprs(ALL)
groups <- factor(c(rep("A", 64), rep("B", 64)))
results <- resamplingPvalues(mat, groups, alternative="greater")
head(sort(results$p.values))
tmp <- humanid::batch_t_test(mat,(1:length(groups))[groups == 'A'],(1:length(groups))[groups == 'B'])
tmp$gene <- rownames(tmp)


library(BioNet)
library(DLBCL)
data(exprLym)
library(genefilter)
library(impute)
# we use rowttest from the package genefilter to analyse differential expression between the ABC and GCB subtype:
expressions <- impute.knn(exprs(exprLym))$data
t.test <- rowttests(expressions, fac=exprLym$Subgroup)
t.test[1:10, ]
t.test[order(t.test$p.value)[1:10], ]
tmp <- humanid::batch_t_test(expressions,(1:length(fac))[fac == 'GCB'],(1:length(fac))[fac == 'ABC'])






上一篇:芯片数据分析 【第三讲:对表达量矩阵用GSEA软件做分析 】
下一篇:芯片数据分析 【第五讲:对差异基因结果做GO/KEGG超几何分...
回复

使用道具 举报

634

主题

1182

帖子

4022

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4022
发表于 2017-2-9 09:18:03 | 显示全部楼层
渊梦无痕 发表于 2016-12-15 03:36
有个问题请教一下,就是我们从geo数据库下载的矩阵数据,在矩阵数据excel表头例如title,characteristics_c ...

下载原始raw去分析,必然没有分组信息,需要自己构造~
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 0 反对 1

使用道具 举报

4

主题

37

帖子

701

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
701
发表于 2016-12-15 03:36:20 | 显示全部楼层
有个问题请教一下,就是我们从geo数据库下载的矩阵数据,在矩阵数据excel表头例如title,characteristics_ch1等能看到具体分组信息。但是我们如果下载原始raw去分析,然后通过rma mas算法等处理后的得到表达矩阵看不到分组信息,难度我们要下载矩阵数据查看分组信息,然后在构建分组信息,这个如何简便方法解决?
回复 支持 反对

使用道具 举报

13

主题

33

帖子

241

积分

版主

Rank: 7Rank: 7Rank: 7

积分
241
 楼主| 发表于 2016-12-16 15:21:21 | 显示全部楼层
渊梦无痕 发表于 2016-12-15 03:36
有个问题请教一下,就是我们从geo数据库下载的矩阵数据,在矩阵数据excel表头例如title,characteristics_c ...

不用的,可以用pData这个函数,你可以自己查下;可以在R中完成的。
回复 支持 反对

使用道具 举报

4

主题

37

帖子

701

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
701
发表于 2016-12-26 01:13:59 | 显示全部楼层
Andyhigh 发表于 2016-12-16 15:21
不用的,可以用pData这个函数,你可以自己查下;可以在R中完成的。

我把的代码贴出来给你看下
ibrary(affy)
library(tcltk)
filters<-matrix(c("CEL file", ".[Cc][Ee][Ll]", "All",".*"), ncol = 2,byrow = T)
cel.files <-tk_choose.files(caption = "Select CELs", multi =TRUE,filters= filters, index = 1)
data.raw<-ReadAffy(filenames = cel.files)       #读取cel文件   GSE79973为例子
eset<-rma(data.raw)
eset
pData(eset)
pData(data.raw)
View(pData(eset))
全部没有看到之前的之前表头文件 根据特征值去提取分组信息。
回复 支持 反对

使用道具 举报

4

主题

37

帖子

701

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
701
发表于 2016-12-26 01:15:29 | 显示全部楼层
渊梦无痕 发表于 2016-12-26 01:13
我把的代码贴出来给你看下
ibrary(affy)
library(tcltk)

直接从GEO数据库下的矩阵文件,在基因表达矩阵表格上方有很多相关数据的,从那里能看到分组信息。
回复 支持 反对

使用道具 举报

3

主题

6

帖子

96

积分

注册会员

Rank: 2

积分
96
发表于 2018-5-15 15:05:36 | 显示全部楼层
【新手求助】所谓的表达差异,是将目标基因与参考基因做比较?怎么得出是上调还是下调?
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2018-12-12 05:28 , Processed in 0.035551 second(s), 28 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.