用DESeq2包来对RNA-seq数据进行差异分析
差异分析的套路都是差不多的,大部分设计思想都是继承limma这个包,DESeq2也不例外。
DESeq2是DESeq包的更新版本,看样子应该不会有DESeq3了,哈哈,它的设计思想就是针对count类型的数据。
可以是任意features的count数据,比如对各个基因的count,或者外显子,或者CHIP-seq的一些feature,都可以用来做差异分析。
使用这个包也是需要三个数据:
- 表达矩阵
- 分组矩阵
- 差异比较矩阵
总结起来三个步骤,我下面会一一讲解
- 重点就是构造一个dds的对象,
- 然后直接用DESeq函数进行normalization处理即可,
- 处理之后用results函数来提取差异比较结果。
读取自己的数据
一般我们会自己读取表达矩阵和分组信息,下面是一个例子:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 | options(warn=-1) suppressMessages(library(DESeq2)) suppressMessages(library(limma)) suppressMessages(library(pasilla)) data(pasillaGenes) exprSet=counts(pasillaGenes) head(exprSet)##表达矩阵如下 ## treated1fb treated2fb treated3fb untreated1fb untreated2fb ## FBgn000000300100 ## FBgn0000008 78 46 43 47 89 ## FBgn000001420000 ## FBgn000001510101 ## FBgn0000017 3187 1672 1859 2445 4615 ## FBgn0000018369150176288383 ## untreated3fb untreated4fb ## FBgn000000300 ## FBgn0000008 53 27 ## FBgn000001410 ## FBgn000001512 ## FBgn0000017 2063 1711 ## FBgn0000018135174 (group_list=pasillaGenes$condition) ## [1] treated treated treated untreated untreated untreated untreated ## Levels: treated untreated ##这是分组信息,7个样本,3个处理的,4个未处理的对照! |
第一步:构建dds对象
那么根据这两个数据就可以构造dds的对象了
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 | colData <- data.frame(row.names=colnames(exprSet), group_list=group_list ) ## 这是一个复杂的方法构造这个对象! dds <- DESeqDataSetFromMatrix(countData = exprSet, colData = colData, design = ~ group_list) ## design 其实也是一个对象,还可以通过design(dds)来继续修改分组信息,但是一般没有必要。 dds ## class: DESeqDataSet ## dim: 14470 7 ## exptData(0): ## assays(1): counts ## rownames(14470): FBgn0000003 FBgn0000008 ... FBgn0261574 ## FBgn0261575 ## rowData metadata column names(0): ## colnames(7): treated1fb treated2fb ... untreated3fb untreated4fb ## colData names(1): group_list |
可以看到我们构造的dds对象有7个样本,共14470features
从基因名可以看出,是果蝇的测序数据
我们也可以直接从expressionSet这个对象构建dds对象!
1 2 3 4 5 6 | library(airway) data(airway) suppressMessages(library(DESeq2)) dds<- DESeqDataSet(airway, design = ~ cell+ dex) design(dds) <- ~ dex ## 构造好的对象还可以更改分组信息 |
第二步:做normalization
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 | suppressMessages(dds2 <- DESeq(dds)) ##直接用DESeq函数即可 ## 下面的代码如果你不感兴趣不需要运行,免得误导你 ## 就是看看normalization前面的数据分布差异 rld <- rlogTransformation(dds2)## 得到经过DESeq2软件normlization的表达矩阵! exprSet_new=assay(rld) par(cex = 0.7) n.sample=ncol(exprSet) if(n.sample>40) par(cex = 0.5) cols <- rainbow(n.sample*1.2) par(mfrow=c(2,2)) boxplot(exprSet, col = cols,main="expression value",las=2) boxplot(exprSet_new, col = cols,main="expression value",las=2) hist(exprSet) hist(exprSet_new) |
第三步:提取差异分析结果
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 | resultsNames(dds2) ## [1] "Intercept" "group_listtreated" "group_listuntreated" ## 其实如果只分成了两组,并没有必要指定这个比较矩阵! res <-results(dds2, contrast=c("group_list","treated","untreated")) ## 提取你想要的差异分析结果,我们这里是trt组对untrt组进行比较 resOrdered <- res[order(res$padj),] resOrdered=as.data.frame(resOrdered) head(resOrdered) ##baseMean log2FoldChange lfcSEstatpvalue ## FBgn0039155453.2753-4.281830 0.1919977 -22.30146 3.576174e-110 ## FBgn0029167 2165.0445-2.182745 0.1080670 -20.198071.017931e-90 ## FBgn0035085366.8279-2.436860 0.1505280 -16.188756.054219e-59 ## FBgn0029896257.9027-2.511257 0.1823764 -13.769643.881667e-43 ## FBgn0034736118.4074-3.166392 0.2375506 -13.329331.562878e-40 ## FBgn0040091610.6035-1.526400 0.1278555 -11.938487.457520e-33 ##padj ## FBgn0039155 2.764025e-106 ## FBgn00291673.933795e-87 ## FBgn00350851.559769e-55 ## FBgn00298967.500351e-40 ## FBgn00347362.415897e-37 ## FBgn00400919.606528e-30 |
差异分析结果很容易看懂啦!
原文来自:www.bio-info-trainee.com
相关推荐:
- NCBI在线BLAST使用方法与结果详解 2941
- 神经网络术语:Epoch、Batch Size和迭代 527
- Consed的安装与使用教程 465
- 陈连福的NGS生物信息学培训教材V2.1 277
- WGCNA分析使用教程 272
最新创建圈子
-
原料药研发及国内外注册申报
2019-01-25 10:41圈主:caolianhui 帖子:33 -
制药工程交流
2019-01-25 10:40圈主:polysciences 帖子:30 -
健康管理
2019-01-25 10:40圈主:neuromics 帖子:20 -
发酵技术
2019-01-25 10:39圈主:fitzgerald 帖子:17 -
医学肿瘤学临床试验
2019-01-25 10:39圈主:bma 帖子:58