chensing的个人博客分享 http://blog.sciencenet.cn/u/chensing

博文

R/Bioconductor与转录组数据分析(6)

已有 5282 次阅读 2018-12-6 09:58 |系统分类:科研笔记| 转录组数据分析

表达差异分析

1、载入需要的包

library(limma)

library(edgeR)

2、加载数据

gse119382.express<- read.table("GSE119382expr.txt",header = T)

exprSet <- gse119382.express

rownames(exprSet) <- exprSet[,1]

exprSet <- exprSet[,-1]

3、描述分组

group_list <- c("Drought","Drought","Drought","Water","Water","Water")

4、生成实验设计矩阵:

design <- model.matrix(~0+factor(group_list))

colnames(design)=levels(factor(group_list))

rownames(design)=colnames(exprSet)

结果

design

                          Drought Water

GSM3373992_D_1       1     0

GSM3373993_D_2       1     0

GSM3373994_D_3       1     0

GSM3373998_W_1       0     1

GSM3373999_W_2       0     1

GSM3374000_W_3       0     1

attr(,"assign")

[1] 1 1

attr(,"contrasts")

attr(,"contrasts")$`factor(group_list)`

[1] "contr.treatment"

5、制作差异比较矩阵

contrast.matrix <- makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)

6、根据分组对表达矩阵进行标准化(RNA-seq数据)

exprSet.norm <- voom(exprSet,design,normalize = "quantile",plot = T)

image.png

7、用limma进行表达差异分析

线性拟合

fit <- lmFit(exprSet.norm,design)

fit2 <- contrasts.fit(fit,contrast.matrix)

经验贝叶斯统计

fit2 <- eBayes(fit2)

按p值排序输出

topTableOutput <- topTable(fit2,coef = 1, n = Inf)

过滤NA值

theDEG <- na.omit(topTableOutput)

8、显示表达差异显著的基因数(adj.p-Value < 0.05)

> results <- decideTests(fit2)

> vennDiagram(results)

image.png



https://blog.sciencenet.cn/blog-62701-1150107.html

上一篇:R/Bioconductor与转录组数据分析(5)
下一篇:生态位(ecological niche)——生态学的核心概念小记
收藏 IP: 222.197.71.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-10-4 04:13

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部