两个表达量矩阵去除批次效应之前是否需要归一化
前面我们视频号直播了一个表达量芯片数据处理，详见：这样去除表达量芯片的批次效应可能不妥，这个物信息学数据挖掘的标题是：《Novel ferroptosis gene biomarkers and immune infiltration profiles in diabetic kidney disease via bioinformatics》，直播回放的视频在：
我们之所以要对两个表达量矩阵做去除批次效应的处理，就是因为两个表达量矩阵的取值范围就不一样，而且每个矩阵内部的每个样品或者每个基因的分布范围也不一样，做去除批次效应的处理就是为了抹去两个矩阵的系统性差异。
批次效应（Batch Effect）是指在生物样本的基因表达数据中，由于实验设计、样本处理、数据采集和处理等非生物学因素导致的样本之间的差异。这些差异可能掩盖或模糊了生物学上真实的变异，因此需要通过去除批次效应来揭示数据中真实的生物学信息。以下是去除批次效应处理的具体解释：
- 取值范围不同：
- 不同的表达量矩阵可能由于实验条件、测量技术或数据标准化流程的差异，导致每个矩阵的基因表达量取值范围不同。例如，一个矩阵的表达量可能在1-10,000，而另一个矩阵的表达量可能在0.1-100。
- 矩阵内部样本或基因分布差异：
- 即使在同一个矩阵内部，不同样本或基因也可能表现出不同的表达量分布特征，如均值、方差、偏度等统计特性。
- 系统性差异：
- 批次效应导致的系统性差异是指由于批次因素引起的一致性偏差，这些偏差可能在不同批次的样本之间导致可重复性问题，影响后续分析的准确性。
- 去除批次效应的目的：
- 抹去系统性差异：通过各种统计和计算方法，如主成分分析（PCA）、多变量回归模型、批次校正算法等，来调整和消除批次效应的影响。
- 增强可比性：去除批次效应后，不同批次、不同平台甚至不同实验室的数据可以进行比较和综合分析，提高了数据的可比性。
- 揭示真实变异：去除批次效应有助于更准确地识别生物学上真实的变异，如疾病状态、药物响应等。
- 常用方法：
- 数据标准化：如使用Z分数（Z-score normalization）或量化归一化（quantile normalization）等方法，使不同数据集的表达量数据具有可比性。
- 批次校正算法：如Combat、MNN（Minimum Covariance Determinant）等，这些算法可以识别并调整批次效应，减少其对数据分析的影响。
- 多变量分析：利用PCA、因子分析等方法识别批次效应，并在后续分析中考虑或排除这些效应。
- 后续分析：
- 在去除批次效应后，可以进行差异表达分析、聚类分析、路径分析等，以探索样本之间的生物学关系和基因功能。
总之，去除批次效应是基因表达数据分析中的重要步骤，它有助于提高数据质量，确保研究结果的可靠性和生物学意义。
那么，问题就来了，两个表达量矩阵去除批次效应之前是否需要归一化呢？
处理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附近哦，如下所示：
值得注意是的，前面的GSE47185数据集的Diabetic的14个样品居然也是有批次效应的， 但是后面主要是做 control 和 Diabetic的差异分析，所以GSE47185数据集内部的异质性这个时候无伤大雅哈：
> table(meta)
group
gse_number control Diabetic
GSE30122 50 19
GSE47185 0 14
与第二个表达量矩阵合并（基于基于cel文件）
同样的，读取两个表达量矩阵后有使用sva包的ComBat函数即可
两个不同策略后的差异分析结果的对比
同样的对比方式：
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》