本文为群中小伙伴进行的一次差异分析探索的记录。
前段时间拿到一个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有很大差别。
