差异分析不是这样做的……

生信菜鸟团

发布于 2023-01-05 13:04:03

2.3K00

代码可运行

文章被收录于专栏：生信菜鸟团生信菜鸟团

运行总次数：0

代码可运行

那天在Frontiers in Immunology 偶然一瞥，看到一篇新鲜出炉的纯生信文章

大喜？！纯生信文章还能往这里投？赶紧学习学习，然后……我就看到了这张神奇的图⬇

原文对差异分析是这么描述的：Using R software’s limma package , differentially expressed genes (DEGs) between septic shock cohort and control cohort were analyzed, with the following criterion: p<0.05 and |fold change (FC)| > 1.

竟然是直接用fold change来作为阈值的，一般差异分析用的更多的是log2FC,这样数字不会特别离谱。

毕竟，log2FC中的FC即 fold change，表示两个样本/组间表达量的比值，对其取以2为底的对数之后才是log2FC。

举几个例子：基因A 在肿瘤和正常组织中的表达量分别为2和4，那么比值为2，即FC=2，此时log2FC=1；基因B 在肿瘤和正常组织中的表达量均为2，此时比值为1，即FC=1，此时log2FC=0；基因C 在肿瘤和正常组织中的表达量分别为2和64，那么比值为32，即FC=32，此时log2FC=5。最后夸张一下！！基因D 在肿瘤和正常组织中的表达量分别为2和1024，那么比值为512，即FC=512，此时log2FC=9。

这样一算，你大概就能明白上面那张图问题出在哪里了吧~

但是光说不练，纸上谈兵，我们还是自己上手分析一下这个数据集，验证一下自己的猜想——

这里的上下调基因取得是top30的哈~,代码放在下面：

cg = rownames(deg)[deg$change != "stable"]
deg1 <- deg[cg,] 
library(dplyr)
top_up <- deg1 %>%  top_n(n=30,wt=logFC) %>% rownames() 
top_down <- deg1 %>%  top_n(n=-30,wt=logFC) %>% rownames()
top_all <- c(top_up,top_down)
n1 <- n[top_all,]

# 由于 ComplexHeatmap 包并不会对数据进行标准化，为了让图形更好看，我们先手动对数据进行标准化
exp <- apply(n1, 1, scale)
rownames(exp) <- colnames(n1)
exp <- t(exp)

然后是画图参数的设置：

col_fun = colorRamp2(c(-2, 0, 2), c("#2fa1dd", "white", "#f87669"))
top_annotation = HeatmapAnnotation(
  cluster = anno_block(gp = gpar(fill = c("#2fa1dd", "#f87669")),
                       labels = c("normal","case"),
                       labels_gp = gpar(col = "white", fontsize = 12)))
p4 <- pheatmap(exp,name = "expression",
               col = col_fun,
               column_split = group_list,
               top_annotation=top_annotation,
               row_split = rep(c("UP","DOWN"), each = 30),
               fontsize_row=5,
               column_title = NULL,
               show_colnames = F)
p4                
p4 <- as.ggplot(p4)

最后是拼图的基操：

library(patchwork)
p <- p3+p4 ##p3就是火山图
ggsave(p,filename = "../outputs/article.pdf",width = 15,height = 9)

为什么我的图和文章里的图迥然不同呢？

相信大家只要对表达量矩阵有一定的熟悉，就应该知道，有的数据集下载以后，需要先观察探针在每一个样本中的表达量数据，一般数值不大于20的话，说明这个矩阵已经是被取过log的，否则的话是需要先取log再做分析的。

options(timeout = 99999999)

getwd()
rm(list = ls())
library(GEOquery)
library(stringr)

gse <- "GSE26440"
# 1.获得表达矩阵 ----------------------------------------------------------------
eSet1 <- getGEO(gse, 
                destdir = '.', 
                getGPL = F)
#(1)提取表达矩阵exp
exp1 <- exprs(eSet1[[1]])
exp1[1:4,1:4]
range(exp1)
# [1]  0.01 1208.00
exp1 <- log2(exp1+1)
range(exp1)