百日筑基篇——差异基因分析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 筑基篇

  • 1
    点赞
  • 18
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 4
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 4
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

星石传说

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值