两个表达量矩阵去除批次效应之前是否需要归一化

Posted on by

前面我们视频号直播了一个表达量芯片数据处理，详见：这样去除表达量芯片的批次效应可能不妥，这个物信息学数据挖掘的标题是：《Novel ferroptosis gene biomarkers and immune infiltration profiles in diabetic kidney disease via bioinformatics》，直播回放的视频在：

我们之所以要对两个表达量矩阵做去除批次效应的处理，就是因为两个表达量矩阵的取值范围就不一样，而且每个矩阵内部的每个样品或者每个基因的分布范围也不一样，做去除批次效应的处理就是为了抹去两个矩阵的系统性差异。

批次效应（Batch Effect）是指在生物样本的基因表达数据中，由于实验设计、样本处理、数据采集和处理等非生物学因素导致的样本之间的差异。这些差异可能掩盖或模糊了生物学上真实的变异，因此需要通过去除批次效应来揭示数据中真实的生物学信息。以下是去除批次效应处理的具体解释：

  1. 取值范围不同：

    • 不同的表达量矩阵可能由于实验条件、测量技术或数据标准化流程的差异，导致每个矩阵的基因表达量取值范围不同。例如，一个矩阵的表达量可能在1-10,000，而另一个矩阵的表达量可能在0.1-100。

  2. 矩阵内部样本或基因分布差异：

    • 即使在同一个矩阵内部，不同样本或基因也可能表现出不同的表达量分布特征，如均值、方差、偏度等统计特性。

  3. 系统性差异：

    • 批次效应导致的系统性差异是指由于批次因素引起的一致性偏差，这些偏差可能在不同批次的样本之间导致可重复性问题，影响后续分析的准确性。

  4. 去除批次效应的目的：

    • 抹去系统性差异：通过各种统计和计算方法，如主成分分析（PCA）、多变量回归模型、批次校正算法等，来调整和消除批次效应的影响。
    • 增强可比性：去除批次效应后，不同批次、不同平台甚至不同实验室的数据可以进行比较和综合分析，提高了数据的可比性。
    • 揭示真实变异：去除批次效应有助于更准确地识别生物学上真实的变异，如疾病状态、药物响应等。

  5. 常用方法：

    • 数据标准化：如使用Z分数（Z-score normalization）或量化归一化（quantile normalization）等方法，使不同数据集的表达量数据具有可比性。
    • 批次校正算法：如Combat、MNN（Minimum Covariance Determinant）等，这些算法可以识别并调整批次效应，减少其对数据分析的影响。
    • 多变量分析：利用PCA、因子分析等方法识别批次效应，并在后续分析中考虑或排除这些效应。

  6. 后续分析：

    • 在去除批次效应后，可以进行差异表达分析、聚类分析、路径分析等，以探索样本之间的生物学关系和基因功能。

总之，去除批次效应是基因表达数据分析中的重要步骤，它有助于提高数据质量，确保研究结果的可靠性和生物学意义。

那么，问题就来了，两个表达量矩阵去除批次效应之前是否需要归一化呢？

处理GSE47185表达量矩阵

直接使用作者上传的表达量矩阵即可，如下所示的代码，因为这个GSE47185表达量矩阵样品数量非常多，分组很复杂，但是就选择了第一个数据集的Diabetic的14个样品，全部的代码如下所示：

library(AnnoProbe)
library(GEOquery) 
library(ggplot2)
library(ggstatsplot)
library(patchwork)
library(reshape2)
library(stringr)
getOption('timeout')
options(timeout=10000)

gse_number <- 'GSE47185' 
list.files() 
if(F){gset <- geoChina(gse_number)}
if(T){ gset = getGEO("GSE47185", destdir = '.', getGPL = F) }
a=gset[[1]] 
dat=exprs(a) #a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数
dim(dat)#看一下dat这个矩阵的维度
dat[1:4,1:4] #查看dat这个矩阵的1至4行和1至4列，逗号前为行，逗号后为列
pd = pData(a)
head(pd)
kp = grepl('Diabetic ', pd$title )
table(kp)
pd=pd[kp,]
dat=dat[,kp]

library(org.Hs.eg.db)
library(clusterProfiler)
ids <- bitr( rownames(dat) ,fromType = "ENTREZID",
 toType = c("SYMBOL"),
 OrgDb = org.Hs.eg.db)
head(ids)
colnames(ids) =c('ID','symbol')

上面的表达量矩阵并不是基因的symbol，所以我们做了一个简单的id的转换哦。很简单的转换代码，如下所示：


if(T){
 #探针基因ID对应以及去冗余
 ids=ids[ids$symbol != '',]
 dat=dat[rownames(dat) %in% ids$ID,]
 ids=ids[match(rownames(dat),ids$ID),]
 head(ids) 
 colnames(ids)=c('probe_id','symbol') 
 ids$probe_id=as.character(ids$probe_id)
 rownames(dat)=ids$probe_id
 dat[1:4,1:4] 

 ids=ids[ids$probe_id %in% rownames(dat),]
 dat[1:4,1:4] 
 dat=dat[ids$probe_id,] 
 probe_exp = dat
 probe_info = ids

 ids$median=apply(dat,1,median) #ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行
 ids=ids[order(ids$symbol,ids$median,decreasing = T),]#对ids$symbol按照ids$median中位数从大到小排列的顺序排序，将对应的行赋值为一个新的ids
 ids=ids[!duplicated(ids$symbol),]#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果s
 dat=dat[ids$probe_id,] #新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat
 rownames(dat)=ids$symbol#把ids的symbol这一列中的每一行给dat作为dat的行名
 dat[1:4,1:4] #保留每个基因ID第一次出现的信息

 fivenum(dat['ACTB',])
 fivenum(dat['GAPDH',])
}

save(gse_number,dat,#group_list,
 # probe_exp,probe_info,pd,
 file = 'step1-output.Rdata')

与第一个表达量矩阵合并（基于zscore表达量矩阵）

只需要读取两个表达量矩阵，然后使用sva包的ComBat函数即可

rm(list = ls()) 
options(stringsAsFactors = F) 
load('../03-GSE47185/GSE47185/step1-output.Rdata')
GSE47185_dat = dat
GSE47185_g = as.factor(rep('Diabetic',ncol(GSE47185_dat)))

load('../01-GSE30122/GSE30122/step1-output.Rdata')
GSE30122_dat = dat
GSE30122_g = group_list
table(GSE30122_g)

ids=intersect(rownames(GSE30122_dat),
 rownames(GSE47185_dat))
merge_dat = cbind(GSE30122_dat[ids,],
 GSE47185_dat[ids,] )
meta=data.frame(
 gse_number = c(rep('GSE30122',length(GSE30122_g)),
 rep('GSE47185',length(GSE47185_g))), 
 group = c(GSE30122_g,GSE47185_g)
)

library(sva)
table(meta)
table(meta$gse_number)
combat_edata1 = ComBat(dat=as.matrix(merge_dat), 
 batch=meta$gse_number, mod=NULL, 
 par.prior=T, prior.plots=F)
combat_edata1[1:4,1:4]
boxplot(combat_edata1)
dat=combat_edata1
group_list=meta$group

boxplot(combat_edata1,col=as.numeric(as.factor(meta$gse_number)))

因为GSE30122数据集本身就是使用了作者提供的zscore后的，所以批次效应抹除后也是总体上的矩阵偏向于0附近哦，如下所示：

image-20240526194543204

值得注意是的，前面的GSE47185数据集的Diabetic的14个样品居然也是有批次效应的， 但是后面主要是做 control 和 Diabetic的差异分析，所以GSE47185数据集内部的异质性这个时候无伤大雅哈：

> table(meta)
 group
gse_number control Diabetic
 GSE30122 50 19
 GSE47185 0 14

与第二个表达量矩阵合并（基于基于cel文件）

同样的，读取两个表达量矩阵后有使用sva包的ComBat函数即可

image-20240526194402620

两个不同策略后的差异分析结果的对比

同样的对比方式：

 zscore_deg
cel_deg down stable up
 down 104 73 0
 stable 46 10358 49
 up 0 53 244

而且也是可以看到， 冲突的基因列表和一致性的基因列表，都是有生物学功能的：

image-20240526195013864

原文里面的做法肯定是不可取的，一个矩阵是zscore的另外一个不是，主要的两个矩阵去去除批次效应总是有点奇怪，详见：这样去除表达量芯片的批次效应可能不妥，这个物信息学数据挖掘的标题是：《Novel ferroptosis gene biomarkers and immune infiltration profiles in diabetic kidney disease via bioinformatics》

​

Comments are closed.