用DESeq2包来对RNA-seq数据进行差异分析
差异分析的套路都是差不多的,大部分设计思想都是继承limma这个包,DESeq2也不例外。
DESeq2是DESeq包的更新版本,看样子应该不会有DESeq3了,哈哈,它的设计思想就是针对count类型的数据。
可以是任意features的count数据,比如对各个基因的count,或者外显子,或者CHIP-seq的一些feature,都可以用来做差异分析。
使用这个包也是需要三个数据:
表达矩阵
分组矩阵
差异比较矩阵
总结起来三个步骤,我下面会一一讲解
重点就是构造一个dds的对象,
然后直接用DESeq函数进行normalization处理即可,
处理之后用results函数来提取差异比较结果
读取自己的数据
一般我们会自己读取表达矩阵和分组信息,下面是一个例子:
options(warn=-1)
suppressMessages(library(DESeq2))
suppressMessages(library(limma))
suppressMessages(library(pasilla))
data(pasillaGenes)
exprSet=counts(pasillaGenes)
head(exprSet)
(group_list=pasillaGenes$condition)
第一步:构建dds对象
那么根据这两个数据就可以构造dds的对象了
colData <- data.frame(row.names=colnames(exprSet),
group_list=group_list
)
dds <- DESeqDataSetFromMatrix(countData = exprSet,
colData = colData,
design = ~ group_list)
dds
可以看到我们构造的dds对象有7个样本,共14470features
从基因名可以看出,是果蝇的测序数据
我们也可以直接从expressionSet这个对象构建dds对象!
library(airway)
data(airway)
suppressMessages(library(DESeq2))
dds <- DESeqDataSet(airway, design = ~ cell+ dex)
design(dds) <- ~ dex
第二步:做normalization
suppressMessages(dds2 <- DESeq(dds))
rld <- rlogTransformation(dds2)
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)
第三步:提取差异分析结果
resultsNames(dds2)
res <- results(dds2, contrast=c("group_list","treated","untreated"))
resOrdered <- res[order(res$padj),]
resOrdered=as.data.frame(resOrdered)
head(resOrdered)
差异分析结果很容易看懂啦!
致敬
http://www.bio-info-trainee.com/bioconductor_China/software/DESeq2.html
Original: https://blog.csdn.net/weixin_46587777/article/details/124922973
Author: 皮肤科大白
Title: 用DESeq2包来对RNA-seq数据进行差异分析
原创文章受到原创版权保护。转载请注明出处:https://www.johngo689.com/697108/
转载文章受原作者版权保护。转载请注明原作者出处!