百日筑基篇——差异基因分析DESeq2包（R语言初识六）

百日筑基篇——差异基因分析DESeq2包（R语言初识六）

---

前言

今天，整理一下如何使用R语言进行基因差异分析。主要会讲述有关的几个扩展包以及可视化，希望以这种方式巩固所学。

一、差异基因分析是什么?

差异基因分析是生物信息学中常用的一种分析方法，用于比较不同样本之间基因表达水平的差异。
可以帮助我们识别与特定生物学过程相关的基因，并进一步理解这些基因在不同条件下的调控机制。

DESeq2使用负二项分布模型来估计基因表达的离散性，并通过广义线性模型（GLM）和Wald统计检验来进行差异分析

二. 基本步骤

1. 数据预处理

对原始基因表达数据进行预处理。这包括数据清洗、去除低表达基因和归一化等步骤

代码如下（示例）：

#使用DESeq2包
library(DESeq2)
# 读取原始基因表达数据
gene <- read.table("gene.txt",header = TRUE,sep="\t",row.names = 1)
#剔除全为0值的行    （示例，实际情况具体分析)
all <- apply(gene,1,function(x){all(x==0)})
gene <- gene[!all,]

#注意要保证表达矩阵中的样本顺序和这里的分组顺序是一一对应的
coldata <- data.frame(group = factor(rep(c('control', 'treat'), each = 8)))

2. 创建DESeq2对象

dds <- DESeqDataSetFromMatrix(countData = gene,colData = coldata,design = ~group)

在上述代码中，countData参数用于指定基因表达数据，colData参数用于指定样本信息数据框，design参数用于指定差异分析的设计

3. 差异分析

运行这段代码后，res将包含差异表达基因的结果对象，其中包括调整后的p值、差异倍数的对数（log2FoldChange）等信息。

#进行差异表达分析。
dds <- DESeq(dds, parallel = FALSE)

#提取差异表达基因
res <- results(dds, contrast = c('group', 'treat', 'control'), pAdjustMethod = 'fdr', alpha = 0.05)

dds是之前创建的DESeq2对象；
contrast参数用于指定对比组，其中group是你感兴趣的条件列，treat是你想要比较的实验组，control是你想要比较的对照组；
pAdjustMethod参数用于指定p值校正方法，这里使用了fdr方法；
alpha参数用于设置显著性阈值，这里设置为0.05。

#查看差异表达基因的结果对象
summary(res)
#返回:
out of 55818 with nonzero total read count
adjusted p-value < 0.05
LFC > 0 (up)       : 77, 0.14%
LFC < 0 (down)     : 97, 0.17%
outliers [1]       : 0, 0%
low counts [2]     : 23888, 43%
(mean count < 4)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

4. 绘制MA图观察

MA图显示基因表达水平的平均差异（M值）与表达量的平均数（A值）之间的关系。在MA图中，每个点代表一个基因，x轴表示A值（通常是两个条件的平均表达量的对数比），y轴表示M值（通常是两个条件的差异的对数比）。
差异表达的基因呈现出偏离中心线的趋势或聚集，可以帮助快速发现和评估差异表达的基因。

plotMA(res, alpha = 0.05, ylim = c(-3, 3))

#数据存储
deseq_res <- as.data.frame(res[order(res$padj), ])
deseq_res$gene_id <- rownames(deseq_res)
write.table(deseq_res[c(7, 1:6)], 'DESeq2-test.txt', row.names = FALSE, sep = '\t', quote = FALSE)

---

三、绘制火山图

#ggplot2对差异基因作图
library(ggplot2)
deseq_res <- read.table('DESeq2-test.txt', sep = '\t',header = TRUE)
#读进来最后的差异基因结果并进行分类
#|log2FC| >= 1 & FDR p-value < 0.05 定为差异
deseq_res[which(deseq_res$padj %in% NA),'sig'] <- 'no diff'
deseq_res[which(deseq_res$log2FoldChange >= 1 & deseq_res$padj < 0.05),'sig'] <- 'up (p.adj < 0.05, log2FC >= 1)'
deseq_res[which(deseq_res$log2FoldChange <= -1 & deseq_res$padj < 0.05),'sig'] <- 'down (p.adj < 0.05, log2FC <= -1)'
deseq_res[which(abs(deseq_res$log2FoldChange) < 1 | deseq_res$padj >= 0.05),'sig'] <- 'no diff'

#也可以获取上调up /下调down 的差异表达基因（padjust < 0.05，并且|log2(foldchange)|>1)
diff_up <-  subset(deseq_res,padj < 0.05 & (log2FoldChange > 1))
write.csv(diff_up,file="diff_up.csv",row.names = F)

diff_down <- subset(deseq_res,padj < 0.05 & (log2FoldChange < -1))
write.csv(diff_down,file="diff_down.csv",row.names = F)
#画火山图
volcano_p <- ggplot(deseq_res, aes(log2FoldChange, -log(padj, 10))) +
  geom_point(aes(color = sig), alpha = 0.6, size = 1) +
  scale_color_manual(values = c('blue2', 'gray30', 'red2')) +
  theme(panel.grid = element_blank(),
        panel.background = element_rect(color = 'black', fill = 'transparent'), 
        legend.position = c(0.26, 0.92)) +
  theme(legend.title = element_blank(), 
        legend.key = element_rect(fill = 'transparent'), 
        legend.background = element_rect(fill = 'transparent')) +
  geom_vline(xintercept = c(-1, 1), color = 'gray', linewidth = 0.25) +
  geom_hline(yintercept = -log(0.05, 10), color = 'gray', linewidth = 0.25) +
  labs(x = 'log2 Fold Change', y = '-log10 p-value', color = NA) +
  xlim(-5, 5)

volcano_p
ggsave('volcano_p.png', volcano_p, width = 5, height = 6)

可以根据需要自主挑选数据绘不同图。

总结

基因差异分析的学习，今天就到这里，理解这背后的统计学原理是十分必要的，同时要一开始就确定想要通过基因差异分析，分析出什么，目标明确。
醉里挑灯看剑，梦回吹角连营。

-2023-7-8 筑基篇