事情起源于，某个~~吃了很多汉堡王~~一起学习的日子，技能树一众学徒一起学习从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数据下载

gset<- getGEO('GSE12452', destdir=".",
              AnnotGPL = F,   
              getGPL = F)

获得表达矩阵及分组信息

eset <- gset[[1]]
expmat <- exprs(eset)
pd <- pData(eset)
group_list <- ifelse(grepl('normal',pd$title),'normal','npc')

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)

先归一化表达矩阵，再做差异分析

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

$log2FoldChange = log2(Sample1/Sample2)$

$= log2(Sample1) - log2(Sample2)$

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大概会落在 $\int _0^4\left(2+2x-\frac{1}{2}x\right)dx$ 里，原始矩阵不同归一化矩阵不同，一顿操作得到的两个FC和它们之间关系肯定也不同，so,大概也并不是很普适(＠_＠;)

plot(abs(correctFC), abs(correctFC)/abs(odered_NMFC))

plot(common_FC, common_FC/odered_NMFC)

目前看起来，以GSE12452为例，表达矩阵被归一化后再做差异分析，FC阈值缩为原来的一半，绝大部分FC₀是FC的1.25倍，极少数在3倍以上，FC₀和FC₀/FC的关系在0.5x以上。

最后，向大家隆重推荐生信技能树的一系列干货！

生信技能树全球公益巡讲：https://mp.weixin.qq.com/s/E9ykuIbc-2Ja9HOY0bn_6g
B站公益74小时生信工程师教学视频合辑：https://mp.weixin.qq.com/s/IyFK7l_WBAiUgqQi8O7Hxw
招学徒：https://mp.weixin.qq.com/s/KgbilzXnFjbKKunuw7NVfw

如果表达矩阵被归一化会发生什么