# 差异分析|DESeq2完成配对样本的差异分析

```# 加载包
library(openxlsx)
library(DESeq2)
library(limma)
library(edgeR)```

1.读入原始数据及分组信息

```rowdata <- read.xlsx("Count_data.xlsx",sheet = 1,rowNames = T)
gset <- rowdata[rowMeans(rowdata)>0,] # 剔除表达量低的基因

group_list <- c(rep("case",5), rep("control",5))
group_list <- factor(group_list,levels = c("control","case"))```

```## 2.1表达矩阵
data = gset

## 2.2分组矩阵
design <- model.matrix(~0+group_list)
rownames(design) = colnames(data)
colnames(design) <- levels(group_list)

## 2.3差比表达矩阵构建，并过滤
DGElist <- DGEList( counts = data, group = group_list)
keep <- rowSums(cpm(DGElist) > 0.5 ) >=2
table(keep)
# FALSE 2443 TRUE 15909

DGElist_QC <- DGElist[keep, ,keep.lib.sizes=FALSE]
dim(DGElist)
# 18352 10 过滤前

dim(DGElist_QC)
# 15909 10过滤后

## 2.4归一化基因表达分布

DGElist_norm <- calcNormFactors(DGElist_QC, method = "TMM")
DGElist_norm\$samples\$norm.factors # 查看每个样品归一化后的系数
DGElist_QC\$samples\$norm.factors #未归一化之前的样本系数，都是1

## 2.5 limma包进行voom函数
v <- voom(DGElist_norm, design, plot = TRUE, normalize = "quantile")
fit <- lmFit(v, design)
cont.matrix <- makeContrasts(contrasts = c('case-control'), levels = design)
fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)

## 2.6提取差异矩阵
nrDEG_limma_voom = topTable(fit2, coef = "case-control", n = Inf)
nrDEG_limma_voom = na.omit(nrDEG_limma_voom)
nrDEG_limma_voom = nrDEG_limma_voom[order(nrDEG_limma_voom\$logFC),]

## 2.7 定义差异基因
nrDEG <- nrDEG_limma_voom
nrDEG\$Group = "notsignificant"

logFC_cutoff <- 0.6 # 定义差异基因的标准，自定义

nrDEG\$Group [which( (nrDEG\$P.Value < 0.05) & (nrDEG\$logFC > logFC_cutoff) )] = "upregulated"
nrDEG\$Group [which( (nrDEG\$P.Value < 0.05) & (nrDEG\$logFC < -logFC_cutoff) )] = "downregulated"

table(nrDEG\$Group)```

OK，尝试使用DESeq2包的非配对差异分析。

```## 3.1表达矩阵
data = apply(gset, 2, as.integer) ## DESeq2分析需要是整数
row.names(data) <- row.names(gset)

## 3.2分组矩阵
condition = group_list
coldata <- data.frame(row.names = colnames(data), condition)

dds <- DESeqDataSetFromMatrix(countData = data,
colData = coldata,
design = ~condition)

dds\$condition<- relevel(dds\$condition, ref = "control") # 指定哪一组作为对照组

## 3.3差异表达矩阵

dds <- DESeq(dds)
nrDEG_DESeq2 <- as.data.frame(results(dds))
nrDEG_DESeq2 = nrDEG_DESeq2[order(nrDEG_DESeq2\$log2FoldChange),]

## 3.4定义差异基因
nrDEG <- nrDEG_DESeq2
nrDEG\$Group = "notsignificant"
logFC_cutoff <- 0.6

nrDEG\$Group[which( (nrDEG\$pvalue < 0.05) & (nrDEG\$log2FoldChange > logFC_cutoff) )] = "upregulated"
nrDEG\$Group[which( (nrDEG\$pvalue < 0.05) & (nrDEG\$log2FoldChange < -logFC_cutoff) )] = "downregulated"

table(nrDEG\$Group)```

```## 4.1表达矩阵

data = apply(gset, 2, as.integer) ## DESeq2分析需要是整数
row.names(data) <- row.names(gset)

## 4.2分组矩阵，配对分析与常规分析最大的区别就在分组矩阵

condition = group_list
# 配对分析要加上这段代码，知道谁和谁是一对，比如1,1是一对，5,5是一对
subject <- factor(c(1,2,3,4,5,1,2,3,4,5))

coldata <- data.frame(row.names = colnames(data), condition)

# 注意在design中加上配对信息
dds <- DESeqDataSetFromMatrix(countData = data,
colData = coldata,
design = ~subject +condition)

dds\$condition<- relevel(dds\$condition, ref = "control")

## 4.3差异表达矩阵，还是和常规分析一样

dds <- DESeq(dds)
nrDEG_DESeq2 <- as.data.frame(results(dds))
rld <- rlog(dds)
# 这里我还提取了标准化后的表达矩阵，可以用于后续的热图绘制等等
normal_gset <- assay(rld)
nrDEG_DESeq2 = nrDEG_DESeq2[order(nrDEG_DESeq2\$log2FoldChange),]

## 4.4定义差异基因

nrDEG <- nrDEG_DESeq2
nrDEG\$Group = "notsignificant"

logFC_cutoff <- 0.6
nrDEG\$Group[which( (nrDEG\$pvalue < 0.05) & (nrDEG\$log2FoldChange > logFC_cutoff) )] = "upregulated"
nrDEG\$Group[which( (nrDEG\$pvalue < 0.05) & (nrDEG\$log2FoldChange < -logFC_cutoff) )] = "downregulated"

table(nrDEG\$Group)```

0 条评论

• ### DESeq2差异表达分析

在前文scRNA-seq marker identification(二)，我们我们提到了差异分析，下面我们来详细了解下

• ### DESeq2差异表达分析(二)

DESeq2工作流程的下一步是QC，它包括样本级和基因级的步骤，对计数数据执行QC检查，以帮助我们确保样本/重复 看起来很好。

• ### DESeq2转录组差异表达分析实例

使用library(DESeq2)加载的时候遇到报错 ：载入了名字空间‘rlang’ 0.4.0，但需要的是>= 0.4.2 解决办法：将rlang包手动删除...

• ### 使用DESeq2进行两组间的差异分析

DESeq2 接受raw count的定量表格，然后根据样本分组进行差异分析，具体步骤如下

• ### 基因芯片数据分析（八）：DESeq2差异分析实战案例

基因芯片数据分析（六）：DESeq2包的基本原理 我们接下来通过一个案例介绍利用edgeR进行差异分析。

• ### 转录组分析 | 使用DESeq2进行基因差异表达分析

通过RSEM我们获取了样本中每个基因的counts和表达量，接下来使用tximport校正不同样本间基因长度的差异。

• ### 不编程就能完成差异表达分析

上几周，生信技能树的元老级人物果子发了一篇重量级文章数据库一网打尽：不会编程又怎样，还不照样拿课题发文章！ 阅读量逼近4000，还超过了Jimmy的署名文一个全...

• ### 比较不同流程（limma/voom，edgeR，DESeq2 ）差异分析的区别

距离第一次听说生信已经十几年了，现在是邋遢大叔重新开始学代码，精力确实已不像从前，各位入坑还是要乘早。后来约莫在5年前，课题组当时有个RNA-Seq数据，lab...