如果表达矩阵被归一化会发生什么
事情起源于，某个吃了很多汉堡王一起学习的日子，技能树一众学徒一起学习从GEO数据挖掘到limma差异分析等等等。
选的数据集倒是很有、趣。依据原始表达矩阵随便画的热图是这样的?（已经很整齐均一了
t(scale(t(x))) 归一化后，画出来是这样的? （已经美的不真实了
而且呢，logFC也高到惊人
checkFC <- fivenum(DEGdf$logFC)
checkFC
# [1] -13.4823993 -0.5407014 0.1940622 0.6915652 13.0133092
所以其实这个表达矩阵是已经归一化过的。
为了探究事物的本质和宇宙奥义，拿一个优秀的、表现无异常的数据集，分别求由原始矩阵得出的logFC, 和由归一化矩阵的logFC, 看看普通FC和诡异FC之间存不存在规律/线性关系。
用来反推的数据集
GSE12452， 这里根据生信技能树GEO视频全R代码流程：
常规操作
(1)GEO数据下载
gset<- getGEO('GSE12452', destdir=".",
AnnotGPL = F,
getGPL = F)
(2)获得表达矩阵及分组信息
eset <- gset[[1]]
expmat <- exprs(eset)
pd <- pData(eset)
group_list <- ifelse(grepl('normal',pd$title),'normal','npc')
(3) limma差异分析
library(limma)
design <- model.matrix(~0+factor(group_list))
colnames(design) <- levels(factor(group_list))
rownames(design) <- colnames(expmat)
contrast.matrix<-makeContrasts("npc-normal",levels = design)
contrast.matrix
# Contrasts
# Levels npc-normal
# normal -1
# npc 1
jimmy大神把拟合过程包在了一个函数里：
deg = function(expmat,design,contrast.matrix){
##step1
fit <- lmFit(expmat,design)
##step2
fit2 <- contrasts.fit(fit, contrast.matrix)
##这一步很重要，大家可以自行看看效果
fit2 <- eBayes(fit2) ## default no trend !!!
##eBayes() with trend=TRUE
##step3
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput)
#write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
head(nrDEG)
return(nrDEG)
}
差异分析结果
DEGs <- deg(expmat,design,contrast.matrix)
先归一化表达矩阵，再做差异分析
这里选择zcore进行归一化：
NM_expmat <- t(scale(t(expmat)))
NMdesign <- model.matrix(~0+factor(group_list))
colnames(NMdesign) <- levels(factor(group_list))
rownames(NMdesign) <- colnames(NM_expmat)
NMcontrast.matrix<-makeContrasts("npc-normal",levels = NMdesign)
NMcontrast.matrix
# Contrasts
# Levels npc-normal
# normal -1
# npc 1
## 其实设计矩阵、比较矩阵和前面的是一样的
差异分析结果
NM_DEGs <- deg(NM_expmat,NMdesign,NMcontrast.matrix)
比较两个FC
画图
common_FC <- DEGs$logFC
weird_FC <- NM_DEGs$logFC
plot(common_FC , weird_FC)
冷静一下，这不是染色体。
所以这意味着，一些logFC值的正负反过来了，意味着，上调下调反过来了。
检验正负变化
table(common_FC > 0)
# FALSE TRUE
# 24849 19300
table(weird_FC > 0)
# FALSE TRUE
# 24849 19300
神奇的是，两个FC中正负值的数量是一样的，难道没变吗？
再用负负得正检验一下
test <- (common_FC * weird_FC)
table(test > 0)
# FALSE TRUE
# 21382 22767
emmmm这说明有21382个FC值的正负变了….吗?
再确认一下：
乍一看没什么问题，顺序是一致的，所以直接画图，乘除没问题。
但，问题就是，顺序居然是一致的？
再检查一下发现自己太年轻，顺序怎么可能一致！！！！
table(row.names(DEGs)==row.names(NM_DEGs))
# FALSE TRUE
# 43984 165
NEVER BELIEVE WHAT YOU SEE - by 鲁迅
重新排序后：
order1 <- row.names(DEGs)
odNM_DEGs <- NM_DEGs[order1,]
table(row.names(DEGs)==row.names(odNM_DEGs))
# TRUE
# 44149
fntest <- (DEGs$logFC * odNM_DEGs$logFC)
table(fntest > 0)
# TRUE
# 44149
所以并没有正负变换、上调下调变换这么恐怖的事….
手动根据公式算FC
mean_210297_s_at <- mean(expmat[,design[,2]==1]['210297_s_at',])
ckmean_210297_s_at <- mean(expmat[,design[,2]==0]['210297_s_at',])
comFC_210297_s_at <- mean_210297_s_at - ckmean_210297_s_at
comFC_210297_s_at
# [1] -3.566486
DEGs["210297_s_at",1]
# [1] -3.566486
mean_210297_s_at <- mean(NM_expmat[,design[,2]==1]['210297_s_at',])
ckmean_210297_s_at <- mean(NM_expmat[,design[,2]==0]['210297_s_at',])
wrdFC_210297_s_at <- mean_210297_s_at - ckmean_210297_s_at
wrdFC_210297_s_at
# [1] -1.36838
odNM_DEGs["210297_s_at",1]
# [1] -1.36838
都是没问题的。
画图
common_FC <- DEGs$logFC
odered_NMFC <- odNM_DEGs$logFC
plot(common_FC , odered_NMFC)
所以正确的应该是染色质。（雾
plot(abs(correctFC), abs(odered_NMFC))
经归一化的表达矩阵做差异分析，得到的FC大概会落在
里，原始矩阵不同归一化矩阵不同，一顿操作得到的两个FC和它们之间关系肯定也不同，so,大概也并不是很普适(＠_＠;)
plot(abs(correctFC), abs(correctFC)/abs(odered_NMFC))
plot(common_FC, common_FC/odered_NMFC)
目前看起来，以GSE12452为例，表达矩阵被归一化后再做差异分析，FC阈值缩为原来的一半，绝大部分FC0是FC的1.25倍，极少数在3倍以上，FC0和FC0/FC的关系在0.5x以上。
A
E
作者
假装
挖金矿的狮毛
编辑
小洁