RNA-SEQ（三）：DEseq2包进行差异分析。

目的：DEseq2是基于负二项分布模型，寻找组间显著表达变化的基因，以解释基因表达水平的变化。

内容：1. 将样本分为对照组和实验组。

2. 输入read counts矩阵进行差异分析。

数据：RNASEQ上游分析得到的read count矩阵文件，以及分组信息的表格文件。

工具：Rstudio。

步骤：

1. 确认对照组和处理组：D6 VS D9，写分组信息储存为csv文件，如D6 VS D9的分组文件如下：

image.png

2. counts矩阵包含所有样本，也要挑出比对的两个样本，生成新的矩阵，文件如下。将矩阵文件与分组文件合并，进行差异分析，脚本如下：

###多样本差异分析Deseq2
if(!requireNamespace("BiocManager",quietly=TRUE))
install.packages("BiocManager")
BiocManager::install("DESeq2")

library(tidyverse)
library(DESeq2)#差异分析包
#import data
setwd("E:/hpcs/f")
#打开counts文件
mycounts<-read.table("counts_matrix_symble.txt", header=TRUE)#header行名
head(mycounts)

#删除重复基因
rows <- rownames(unique(mycounts['Geneid']))
mycounts <- mycounts[rows,]

#修改行名
rownames(mycounts)<-mycounts[,1]
mycounts <- mycounts[,-1]
head(mycounts)

#载入分组信息
colData <-read.csv("D6_VS_D9_condition.csv")
head(colData)

#counts和olData合并
dds <- DESeqDataSetFromMatrix(mycounts, colData, design= ~ condition)

dds <- DESeq(dds)
dds
sizeFactors(dds)
#慎用，control vs treatment
#res <- results(dds, contrast=c("condition", "control", "treatment"))
# treatment vs control 进行差异分析
res = results(dds)
res <- results(dds, contrast=c("condition", "treatment", "control"))

res <- res[order(res$pvalue),]
head(res)
summary(res)
#查看数据类型
class(res)
#保存
write.csv(res,file="DESEQ2/D6_VS_D9_results.csv")
table(res$padj<0.05)#筛选(res$padj<0.05,na.rm=TRUE)
#筛选想要的列
diff_gene_deseq2 <-subset(res, padj < 0.05 & abs(log2FoldChange) > 1.5)##diff_gene_deseq2 <-subset(res,padj < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1))
dim(diff_gene_deseq2)  
head(diff_gene_deseq2)
write.csv(diff_gene_deseq2,file= "DESEQ2/D6_VS_D9_diffExpression.csv")

####进一步筛选
diff_gene_deseq2 <- na.omit(diff_gene_deseq2)#返回不含NA的值
###取出两列log2FoldChange和padj
nrDEG <- diff_gene_deseq2[,c(2,6)]
View(nrDEG)
###修改列名
colnames(nrDEG) <- c('log2FoldChange','pvalue')
View( nrDEG)

write.csv(nrDEG,file= "DESEQ2/D6_VS_D9__diffExpression_logFC.csv")

结果：得到差异基因的列表，如下，可以根据log2foldchange和padj的值筛选差异表达的基因。一般设置padj<0.05或者0.01, log2foldchange绝对值大于1/1.5/2，具体的根据自己实验进行设置。

image.png