这个差异表达分析一看就是错误的？

生信技能树

发布于 2025-05-18 01:57:07

13300

代码可运行

文章被收录于专栏：生信技能树生信技能树

运行总次数：0

代码可运行

今天是数据挖掘分析主题，曾老板几天前发来了一个数据集，让我看一下，他说：这个数据的文献中的差异分析一看就是错误的！来看看老板的火眼金睛怎么练出来的~

数据介绍

这个数据对应的文献于2019年8月6号发表在Cell Rep.杂志上，文献标题为《Modeling Human Cancer-induced Cachexia》。对应的10个样本的转录组测序数据已经上传到了GEO数据库，链接：

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE133523

具体的样本情况如下：5 个实验组 vs 5个对照组

GSM3911258 control patient muscle 1
GSM3911259 control patient muscle 2
GSM3911260 control patient muscle 3
GSM3911261 control patient muscle 4
GSM3911262 control patient muscle 5
GSM3911263 cachectic patient muscle 1
GSM3911264 cachectic patient muscle 2
GSM3911265 cachectic patient muscle 3
GSM3911266 cachectic patient muscle 4
GSM3911267 cachectic patient muscle 5

文献中做完差异分析后，使用热图展示了差异结果，就是下面这个热图，曾老板的火眼金睛一下子就看出来它有问题！为什么他看出来了不同寻常？从哪里看出来的？

Figure 1. Comparing Xenograft Models with Cachectic Pancreatic Cancer Patients. (B) RNA-seq of n = 5 male non-cancer patients and 5 male cachectic PDA patients.

差异结果：Interestingly, when similar gene expression profiling was performed with cachectic PDA patients, only 141 DEGs were identified compared with weight-stable patients without cancer (Figure 1B; Table S1)

使用 limma 算法得到了141个显著差异表达基因！

胰腺癌恶病质（Cachectic Pancreatic Cancer）

看一眼数据对应的疾病背景：

胰腺癌恶病质是一种复杂的多因素综合征，其特征是在厌食的背景下，脂肪组织和骨骼肌的持续消耗，导致功能逐渐受损。根据国际共识，癌症恶病质被定义为一种多因素综合征，其特征是骨骼肌质量的持续丢失，这种丢失不能通过传统的营养支持完全逆转，并导致功能逐渐受损。诊断标准包括在过去6个月内体重下降超过5%，或体重指数（BMI）小于20 kg/m²的个体体重下降超过2%，或存在肌肉减少症（sarcopenia）且体重下降超过 2% 。

胰腺癌恶病质的临床表现包括：

体重下降：超过5%的体重下降是常见的；
肌肉消耗：骨骼肌质量的持续丢失；
厌食：食欲减退，食物摄入减少；
功能受损：体力下降，活动能力减弱；
代谢紊乱：包括胰岛素抵抗、脂肪组织棕色化、氧化应激增加等。

数据下载与预处理

先把数据下载下来看看数据质量如何：

rm(list = ls())#清空当前的工作环境
options(scipen = 20)#不以科学计数法显示
library(data.table)
library(tinyarray)
library(stringr)
library(GEOquery)
library(AnnoProbe)
library(org.Hs.eg.db)
# 创建目录
getwd()
gse <- "GSE133523"
dir.create(gse)

# 数据count下载以及基因symbol转换
if(T){
# load counts table from GEO
#urld <- "https://www.ncbi.nlm.nih.gov/geo/download/?type=rnaseq_counts&acc=GSE224401&format=file&file=GSE224401_raw_counts_GRCh38.p13_NCBI.tsv.gz"
  urld <- "https://www.ncbi.nlm.nih.gov/geo/download/?type=rnaseq_counts&"
  gse <- "GSE133523"
  path <- paste0(urld, "acc=", gse, "&format=file&file=",gse,"_raw_counts_GRCh38.p13_NCBI.tsv.gz");
  file <- paste0(gse,"/",gse,"_raw_counts_GRCh38.p13_NCBI.tsv.gz")
  path
  file
  download.file(url = path, destfile = file)

  tbl <- as.matrix(data.table::fread(file, header=T, colClasses="integer"), rownames=1)
  dim(tbl)
  ensembl_matrix=as.data.frame(tbl) 

  keytypes(org.Hs.eg.db)
  e2s<-AnnotationDbi::select(org.Hs.eg.db,
                             keys= rownames(ensembl_matrix),
                             columns="SYMBOL",
                             keytype = "ENTREZID")
  head(e2s)
  ids = na.omit(e2s)
  ids=ids[!duplicated(ids$SYMBOL),]
  ids=ids[!duplicated(ids$ENTREZID),]
  head(ids)
  symbol_matrix= ensembl_matrix[match(ids$ENTREZID,rownames(ensembl_matrix)),]
  rownames(symbol_matrix) = ids$SYMBOL
  symbol_matrix[1:4,1:4]
}


# 获取样本分组信息
if(T){
  gse
  gset = getGEO(gse, destdir = './GSE133523/', getGPL = F) 
  pd = pData(gset[[1]]) 
  pd[1:4,1:5]
  com <- intersect(rownames(pd), colnames(symbol_matrix))
  symbol_matrix <- symbol_matrix[, com]

  pd=pd[com,]
  pd$title
  group_list=ifelse(grepl('control',pd$title ), 'control','case' )
}

colnames(symbol_matrix) 
group_list
table(group_list)
group_list = factor(group_list,levels = c('control','case' )) 
dat=log10(edgeR::cpm(symbol_matrix)+1)
save(symbol_matrix,dat,group_list,file = 'GSE133523/step1-output.Rdata')

这样数据就处理好啦：

先看一下数据质量：

rm(list = ls())  ## 魔幻操作，一键清空~
library("FactoMineR")#画主成分分析图需要加载这两个包
library("factoextra")  
dir.create("1-QC")
list.files("1-QC/")

load(file = 'GSE133523/step1-output.Rdata')
table(group_list)
# 每次都要检测数据
dat[1:4,1:4] 
dim(dat)
#apply按行（'1'是按行取，'2'是按列取）取每一行的方差，从小到大排序，取最大的1000个
cg=names(tail(sort(apply(dat,1,sd)),1000))
## 下面是画PCA的必须操作，需要看说明书。
dat=t(dat[cg,])#画PCA图时要求是行名时样本名，列名时探针名，因此此时需要转换
dat=as.data.frame(dat)#将matrix转换为data.frame 
dat.pca <- PCA(dat , graph = FALSE)#现在dat最后一列是group_list，需要重新赋值给一个dat.pca,这个矩阵是不含有分组信息的
fviz_pca_ind(dat.pca, geom.ind = "point", # show points only (nbut not "text")
             col.ind =  group_list, addEllipses = T, legend.title = "Groups")
ggsave('1-QC/all_samples_PCA.pdf',height = 5,width = 5)

PCA分析结果：

差异分析

然后是简单的差异分析，差异基因热图如下：

结合上面的PCA图，可以看到 实验组的的组内样本异质性比较大，差异基因热图显示出不同的基因表达模式完全不均一！

也就是说根本不可能达到文章中那样在某一组中的所有的表达都高或者都低那么均一的情况！

那曾老板是咋发现的呢，没有做数据分析就发现了吗？

此外，文章还提供了一个差异表达基因结果列表：这两列值非常诡异！

下面这个是我的结果：

我问了一下kimi：校正后的pvalue有很多相同的值是什么原因呢？

在校正后的p值（p-value）中出现很多相同的值，可能有以下几个原因：

1. 数据的离散性

如果原始数据是离散的（例如计数数据），那么在进行统计检验时，计算出的p值可能会集中在某些特定的值上。例如，在进行卡方检验时，如果数据的类别较少，计算出的p值可能会重复。
例子：假设你有一个简单的2×2列联表，数据的分布使得卡方统计量只能取几个有限的值，那么对应的p值也会是有限的几个值。

2. 校正方法的特性

一些校正方法（如Bonferroni校正）会将原始p值乘以一个固定的倍数（通常是检验的次数）。如果原始p值中有相同的值，校正后的p值也会相同。
例子：假设你有多个假设检验的原始p值分别为0.01、0.01、0.02，进行Bonferroni校正时，校正后的p值分别为0.01×n、0.01×n、0.02×n（n为检验次数），其中0.01×n的值会重复。

3. p值的精度限制

在计算p值时，计算机的浮点数精度是有限的。当p值非常小（如小于10⁻¹⁵）时，计算机可能无法区分非常接近的值，从而导致校正后的p值看起来相同。
例子：假设两个原始p值分别为10⁻¹⁶和10⁻¹⁷，在计算机中可能都被近似为0，校正后也会相同。

4. 数据的重复性或相似性

如果你的数据中有大量重复的观测值，或者不同组之间的差异非常小，那么统计检验的结果可能会产生相同的p值。
例子：在进行t检验时，如果两组数据的均值非常接近，且方差较小，计算出的p值可能会很高且相同。

5. 校正方法的分组特性

一些校正方法（如Benjamini-Hochberg校正）会根据p值的排序进行校正。如果原始p值中有许多相同的值，校正后的p值也可能会相同。
例子：假设原始p值为0.01、0.01、0.02、0.02，按照Benjamini-Hochberg校正方法，校正后的p值可能会因为排序和校正公式而产生相同的值。