DESeq2 差異分析

author：

小杜

的生信筆記

DEseq2官方網址：

http：//www。

bioconductor。org/packag

es/release/bioc/html/DESeq2。html

DESeq2是最常用的差異分析的方法，在上一個教程我們分享了limma差異分析 | 簡書，差異分析——limma包 - 知乎 等平臺，我們後續也會分享常用的差異分析方法，並進行一個比較。

DESeq2包是為高維計數資料的歸一化、視覺化和

差分分析

而設計的。 它利用經驗貝葉斯技術估計對數摺疊變化和

離散

的先驗，並計算這些量的後驗估計。

1 匯入所需包

#if （！requireNamespace（“BiocManager”， quietly = TRUE））

install。packages（“BiocManager”）

#BiocManager：：install（“DESeq2”）

library（DESeq2）

library（dplyr）

2 匯入Counts

資料矩陣

## 匯入Counts資料矩陣

countdata <- read。table（“Countdata。txt”， header = T，row。names = 1）

countdata［1：10，1：6］

> countdata［1：10，1：6］

Control_01 Control_02 Control_03 Treat_01 Treat_02 Treat_03

gene01 11 11 11 10 11 10

gene02 14 14 14 13 14 13

gene03 14 14 14 13 14 13

gene04 7 7 7 6 7 6

gene05 11 11 11 11 11 10

gene06 14 14 14 13 14 12

gene07 13 13 13 12 13 12

gene08 11 11 10 10 11 9

gene09 11 11 11 10 11 9

gene10 10 12 12 12 12 10

3

資料過濾

## 過濾在所有重複樣本中小於1的基因

countdata = countdata［rowMeans（countdata） > 1，］

4 資料樣本資訊

樣本註釋資訊自己在本地後製作後進行匯入

coldata <- read。csv（“countdata。csv”， header = T， row。names = 1）

head（coldata）

> head（coldata）

condition

Control_01 CK

Control_02 CK

Control_03 CK

Treat_01 Treat

Treat_02 Treat

Treat_03 Treat

檢查資料Counts檔案與coldata資料是否匹配

all（rownames（coldata） %in% colnames（countdata））

all（rownames（coldata） == colnames（countdata））

當返回

TRUE

時，表明兩個數匹配。

> all（rownames（coldata） %in% colnames（countdata）） #

［1］ TRUE

> all（rownames（coldata） == colnames（countdata））

［1］ TRUE

5 差異分析

## 製作差異矩陣

dds <- DESeqDataSetFromMatrix（countData = countdata，colData = coldata，design = ~ condition）

dim（dds）

## 過濾

dds <- dds［rowSums（counts（dds）） > 1，］

nrow（dds）

## 差異比較

dep <- DESeq（dds）

res <- results（dep）

diff = res

diff <- na。omit（diff） ## 去除NA

dim（diff）

write。csv（diff，“all_diff。csv”） # 匯出所有的差異檔案

DESeq輸出的差異檔案

> head（diff）

log2 fold change （MLE）： condition Treat vs CK

Wald test p-value： condition Treat vs CK

DataFrame with 6 rows and 6 columns

baseMean log2FoldChange lfcSE stat pvalue padj

gene01 10。6594 -0。01071206 0。788745 -0。0135811 0。989164 0。999747

gene02 13。6626 0。00959682 0。713528 0。0134498 0。989269 0。999747

gene03 13。6626 0。00959682 0。713528 0。0134498 0。989269 0。999747

gene05 10。8224 0。03339267 0。783215 0。0426354 0。965992 0。999747

gene06 13。4731 -0。03001613 0。716954 -0。0418662 0。966605 0。999747

gene07 12。6616 0。00389414 0。735478 0。0052947 0。995775 0。999747

篩選差異基因，使用P值和LogFC進行篩選，這裡我們就不多做解釋，詳情可以看我們前面的推文。

# 設定篩選標準

foldChange = 1

padj = 0。05

#

diffsig <- diff［（diff$pvalue < padj & abs（diff$log2FoldChange） > foldChange），］

dim（diffsig）

write。csv（diffsig， “All_diffsig。csv”）

注：上調和下調的差異基因的篩選，我們在這裡不在講述，請檢視我們前面的推文即可

6 差異基因

熱圖

6。1 標準化Count值

我們後面的熱圖使用Count值進行分析，以及後續的分析也可以使用count值進行分析。我們的Count值是基因的原始表達量，因此需要進行標準化出來，我們這裡使用

vst

函式進行標準化處理。

## 標準化（標準化Counts值）

vsd <- vst（dds， blind = FALSE）

normalizeExp <- assay（vsd）

繪製差異基因熱圖

### 差異表達熱圖

df <- read。csv（“All_diffsig。csv”， header = T）

head（df）

df02 <- as。character（df$X）

## 標準化（標準化Counts值）

vsd <- vst（dds， blind = FALSE）

normalizeExp <- assay（vsd）

## 差異基因的Count

diff_expr <- normalizeExp［df02，］

head（diff_expr）

## 差異熱圖

library（pheatmap）

annotation_col <- data。frame（Group = factor（c（rep（“CK”，3）， rep（“Treat”，3））））

rownames（annotation_col） <- colnames（diff_expr）

pdf（“差異基因熱圖。pdf”， height = 8， width = 12）

p <- pheatmap（diff_expr，

annotation_col =

annotation_col

，

color = colorRampPalette（c（“green”，“black”，“red”））（50），

cluster_cols

= F，

show_rownames = F，

show_colnames = F，

scale = “row”， ## none， row， column

fontsize = 12，

fontsize_row = 12，

fontsize_col = 6，

border = FALSE）

print（p）

dev。off（）

“小杜的生信筆記” 公眾號、知乎、簡書、B站平臺，主要發表或收錄生物資訊學的教程，以及基於R的分析和視覺化（包括資料分析，圖形繪製等）；分享感興趣的文獻和學習資料！

01、簡書：差異分析|使用limma包 - 簡書 （jianshu。com）

02、知乎：差異分析——limma包 - 知乎 （zhihu。com）

03、公眾號：R語言差異分析 | limma包 （qq。com）