轉(zhuǎn)錄組測序的目的之一，就是探索組間顯著表達(dá)變化的基因——差異表達(dá)基因。
那么，如何基于轉(zhuǎn)錄組測序獲得的定量表達(dá)值，識(shí)別差異表達(dá)變化的基因或其它非編碼RNA分子呢？
DESeq2是目前文獻(xiàn)中尋找差異基因使用頻率最高的R包之一。

DEseq2安裝

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

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

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

差異分析

DESeq2是一種基于負(fù)二項(xiàng)式分布的方法，使用局部回歸推算均值和方差，通過離散度和倍數(shù)變化的收縮率估計(jì)以提高穩(wěn)定性。定量分析關(guān)注的更多是差異表達(dá)的“強(qiáng)度”，而非“存在”。

輸入數(shù)據(jù)

表達(dá)量矩陣（我們就是用之前使用RNA-seq數(shù)據(jù)分析—HTseq定量的結(jié)果即可），若有自己的數(shù)據(jù)要求如下：

• 輸入數(shù)據(jù)是由整數(shù)組成的矩陣。

• 矩陣是沒有標(biāo)準(zhǔn)化的。

  
  

篩選差異基因

差異分析
DESeq2包分析差異表達(dá)基因簡單來說只有三步：構(gòu)建dds矩陣，標(biāo)準(zhǔn)化，以及進(jìn)行差異分析，我利用以下代碼實(shí)現(xiàn)：

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

library(DESeq2)

#讀入輸入數(shù)據(jù)，設(shè)計(jì)差異分組
rawdata <- read.table("rawread.txt",header=T,sep = "\t",row.names = 1)
diffcount <- rawdata[,1:6]  ##由于我的數(shù)據(jù)有三個(gè)條件，所以我把需要進(jìn)行差異分析的兩組數(shù)據(jù)挑選出來，若分析數(shù)據(jù)僅包含兩種條件則可不運(yùn)行
sampleTable <- data.frame(condition = factor(rep(c("Control", "shBmal1"), each = 3))) ##可以依據(jù)實(shí)際情況修改各個(gè)條件的名稱及數(shù)量

#構(gòu)建dds矩陣
dds <- DESeqDataSetFromMatrix(countData = round(diffcount), colData = sampleTable, design = ~condition)
#對(duì)原始dds進(jìn)行normalize
dds <- DESeq(dds)
res <- results(dds)# 將結(jié)果用results()函數(shù)來獲取，賦值給res變量

#保存全部差異結(jié)果
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)

#依據(jù)提取符合閾值的差異結(jié)果
table(diff_res$padj<0.05) #取FDR小于0.05的結(jié)果，閾值可依據(jù)實(shí)際需要調(diào)整
diff_FDR <- diff_res[order(diff_res$padj),]
diff_gene_deseq2_FDR <-subset(diff_FDR,padj < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1)) #閾值可依據(jù)實(shí)際需要調(diào)整
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的結(jié)果
diff_p <- diff_res[order(diff_res$pvalue),]
diff_gene_deseq2_p <-subset(diff_p,pvalue < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1)) #閾值可依據(jù)實(shí)際需要調(diào)整
write.table(diff_gene_deseq2_p,file= "diff_gene_deseq2_pvalue.xls",sep = '\t',row.names = F)

可視化

差異結(jié)果的可視化以火山圖形式呈現(xiàn)：

#繪制火山圖(FDR篩選)
library(ggplot2)

cut_off_stat = 0.05 #設(shè)置統(tǒng)計(jì)值閾值
cut_off_logFC = 1 #設(shè)置表達(dá)量閾值,可修改
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 #設(shè)置統(tǒng)計(jì)值閾值
cut_off_logFC = 1 #設(shè)置表達(dá)量閾值,可修改
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 #設(shè)置統(tǒng)計(jì)值閾值
cut_off_logFC = 1 #設(shè)置表達(dá)量閾值,可修改
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()