差异分析|DESeq2完成配对样本的差异分析

简介: 差异分析|DESeq2完成配对样本的差异分析



本文为群中小伙伴进行的一次差异分析探索的记录。

前段时间拿到一个RNA-seq测序数据(病人的癌和癌旁样本,共5对)及公司做的差异分析结果(1200+差异基因),公司告知用的是配对样本的DESeq分析。

考虑到平时limma和DESeq2包进行差异分析时没有特别注明是否配对,这配对和非配对有啥区别呢?

于是分别尝试使用limma和DESeq2包的非配对分析,发现得到的差异基因和公司的差距很大。我查了好多关于RNA-seq配对分析的资料,发现几乎没有这方面的帖子。

询问公司DESeq配对分析的代码,公司说保密不能给,此外公司还告知现在的配对样本的分析都改用了DESeq2。好吧,那就只能自己动手,以下为探索过程的一个记录。


一 载入R包,数据


# 加载包
library(openxlsx)
library(DESeq2)
library(limma)
library(edgeR)


1.读入原始数据及分组信息

rowdata <- read.xlsx("Count_data.xlsx",sheet = 1,rowNames = T)
gset <- rowdata[rowMeans(rowdata)>0,] # 剔除表达量低的基因
group_list <- c(rep("case",5), rep("control",5))
group_list <- factor(group_list,levels = c("control","case"))




二 limma包进行差异分析

## 2.1表达矩阵
data = gset
## 2.2分组矩阵
design <- model.matrix(~0+group_list)
rownames(design) = colnames(data)
colnames(design) <- levels(group_list)
## 2.3差比表达矩阵构建,并过滤
DGElist <- DGEList( counts = data, group = group_list)
keep <- rowSums(cpm(DGElist) > 0.5 ) >=2
table(keep)
# FALSE 2443 TRUE 15909
DGElist_QC <- DGElist[keep, ,keep.lib.sizes=FALSE]
dim(DGElist)
# 18352 10 过滤前
dim(DGElist_QC)
# 15909 10过滤后
## 2.4归一化基因表达分布
DGElist_norm <- calcNormFactors(DGElist_QC, method = "TMM")
DGElist_norm$samples$norm.factors # 查看每个样品归一化后的系数
DGElist_QC$samples$norm.factors #未归一化之前的样本系数,都是1
## 2.5 limma包进行voom函数
v <- voom(DGElist_norm, design, plot = TRUE, normalize = "quantile")
fit <- lmFit(v, design)
cont.matrix <- makeContrasts(contrasts = c('case-control'), levels = design)
fit2 <- contrasts.fit(fit, cont.matrix)
fit2 <- eBayes(fit2)
## 2.6提取差异矩阵
nrDEG_limma_voom = topTable(fit2, coef = "case-control", n = Inf)
nrDEG_limma_voom = na.omit(nrDEG_limma_voom)
nrDEG_limma_voom = nrDEG_limma_voom[order(nrDEG_limma_voom$logFC),]
## 2.7 定义差异基因
nrDEG <- nrDEG_limma_voom
nrDEG$Group = "notsignificant"
logFC_cutoff <- 0.6 # 定义差异基因的标准,自定义
nrDEG$Group [which( (nrDEG$P.Value < 0.05) & (nrDEG$logFC > logFC_cutoff) )] = "upregulated"
nrDEG$Group [which( (nrDEG$P.Value < 0.05) & (nrDEG$logFC < -logFC_cutoff) )] = "downregulated"
table(nrDEG$Group)

可以看到只有67个下调的33个上调的,火山图不好看,而且根本没法继续做GO和KEGG分析。

OK,尝试使用DESeq2包的非配对差异分析。




三 DESeq2进行差异分析


## 3.1表达矩阵
data = apply(gset, 2, as.integer) ## DESeq2分析需要是整数
row.names(data) <- row.names(gset)
## 3.2分组矩阵
condition = group_list
coldata <- data.frame(row.names = colnames(data), condition)
dds <- DESeqDataSetFromMatrix(countData = data,
                             colData = coldata,
                             design = ~condition)
dds$condition<- relevel(dds$condition, ref = "control") # 指定哪一组作为对照组
## 3.3差异表达矩阵
dds <- DESeq(dds)  
nrDEG_DESeq2 <- as.data.frame(results(dds))
nrDEG_DESeq2 = nrDEG_DESeq2[order(nrDEG_DESeq2$log2FoldChange),]
## 3.4定义差异基因
nrDEG <- nrDEG_DESeq2
nrDEG$Group = "notsignificant"
logFC_cutoff <- 0.6
nrDEG$Group[which( (nrDEG$pvalue < 0.05) & (nrDEG$log2FoldChange > logFC_cutoff) )] = "upregulated"
nrDEG$Group[which( (nrDEG$pvalue < 0.05) & (nrDEG$log2FoldChange < -logFC_cutoff) )] = "downregulated"
table(nrDEG$Group)

可以看到常规的DESeq2分析比limma voom分析多了一些差异基因,但是和公司给的1200+的差异基因还是差远了。


发现差异之后开始了检索和求助之旅,查了很多帖子,也求助了一些大神,似乎很少人注意过DESeq2包做配对的差异分析。

讨论群中小伙伴贴了limma进行配对分析的方式,于是我去查阅了DESeq2的说明书,可以看到说明书就这么一句话:

剩下的事情就简单了,依此修改后,DESeq2包成功做出了配对差异分析,复现了公司的结果。好了,下面就是使用DESeq2包完成配对差异分析的代码了,自取!


四 DESeq2配对差异分析


## 4.1表达矩阵
data = apply(gset, 2, as.integer) ## DESeq2分析需要是整数
row.names(data) <- row.names(gset)
## 4.2分组矩阵,配对分析与常规分析最大的区别就在分组矩阵
condition = group_list
# 配对分析要加上这段代码,知道谁和谁是一对,比如1,1是一对,5,5是一对
subject <- factor(c(1,2,3,4,5,1,2,3,4,5))  
coldata <- data.frame(row.names = colnames(data), condition)
# 注意在design中加上配对信息
dds <- DESeqDataSetFromMatrix(countData = data,
                             colData = coldata,
                             design = ~subject +condition)
dds$condition<- relevel(dds$condition, ref = "control")
## 4.3差异表达矩阵,还是和常规分析一样
dds <- DESeq(dds)
nrDEG_DESeq2 <- as.data.frame(results(dds))
rld <- rlog(dds)
# 这里我还提取了标准化后的表达矩阵,可以用于后续的热图绘制等等
normal_gset <- assay(rld)
nrDEG_DESeq2 = nrDEG_DESeq2[order(nrDEG_DESeq2$log2FoldChange),]
## 4.4定义差异基因
nrDEG <- nrDEG_DESeq2
nrDEG$Group = "notsignificant"
logFC_cutoff <- 0.6
nrDEG$Group[which( (nrDEG$pvalue < 0.05) & (nrDEG$log2FoldChange > logFC_cutoff) )] = "upregulated"
nrDEG$Group[which( (nrDEG$pvalue < 0.05) & (nrDEG$log2FoldChange < -logFC_cutoff) )] = "downregulated"
table(nrDEG$Group)

此时终于成功得到了1200+的差异基因,通过对比公司的分析结果和我做的结果,几乎完成了复现。有一点区别在于我做了一些低表达基因的过滤,导致log2Foldchang和pvalue略微有些变化,但这区别可以忽略不计。

总结来说,由于算法的不同,不同差异分析的R包得到的差异基因数量不完全一致。重要的是,针对配对的样本,如果不进行配对分析而用常规的差异分析,这样的结果可能会大不相同。因此,在分析数据的时候,一定要明白实验设计


最后,我还发现有意思的一个情况。在进行clusterProfilerR包的GSEA分析的时候,我用非配对分析得到的log2FoldChang出的GSEA结果,和配对分析得到的log2FoldChang出的GSEA结果几乎是一致的,尽管两种分析方法得到的log2FoldChang有很大差别。


相关文章
|
28天前
偏微分方程有了基础模型:样本需求数量级减少,14项任务表现最佳
【6月更文挑战第16天】研究人员提出Poseidon模型,减少求解偏微分方程(PDEs)的样本需求,提升效率。在15个挑战任务中,该模型在14项表现最优。基于scOT的多尺度架构, Poseidon降低了计算成本,但仍有泛化和资源限制。[论文链接](https://arxiv.org/pdf/2405.19101)**
60 4
|
2月前
|
机器学习/深度学习 数据可视化
数据分享|R语言生存分析模型因果分析:非参数估计、IP加权风险模型、结构嵌套加速失效(AFT)模型分析流行病学随访研究数据
数据分享|R语言生存分析模型因果分析:非参数估计、IP加权风险模型、结构嵌套加速失效(AFT)模型分析流行病学随访研究数据
|
2月前
R语言估计多元标记的潜过程混合效应模型(lcmm)分析心理测试的认知过程
R语言估计多元标记的潜过程混合效应模型(lcmm)分析心理测试的认知过程
|
2月前
R语言向量误差修正模型 (VECMs)分析长期利率和通胀率影响关系
R语言向量误差修正模型 (VECMs)分析长期利率和通胀率影响关系
|
2月前
【SPSS】两独立样本的极端反应检验和两配对样本的非参数检验详细操作教程(附案例实战)
【SPSS】两独立样本的极端反应检验和两配对样本的非参数检验详细操作教程(附案例实战)
94 0
|
2月前
|
资源调度 算法 数据挖掘
变异系数法:一种强大的数据离散度度量工具
变异系数法:一种强大的数据离散度度量工具
119 0
变异系数法:一种强大的数据离散度度量工具
|
11月前
|
JSON 算法 数据格式
优化cv2.findContours()函数提取的目标边界点,使语义分割进行远监督辅助标注
可以看到cv2.findContours()函数可以将目标的所有边界点都进行导出来,但是他的点存在一个问题,太过密集,如果我们想将语义分割的结果重新导出成labelme格式的json文件进行修正时,这就会存在点太密集没有办法进行修改,这里展示一个示例:没有对导出的结果进行修正,在labelme中的效果图。
118 0
|
数据可视化 Serverless Go
scRNA分析|单细胞GSVA + limma差异分析-celltype分组?样本分组?
scRNA分析|单细胞GSVA + limma差异分析-celltype分组?样本分组?
811 0
|
Linux Windows Perl
没有生物学重复的转录组数据怎么进行差异分析?
设置生物学重复这个环节也是你实验设计很重要的一part,设置的好对你下游分析也有利,通常我们做转录组测序,需要的样本量每组至少为3个生物学重复,这个处理起来就很合理,并且现在流行的差异分析软件DEseq2,limma,edgeR等等都是针对有重复的数据去做的,但有时候会不幸碰到样品测序失败不能用,导致每组就给你剩一个重复时候该怎么办,之前我有批数据就是这样,但是办法总比困难多不能放过任何实验数据,搜了搜其实还是有一些方法可以去解决的,在这里介绍下我搜到的几种方法。
706 0
|
数据挖掘 Serverless
Robust火山图:一种含离群值的代谢组数据差异分析方法
代谢组学中差异代谢物的识别仍然是一个巨大的挑战,并在代谢组学数据分析中发挥着突出的作用。由于分析、实验和生物的模糊性,代谢组学数据集经常包含异常值,但目前可用的差异代谢物识别技术对异常值很敏感。作者这里提出了一种基于权重的具有稳健性火山图方法,助于从含有离群值的代谢组数据中更加准确鉴定差异代谢物。
165 0