转录组测序的目的之一，就是探索组间显著表达变化的基因——差异表达基因。
那么，如何基于转录组测序获得的定量表达值，识别差异表达变化的基因或其它非编码RNA分子呢？
DESeq2是目前文献中寻找差异基因使用频率最高的R包之一。

DEseq2安装

#若未安装“BiocManager”
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("DESeq2")

#若已安装“BiocManager”
library(BiocManager)
BiocManager::install("DESeq2")

安装完成直接library使用即可。

差异分析

DESeq2是一种基于负二项式分布的方法，使用局部回归推算均值和方差，通过离散度和倍数变化的收缩率估计以提高稳定性。定量分析关注的更多是差异表达的“强度”，而非“存在”。

输入数据

表达量矩阵（我们就是用之前使用RNA-seq数据分析—HTseq定量的结果即可），若有自己的数据要求如下：

• 输入数据是由整数组成的矩阵。

• 矩阵是没有标准化的。

  
  

筛选差异基因

差异分析
DESeq2包分析差异表达基因简单来说只有三步：构建dds矩阵，标准化，以及进行差异分析，我利用以下代码实现：

setwd("F:/RNA-Seq/GSE176393（SRP323246）/5-Diff")

library(DESeq2)

#读入输入数据，设计差异分组
rawdata <- read.table("rawread.txt",header=T,sep = "\t",row.names = 1)
diffcount <- rawdata[,1:6]  ##由于我的数据有三个条件，所以我把需要进行差异分析的两组数据挑选出来，若分析数据仅包含两种条件则可不运行
sampleTable <- data.frame(condition = factor(rep(c("Control", "shBmal1"), each = 3))) ##可以依据实际情况修改各个条件的名称及数量

#构建dds矩阵
dds <- DESeqDataSetFromMatrix(countData = round(diffcount), colData = sampleTable, design = ~condition)
#对原始dds进行normalize
dds <- DESeq(dds)
res <- results(dds)# 将结果用results()函数来获取，赋值给res变量

#保存全部差异结果
diff_res <- as.data.frame(res)
diff_res$gene_id <- rownames(diff_res)
diff_res<-diff_res[, colnames(diff_res)[c(7,1:6)]]
write.table(diff_res,file = 'All_DESeq2_results.xls',sep = '\t',row.names = FALSE)

#依据提取符合阈值的差异结果
table(diff_res$padj<0.05) #取FDR小于0.05的结果，阈值可依据实际需要调整
diff_FDR <- diff_res[order(diff_res$padj),]
diff_gene_deseq2_FDR <-subset(diff_FDR,padj < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1)) #阈值可依据实际需要调整
write.table(diff_gene_deseq2_FDR,file= "diff_gene_deseq2_FDR.xls",sep = '\t',row.names = F)

table(diff_res$pvalue<0.05) #取P值小于0.05的结果
diff_p <- diff_res[order(diff_res$pvalue),]
diff_gene_deseq2_p <-subset(diff_p,pvalue < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1)) #阈值可依据实际需要调整
write.table(diff_gene_deseq2_p,file= "diff_gene_deseq2_pvalue.xls",sep = '\t',row.names = F)

可视化

差异结果的可视化以火山图形式呈现：

#绘制火山图(FDR筛选)
library(ggplot2)

cut_off_stat = 0.05 #设置统计值阈值
cut_off_logFC = 1 #设置表达量阈值,可修改
diff_res[is.na(diff_res)] <- 0
diff_res$change = ifelse(diff_res$padj < cut_off_stat & abs(diff_res$log2FoldChange) >= cut_off_logFC, 
                        ifelse(diff_res$log2FoldChange> cut_off_logFC ,'Up','Down'),
                        'Stable') 
pdf("volcanol_FDR.pdf")
ggplot(diff_res, 
  aes(x = log2FoldChange, 
      y = -log10(padj), 
      colour=change)) +
  geom_point(alpha=0.4, size=3.5) +
  scale_color_manual(values=c("#546de5", "#d2dae2","#ff4757"))+
  geom_vline(xintercept=c(-1,1),lty=4,col="black",lwd=0.8) +
  geom_hline(yintercept = -log10(cut_off_stat),lty=4,col="black",lwd=0.8) +
  labs(x="log2(fold change)",
       y="-log10 (FDR)")+
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5), 
        legend.position="right", 
        legend.title = element_blank()
  )
dev.off()

#绘制火山图(pvalue筛选)

library(ggplot2)

cut_off_stat = 0.05 #设置统计值阈值
cut_off_logFC = 1 #设置表达量阈值,可修改
diff_res[is.na(diff_res)] <- 0
diff_res$change = ifelse(diff_res$pvalue < cut_off_stat & abs(diff_res$log2FoldChange) >= cut_off_logFC, 
                         ifelse(diff_res$log2FoldChange> cut_off_logFC ,'Up','Down'),
                         'Stable')
pdf("volcanol_pvalue.pdf")
ggplot(diff_res, 
       aes(x = log2FoldChange, 
           y = -log10(pvalue), 
           colour=change)) +
  geom_point(alpha=0.4, size=3.5) +
  scale_color_manual(values=c("#546de5", "#d2dae2","#ff4757"))+
  geom_vline(xintercept=c(-1,1),lty=4,col="black",lwd=0.8) +
  geom_hline(yintercept = -log10(cut_off_stat),lty=4,col="black",lwd=0.8) +
  labs(x="log2(fold change)",
       y="-log10 (p-value)")+
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5), 
        legend.position="right", 
        legend.title = element_blank()
  )
dev.off()

##绘制带有基因名称的火山图
library(ggplot2) 
library(ggrepel)

cut_off_stat = 0.05 #设置统计值阈值
cut_off_logFC = 1 #设置表达量阈值,可修改
diff_res[is.na(diff_res)] <- 0
diff_res$change = ifelse(diff_res$pvalue < cut_off_stat & abs(diff_res$log2FoldChange) >= cut_off_logFC, 
                         ifelse(diff_res$log2FoldChange> cut_off_logFC ,'Up','Down'),
                         'Stable')
diff_res$label = ifelse(diff_res$padj < cut_off_stat & abs(diff_res$log2FoldChange) >= 4, as.character(diff_res$gene_id),"")
pdf("volcanol_gene.pdf")
ggplot(diff_res, 
       aes(x = log2FoldChange, 
           y = -log10(padj), 
           colour=change)) +
  geom_point(alpha=0.4, size=3.5) +
  scale_color_manual(values=c("#546de5", "#d2dae2","#ff4757"))+
  geom_vline(xintercept=c(-1,1),lty=4,col="black",lwd=0.8) +
  geom_hline(yintercept = -log10(cut_off_stat),lty=4,col="black",lwd=0.8) +
  labs(x="log2(fold change)",
       y="-log10 (FDR)")+
  theme_bw()+
  theme(plot.title = element_text(hjust = 0.5), 
        legend.position="right", 
        legend.title = element_blank()
  )+
  geom_text_repel(data = diff_res, aes(x = diff_res$log2FoldChange, 
                                      y = -log10(diff_res$padj), 
                                      label = label),
                  size = 3,box.padding = unit(0.8, "lines"),
                  point.padding = unit(0.8, "lines"), 
                  show.legend = FALSE)
dev.off()