再小的差异也能被gsea出来

再小的差异也能被gsea出来

我们的马拉松授课专注于表达量矩阵的数据处理技巧传授，包括表达量芯片，转录组 测序， 单细胞转录组，都是一脉相承的。

每个知识点都有对应的练习题安排给学员来考验大家是否掌握差异分析和富集分析的精髓，其中表达量芯片环节大家完成作业还是比较积极的，转录组测序也还行， 有一半的小伙伴还是可以交作业的！比如其中一个学员就读了这个文章：《KIAA1522 is a new biomarker of promoting the tumorigenesis and distant metastasis of colorectal carcinoma》，对应的数据集是：GSE174449，做了处理，然后提出来了一个很有意思的问题， 就是文献里面的差异基因和通路，在可视化里面非常明显的差异，但实际上处理这个GSE174449时候看到的变化倍数还有gsea打分都很低很低！！！

在可视化里面非常明显的差异

让我们看看学员提出来的问题该如何解决吧：.cgi?acc=GSE174449

代码语言：javascript代码运行次数：0运行复制

 
GSM511657 SW480 cells, guo_1
GSM511658 SW480 cells, guo_2
GSM511659 SW480 cells, guo_
GSM511660 SW480 cells, C_1
GSM511661 SW480 cells, C_2
GSM511662 SW480 cells, C_ 

首先呢，下载作者提供的 GSE174449_gene_gz 矩阵文件，然后读取：

代码语言：javascript代码运行次数：0运行复制

data<::fread("GSE174449_gene_gz",
                         = F) 
data$V1=data$gene_name
head(data)
data=data[!duplicated(data$V1),]
mat<-data[,grepl('count',colnames(data)) ]
rownames(mat)=data$V1
mat[1:4,1:4]
keep_feature <- rowSums (mat > 1) > 1 ;table(keep_feature)
symbol_matrix <- mat[keep_feature, ]  
rownames(symbol_matrix)=rownames(mat)[keep_feature]
symbol_matrix[1:4,1:4]

可以看到是很标准的转录组count矩阵：

代码语言：javascript代码运行次数：0运行复制

> symbol_matrix[1:4,1:4]
           C-1_count C-2_count C-_count guo-1_count
AC119427.1         66         50         5          49
MADD             177       1164       1181        1771
MAST2            1288       1075       110        1574
TRBJ2-5             2          1          1           4

所以对它走金标准算法（DESeq2，edgeR，limma-voom）差异分析即可，然后我们把文献里面的热图里面的基因列表去差异分析结果里面看看；

代码语言：javascript代码运行次数：0运行复制

load('deg/DEG_deseq2.Rdata')
head(DEG_deseq2)
DEG =DEG_deseq2
DEG$name = rownames(DEG) 
cg = c("MTOR", "MLST8", "RPS6KB1", "EIF4B", "DLL1", "JAG1", "JAG2", "ADAM17", 
                "OTCH1", "OTCH2", "DTX4", "DTX", "DTXL", "CST", "APH1A", 
                "PSE1", "KAT2B", "CREBBP", "EP00", "MAML", "MAML1", "RBPJ", 
                "SHH", "SMO", "PTCH1")
tmp = DEG[DEG$name %in% cg,] 

可以看到，基本上没什么基因是统计学显著的，因为变化倍数都不怎么样，而且呢，矫正后的p值也只有一半是显著的：

代码语言：javascript代码运行次数：0运行复制

> tmp[,c(1,2,6)]
         baseMean log2FoldChange         padj
SHH     128.8666     0.86794194 1.956649e-06
OTCH1  1496.0912     0.41581512 1.11567e-04
CREBBP   941.122     0.4461568 .887181e-04
PTCH1   227.824     0.64465016 1.449410e-0
JAG1     660.459     0.51754889 1.65606e-0
MAML1   116.0055     0.4762947 4.6285e-0
EP00   1418.2984     0.82901 9.04788e-0
MAML    158.706     0.57017870 2.11556e-02
RPS6KB1  848.9160     0.27926560 2.67459e-02
OTCH2  2202.0024     0.261104 5.208877e-02
KAT2B    109.1602     0.6052208 6.965226e-02
ADAM17  109.1085     0.24147656 8.159929e-02
PSE1   146.2611     0.1857551 1.655e-01
DTX     25.861     0.2871524 1.76214e-01
RBPJ     825.566     0.12187140 4.456628e-01
DTXL    550.7975     0.1428291 4.94770e-01
CST    86.5494     0.0979721 5.28496e-01
MLST8    522.9058     0.11152568 5.40178e-01
EIF4B   7866.1149     0.2051489 5.679472e-01
APH1A    91.165     0.107740 5.786266e-01
MTOR    164.9514     0.0988424 6.021265e-01
JAG2    1175.8502     0.14259294 6.47921e-01
DTX4     218.9629     0.1529427 6.69121e-01
SMO      22.057     0.1145572 6.75498e-01
DLL1     124.0805     0.14822158 7.152226e-01

那么问题来了，这些基因为什么被挑出来了呢，而且为什么在热图里面就看起来在两分组差异很明显呢！如下所示的热图可视化：

代码语言：javascript代码运行次数：0运行复制

library(pheatmap) 
n=t(scale(t(dat[cg,]))) # 'scale'可以对log-ratio数值进行归一化
n[n>2]=2 
n[n< -2]= -2
n[1:4,1:4]
pheatmap(n,show_colnames =F,show_rownames = F)
ac=data.frame(group=group_list) 
rownames(ac)=colnames(n) 
pheatmap(n,show_colnames =F,show_rownames = T,  
         annotation_col=ac) 

秘诀就在于样品之间的zscore，同一个基因仅仅是在不同样品进行比较，如果它在2分组的总计6个样品被zscore了，那么它哪怕是再微小的差异也会被可视化显示出来，因为不同基因完全不需要比较表达量高低。

同一个基因仅仅是在不同样品进行比较

那么，这些基因如此微不足道，为什么会被作者 挑出来呢，秘诀就是gsea分析，它并不需要预先使用统计学指标筛选上下调基因，而是整体之间看通路的变化，只需要一个通路是有变化的，里面的基因再微小的差异也没有关系！

本文参与 腾讯云自媒体同步曝光计划，分享自。原始发表：2024-12-18，如有侵权请联系 cloudcommunity@tencent 删除技巧可视化芯片数据处理data