hello，今天小云教大家如何使用DESeq2包對(duì)counts表達(dá)矩陣進(jìn)行差異基因分析，并對(duì)分析結(jié)果繪制火山圖實(shí)現(xiàn)差異基因分析的可視化，如果你感興趣的話就和小云一起看下去吧！

首先小云先來簡單介紹一下這個(gè)包：

了解DESeq2

DESeq2是一個(gè)R語言中常用的差異表達(dá)分析包，用于分析高通量測序數(shù)據(jù)中的差異表達(dá)基因。它是Bioconductor項(xiàng)目的一部分，提供了一套統(tǒng)計(jì)方法和功能，用于檢測基因在不同條件下的表達(dá)差異，并找出具有生物學(xué)意義的差異表達(dá)基因。DESeq2包的設(shè)計(jì)目標(biāo)是針對(duì)RNA-seq數(shù)據(jù)進(jìn)行差異表達(dá)分析，并且能夠處理具有低復(fù)制數(shù)、高變異性和有偏差的數(shù)據(jù)。它采用了負(fù)二項(xiàng)分布模型和廣義線性模型（GLM）的方法，通過考慮樣本之間的技術(shù)和生物學(xué)變異，進(jìn)行準(zhǔn)確的差異表達(dá)分析。

DESeq2包提供了一系列功能，包括：

1.?數(shù)據(jù)導(dǎo)入：支持從不同的數(shù)據(jù)源導(dǎo)入RNA-seq數(shù)據(jù)，包括原始計(jì)數(shù)數(shù)據(jù)和已經(jīng)進(jìn)行了基因注釋的數(shù)據(jù)。

2.?數(shù)據(jù)預(yù)處理：包括數(shù)據(jù)規(guī)范化、去除低表達(dá)基因、去除批次效應(yīng)等。

3.?差異表達(dá)分析：使用負(fù)二項(xiàng)分布模型和GLM進(jìn)行統(tǒng)計(jì)分析，計(jì)算差異表達(dá)基因的顯著性。

4.?多組實(shí)驗(yàn)設(shè)計(jì)：支持多組條件下的差異表達(dá)分析，包括兩組比較、多組比較和時(shí)間序列分析。

5.?基因注釋和功能富集分析：支持將差異表達(dá)基因進(jìn)行注釋，并進(jìn)行功能富集分析，以了解差異表達(dá)基因的生物學(xué)功能。

6.?可視化：DESeq2包集成了ggplot2和其他可視化工具，可以方便地繪制差異表達(dá)基因的圖表，如火山圖、MA圖等。

差異基因分析流程

R包與數(shù)據(jù)加載

首先，設(shè)置工作目錄為數(shù)據(jù)文件所在的路徑，并加載DESeq2和其他所需的R包，并導(dǎo)入我們準(zhǔn)備好的counts表達(dá)矩陣：

setwd("/media/desk16/iyun007/kegg")
#BiocManager::install("DESeq2") ?
library(DESeq2)
data1 <- read.csv("Gene Symbol_matrix.csv",header = T, row.names = 1)
dim(data1) ?#59427 ? ?5 ? （COUNT矩陣有59427 個(gè)基因、5個(gè)樣本）

注意，我們導(dǎo)入的數(shù)據(jù)基本格式如下，同學(xué)們要提供符合要求的數(shù)據(jù)集哦！

數(shù)據(jù)預(yù)處理

接下來，我們要構(gòu)建分組矩陣，即根據(jù)數(shù)據(jù)框中樣本的命名規(guī)則，提取出癌癥樣本和正常樣本的列名，并將它們分別存儲(chǔ)在tumor和normal向量中：

#構(gòu)建分組矩陣
tumor <- colnames(data1)[as.integer(substr(colnames(data1),14,15)) < 10] #data1表達(dá)矩陣中的所有癌癥樣本（8個(gè)）
tumor
normal <- colnames(data1)[as.integer(substr(colnames(data1),14,15)) >= 10] #data1表達(dá)矩陣中的所有正常樣本（2個(gè)）
normal

tumor_sample <- data1[,tumor] ?#癌癥樣本的表達(dá)矩陣
View(tumor_sample)
tumor_sample
ncol(tumor_sample)

normal_sample <- data1[,normal] ? #正常樣本的表達(dá)矩陣
View(normal_sample)
ncol(normal_sample)
data1 <- cbind(tumor_sample,normal_sample) ?#合并癌癥和正常樣本矩陣
col_index <- which(colnames(matrix) == "normal_sample")
colnames(matrix)[col_index] <- normal
#colnames(data1)[5] <- "TCGA.AA.3525.11A.01R.A32Z.07"
View(data1)
group_list <- c(rep('tumor',ncol(tumor_sample)),rep('normal',1)) #根據(jù)自己的數(shù)據(jù)集設(shè)置tumor和normal數(shù)據(jù)的分組數(shù)

condition = factor(group_list)
coldata <- data.frame(row.names = colnames(data1), condition) ? #coldata是分組矩陣
View(coldata)

最后合并好的數(shù)據(jù)如下：

分組信息如下：

差異分析

接下來，我們可以：

1.?創(chuàng)建DESeq2數(shù)據(jù)對(duì)象：使用DESeqDataSetFromMatrix()函數(shù)創(chuàng)建DESeq2所需的數(shù)據(jù)對(duì)象，將合并后的數(shù)據(jù)框和分組矩陣作為參數(shù)傳遞給該函數(shù)。

2.?進(jìn)行差異分析：使用DESeq()函數(shù)對(duì)DESeq2數(shù)據(jù)對(duì)象進(jìn)行差異分析。

3.?保存差異表達(dá)基因結(jié)果：將差異表達(dá)基因結(jié)果保存為CSV文件。

#構(gòu)建dds對(duì)象，構(gòu)建差異基因分析所需的數(shù)據(jù)格式
dds <- DESeqDataSetFromMatrix(countData = data1, colData = coldata, design = ~ condition)
#差異分析
dds <- DESeq(dds)
#提取結(jié)果
result <- as.data.frame(results(dds)) ? ? ? ?# results()從DESeq分析中提取出結(jié)果表

DGE <-subset(result,padj < 0.05 & (log2FoldChange > 2 | log2FoldChange < -2)) ? ?#這里我們提取出2倍差異表達(dá)、校正p值<0.05的顯著差異表達(dá)基因
DGE <- DGE[order(DGE$log2FoldChange),]

write.csv(DGE,'DGE.csv',row.names = TRUE)

我們一起來看下分析后的結(jié)果把！

繪制火山圖

接下來，我們可以通過火山圖可視化我們的分析結(jié)果，具體操作如下：

這里小云使用ggplot()函數(shù)創(chuàng)建繪圖對(duì)象，并設(shè)置數(shù)據(jù)源為result數(shù)據(jù)框。然后使用geom_point()函數(shù)繪制散點(diǎn)圖，其中x軸為log2FoldChange，y軸為-log10(padj)，并根據(jù)差異表達(dá)的方向和顯著性進(jìn)行顏色分類。通過調(diào)整繪圖的各種參數(shù)，如標(biāo)題、坐標(biāo)軸標(biāo)簽和圖例，可以使火山圖更具可讀性和美觀性哦。

## 5.提取結(jié)果
result <- as.data.frame(results(dds)) # results()從DESeq分析中提取出結(jié)果表

result <- na.omit(result)
result1 <- result[c(2,6)]

## 6.火山圖
result1$change = ifelse(result1$padj < 0.05 & abs(result1$log2FoldChange) >= 2, ifelse(result1$log2FoldChange> 2 ,'Up','Down'),'Stable')
library('ggplot2')
ggplot(result1, aes(x = log2FoldChange, y = -log10(padj),colour=change)) +
??geom_point(alpha=0.3, size=3.5) +
??scale_color_manual(values=c("#546de5","#d2dae2","#ff4757"))+
??labs(title='Volcano plot',x="log2(fold change)",y="-log10 (padj)")+ # 坐標(biāo)軸# 坐標(biāo)軸和圖標(biāo)題title="Volcano plot",
??theme_bw()+ #去除背景色
??#xlim(-20, 20)+ #設(shè)置坐標(biāo)軸范圍
??#ylim(0,40)+
??# 輔助線
??geom_vline(xintercept=c(-2,2),lty=3,col="black",lwd=0.8) +
??geom_hline(yintercept = -log10(0.05),lty=3,col="black",lwd=0.8) +
??theme(panel.grid = element_blank())+ #去除網(wǎng)格線
??theme(plot.title = element_text(hjust = 0.5,size=24),
?? ? ? ?legend.position="bottom",
?? ? ? ?legend.title = element_blank(),
?? ? ? ?legend.text=element_text(size=18),
?? ? ? ?legend.key.size = unit(1, 'cm'),
?? ? ? ?legend.background = element_rect(fill="gray90", linetype="solid",colour ="gray"),
?? ? ? ?axis.title.x =element_text(size=18),
?? ? ? ?axis.title.y=element_text(size=18),
?? ? ? ?axis.text=element_text(size=14,face = "bold")
??)
ggsave("De_ana.png",width=8,height=6)# 保存為png格式

我們一起來看下最后的結(jié)果圖吧！

怎么樣，你學(xué)會(huì)怎么使用DESeq2包來進(jìn)行差異分析了嘛？是不是很簡單！更多學(xué)習(xí)干活要繼續(xù)關(guān)注小云哦！