DESeq2是一種常用的差異表達(dá)基因分析工具，可用于RNA-seq數(shù)據(jù)的差異表達(dá)分析。下面是DESeq2的詳細(xì)使用步驟和全部腳本示例。
文章參考
Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2 | Genome Biology | Full Text (biomedcentral.com)
bioconda源對(duì)工具包的介紹：
Bioconductor - DESeq2
安裝
下面是在R中安裝DESeq2的詳細(xì)步驟：
-
安裝R和RStudio：
- 如果你還沒(méi)有安裝R，可以在R官方網(wǎng)站下載并安裝最新版本的R。
- 推薦使用RStudio作為R語(yǔ)言的集成開發(fā)環(huán)境。你可以在RStudio官網(wǎng)下載并安裝適合你操作系統(tǒng)的版本。
-
啟動(dòng)R或RStudio：
- 打開R或者RStudio。
-
安裝DESeq2包：
- 在R或RStudio的命令行中輸入以下命令安裝DESeq2包：
if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("DESeq2")這將會(huì)從Bioconductor倉(cāng)庫(kù)安裝DESeq2包及其依賴項(xiàng)。
-
加載DESeq2包：
- 安裝完成后，在R或RStudio中輸入以下命令加載DESeq2包：
library(DESeq2)確保沒(méi)有報(bào)錯(cuò)，這樣DESeq2包就已經(jīng)成功加載了。
安裝完成后，你就可以使用DESeq2進(jìn)行基因表達(dá)差異分析了。記得在分析之前準(zhǔn)備好你的RNA-seq數(shù)據(jù)并按照DESeq2的文檔或教程進(jìn)行分析。
DESeq2 的完整使用步驟和示例分析腳本，以及每個(gè)步驟的輸入、輸出和解釋。
步驟 1: 讀取和整理數(shù)據(jù)
首先，加載必要的 R 包和數(shù)據(jù)文件。數(shù)據(jù)應(yīng)該包括表達(dá)矩陣和樣本信息。
# 讀取 DESeq2 包
library(DESeq2)
# 讀取表達(dá)矩陣
countData <- as.matrix(read.csv("count_matrix.csv", row.names = 1))
# 讀取樣本信息
sampleInfo <- read.csv("sample_info.csv")
# 創(chuàng)建 DESeq2 數(shù)據(jù)對(duì)象
dds <- DESeqDataSetFromMatrix(countData, colData = sampleInfo, design = ~ condition)
-
count_matrix.csv: 包含基因或轉(zhuǎn)錄本的表達(dá)矩陣，行代表基因或轉(zhuǎn)錄本，列代表樣本。 -
sample_info.csv: 包含每個(gè)樣本的信息，例如條件、分組等。
步驟 2: 數(shù)據(jù)標(biāo)準(zhǔn)化和差異表達(dá)分析
使用 DESeq2 對(duì)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化和差異表達(dá)分析。
# 標(biāo)準(zhǔn)化數(shù)據(jù)
dds <- DESeq(dds)
# 進(jìn)行差異表達(dá)分析
res <- results(dds)
步驟 3: 結(jié)果解釋和可視化
對(duì)差異表達(dá)結(jié)果進(jìn)行解釋和可視化。
# 查看差異表達(dá)基因
topGenes <- head(rownames(res[order(res$padj), ]), 10)
# 輸出差異表達(dá)基因
write.csv(res, file = "deseq2_results.csv")
# 繪制差異表達(dá)圖
plotCounts(dds, gene = topGenes, intgroup = "condition")
-
deseq2_results.csv: 包含差異表達(dá)分析結(jié)果的輸出文件。
解釋和注意事項(xiàng)：
-
DESeqDataSetFromMatrix(): 用于創(chuàng)建 DESeq2 數(shù)據(jù)對(duì)象，其中countData是表達(dá)矩陣，sampleInfo包含樣本信息，design參數(shù)指定實(shí)驗(yàn)設(shè)計(jì)。 -
DESeq(): 對(duì)數(shù)據(jù)進(jìn)行歸一化和標(biāo)準(zhǔn)化，準(zhǔn)備進(jìn)行差異表達(dá)分析。 -
results(): 提取差異表達(dá)分析的結(jié)果，包括基因表達(dá)差異統(tǒng)計(jì)信息。 - 結(jié)果包括基因表達(dá)水平的差異統(tǒng)計(jì)指標(biāo)，如 fold change、調(diào)整的 p 值（padj）等。
-
plotCounts(): 用于繪制基因表達(dá)水平的差異示意圖，以更直觀地展示不同條件下基因的表達(dá)情況。
使用案例
以下是三個(gè)使用 DESeq2 工具包的案例，包括完整的腳本以及輸入輸出文件內(nèi)容和格式的詳細(xì)解釋。
案例 1: 基因差異表達(dá)分析
輸入文件:
-
count_matrix.csv: 包含基因表達(dá)計(jì)數(shù)矩陣，行代表基因，列代表樣本。 -
sample_info.csv: 包含每個(gè)樣本的信息，例如條件或組別。
腳本:
# 讀取 DESeq2 包
library(DESeq2)
# 讀取表達(dá)矩陣和樣本信息
countData <- as.matrix(read.csv("count_matrix.csv", row.names = 1))
sampleInfo <- read.csv("sample_info.csv")
# 創(chuàng)建 DESeq2 數(shù)據(jù)對(duì)象
dds <- DESeqDataSetFromMatrix(countData, colData = sampleInfo, design = ~ condition)
# 標(biāo)準(zhǔn)化數(shù)據(jù)和進(jìn)行差異表達(dá)分析
dds <- DESeq(dds)
res <- results(dds)
# 輸出差異表達(dá)基因列表和統(tǒng)計(jì)信息
write.csv(res, file = "deseq2_results.csv")
# 繪制差異表達(dá)基因的表達(dá)圖
topGenes <- head(rownames(res[order(res$padj), ]), 10)
plotCounts(dds, gene = topGenes, intgroup = "condition")
輸出文件:
-
deseq2_results.csv: 包含差異表達(dá)分析結(jié)果的輸出文件。包括基因、fold change、p 值、調(diào)整的 p 值等信息。 - 圖形文件：包含差異表達(dá)基因的表達(dá)圖，顯示不同條件下基因的表達(dá)情況。
案例 2: 多組實(shí)驗(yàn)設(shè)計(jì)的差異分析
輸入文件:
count_matrix.csv-
sample_info_multigroup.csv: 包含多組實(shí)驗(yàn)設(shè)計(jì)的樣本信息。
腳本:
# 讀取 DESeq2 包
library(DESeq2)
# 讀取表達(dá)矩陣和樣本信息
countData <- as.matrix(read.csv("count_matrix.csv", row.names = 1))
sampleInfo <- read.csv("sample_info_multigroup.csv")
# 創(chuàng)建 DESeq2 數(shù)據(jù)對(duì)象（多組實(shí)驗(yàn)設(shè)計(jì)）
dds <- DESeqDataSetFromMatrix(countData, colData = sampleInfo, design = ~ group + condition)
# 標(biāo)準(zhǔn)化數(shù)據(jù)和進(jìn)行差異表達(dá)分析
dds <- DESeq(dds)
res <- results(dds)
# 輸出差異表達(dá)基因列表和統(tǒng)計(jì)信息
write.csv(res, file = "deseq2_results_multigroup.csv")
輸出文件:
-
deseq2_results_multigroup.csv: 包含多組實(shí)驗(yàn)設(shè)計(jì)差異表達(dá)分析結(jié)果的輸出文件。
| gene_id | baseMean | log2FoldChange | lfcSE | stat | pvalue | padj |
| GeneA | 100 | 1.5 | 0.2 | 7.2 | 0.0001 | 0.001 |
| GeneB | 80 | -0.8 | 0.3 | -4.5 | 0.0002 | 0.002 |
| GeneC | 50 | 2.1 | 0.5 | 6 | 0.0003 | 0.003 |
-
gene_id: 基因或轉(zhuǎn)錄本的標(biāo)識(shí)符。 -
baseMean: 平均表達(dá)量。 -
log2FoldChange: 對(duì)數(shù)變換后的 fold change，表示在不同條件之間的表達(dá)倍數(shù)變化。 -
lfcSE: log2 fold change 的標(biāo)準(zhǔn)誤差。 -
stat: 統(tǒng)計(jì)檢驗(yàn)值。 -
pvalue: 未經(jīng)校正的 p 值。 -
padj: 經(jīng)過(guò)多重假設(shè)檢驗(yàn)校正后的調(diào)整 p 值（通常使用 FDR 校正），用于控制假陽(yáng)性發(fā)現(xiàn)率。
案例 3: 時(shí)間序列分析
輸入文件:
-
count_matrix_timeseries.csv: 包含時(shí)間序列實(shí)驗(yàn)的基因表達(dá)計(jì)數(shù)矩陣。
| GeneID | Sample1 | Sample2 | Sample3 | Sample4 |
| GeneA | 10 | 15 | 20 | 25 |
| GeneB | 5 | 8 | 12 | 18 |
| GeneC | 30 | 35 | 40 | 45 |
| ... |
-
sample_info_timeseries.csv: 包含時(shí)間序列實(shí)驗(yàn)的樣本信息，包括時(shí)間點(diǎn)等信息。
| Sample | TimePoint | Treatment |
| Sample1 | 0 | Control |
| Sample2 | 3 | DrugA |
| Sample3 | 6 | DrugA |
| Sample4 | 9 | Control |
| ... |
腳本:
# 讀取 DESeq2 包
library(DESeq2)
# 讀取表達(dá)矩陣和樣本信息
countData <- as.matrix(read.csv("count_matrix_timeseries.csv", row.names = 1))
sampleInfo <- read.csv("sample_info_timeseries.csv")
# 創(chuàng)建 DESeq2 數(shù)據(jù)對(duì)象（時(shí)間序列實(shí)驗(yàn)設(shè)計(jì)）
dds <- DESeqDataSetFromMatrix(countData, colData = sampleInfo, design = ~ time_point)
# 標(biāo)準(zhǔn)化數(shù)據(jù)和進(jìn)行差異表達(dá)分析
dds <- DESeq(dds)
res <- results(dds)
# 輸出差異表達(dá)基因列表和統(tǒng)計(jì)信息
write.csv(res, file = "deseq2_results_timeseries.csv")
輸出文件:
-
deseq2_results_timeseries.csv: 包含時(shí)間序列實(shí)驗(yàn)差異表達(dá)分析結(jié)果的輸出文件。
?DESeq2分析結(jié)果進(jìn)行差異表達(dá)火山圖的繪制
DESeq2 的結(jié)果文件 deseq2_results.csv，文件的格式類似于：
| gene_id | baseMean | log2FoldChange | lfcSE | stat | pvalue | padj |
| GeneA | 100 | 1.5 | 0.2 | 7.2 | 0.0001 | 0.001 |
| GeneB | 80 | -0.8 | 0.3 | -4.5 | 0.0002 | 0.002 |
| GeneC | 50 | 2.1 | 0.5 | 6 | 0.0003 | 0.003 |
使用 R 語(yǔ)言和 ggplot2 庫(kù)來(lái)繪制火山圖的示例代碼：文章來(lái)源：http://www.zghlxwxcb.cn/news/detail-852321.html
# 導(dǎo)入必要的庫(kù)
library(ggplot2)
library(dplyr) # 用于數(shù)據(jù)處理
# 讀取差異表達(dá)結(jié)果文件
results <- read.csv("deseq2_results.csv")
# 設(shè)定顯著性閾值（例如，調(diào)整的 p 值）
threshold <- 0.05
# 根據(jù)顯著性閾值篩選差異表達(dá)基因
significant_genes <- results %>%
filter(padj < threshold)
# 繪制火山圖
ggplot(results, aes(x = log2FoldChange, y = -log10(padj), color = ifelse(padj < threshold, "red", "black"))) +
geom_point(size = 1.5) +
scale_color_manual(values = c("black", "red"), labels = c("Not significant", "Significant")) +
labs(x = "log2(Fold Change)", y = "-log10(adjusted p-value)", title = "Volcano Plot of Differential Expression") +
theme_minimal()
這段代碼將根據(jù)調(diào)整的 p 值（padj）和 fold change（log2FoldChange）繪制火山圖。顯著性基因?qū)⒁约t色標(biāo)記，非顯著性基因?qū)⒁院谏珮?biāo)記。您可以調(diào)整 threshold 的值來(lái)控制顯著性。文章來(lái)源地址http://www.zghlxwxcb.cn/news/detail-852321.html
到了這里，關(guān)于基因表達(dá)差異分析R工具包DESeq2的詳細(xì)使用方法和使用案例的文章就介紹完了。如果您還想了解更多內(nèi)容，請(qǐng)?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！
