用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)##表达矩阵如下

## 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的对象了

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对象！

library(airway)

data(airway)

suppressMessages(library(DESeq2))

dds<- DESeqDataSet(airway, design = ~ cell+ dex)

design(dds) <- ~ dex

## 构造好的对象还可以更改分组信息

第二步：做normalization

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)

第三步：提取差异分析结果

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