事情起源于，某個吃了很多漢堡王一起學(xué)習(xí)的日子，技能樹一眾學(xué)徒一起學(xué)習(xí)從GEO數(shù)據(jù)挖掘到limma差異分析等等等。選的數(shù)據(jù)集倒是很有、趣。依據(jù)原始表達(dá)矩陣隨便畫的熱圖是這樣的??（已經(jīng)很整齊均一了  t(scale(t(x))) 歸一化后，畫出來是這樣的?? （已經(jīng)美的不真實(shí)了
 而且呢，logFC也高到驚人 checkFC <- fivenum(DEGdf$logFC)
checkFC
# [1] -13.4823993 -0.5407014 0.1940622 0.6915652 13.0133092
所以其實(shí)這個表達(dá)矩陣是已經(jīng)歸一化過的。 為了探究事物的本質(zhì)和宇宙奧義，拿一個優(yōu)秀的、表現(xiàn)無異常的數(shù)據(jù)集，分別求由原始矩陣得出的logFC, 和由歸一化矩陣的logFC, 看看普通FC和詭異FC之間存不存在規(guī)律/線性關(guān)系。
用來反推的數(shù)據(jù)集GSE12452， 這里根據(jù)生信技能樹GEO視頻全R代碼流程： 
常規(guī)操作(1)GEO數(shù)據(jù)下載gset<- getGEO('GSE12452', destdir=".",
AnnotGPL = F,
getGPL = F)
(2)獲得表達(dá)矩陣及分組信息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大神把擬合過程包在了一個函數(shù)里： 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)
}
差異分析結(jié)果DEGs <- deg(expmat,design,contrast.matrix)
先歸一化表達(dá)矩陣，再做差異分析這里選擇zcore進(jìn)行歸一化：
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
## 其實(shí)設(shè)計矩陣、比較矩陣和前面的是一樣的
差異分析結(jié)果NM_DEGs <- deg(NM_expmat,NMdesign,NMcontrast.matrix)
比較兩個FC畫圖common_FC <- DEGs$logFC
weird_FC <- NM_DEGs$logFC
plot(common_FC , weird_FC)
  冷靜一下，這不是染色體。 所以這意味著，一些logFC值的正負(fù)反過來了，意味著，上調(diào)下調(diào)反過來了。  檢驗(yàn)正負(fù)變化table(common_FC > 0)
# FALSE TRUE
# 24849 19300
table(weird_FC > 0)
# FALSE TRUE
# 24849 19300
神奇的是，兩個FC中正負(fù)值的數(shù)量是一樣的，難道沒變嗎？ 再用負(fù)負(fù)得正檢驗(yàn)一下 test <- (common_FC * weird_FC)
table(test > 0)
# FALSE TRUE
# 21382 22767
emmmm這說明有21382個FC值的正負(fù)變了….嗎? 再確認(rèn)一下：  乍一看沒什么問題，順序是一致的，所以直接畫圖，乘除沒問題。 但，問題就是，順序居然是一致的？ 再檢查一下發(fā)現(xiàn)自己太年輕，順序怎么可能一致?。。。?/p>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
所以并沒有正負(fù)變換、上調(diào)下調(diào)變換這么恐怖的事…. 手動根據(jù)公式算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)
 所以正確的應(yīng)該是染色質(zhì)。（霧 plot(abs(correctFC), abs(odered_NMFC))
 經(jīng)歸一化的表達(dá)矩陣做差異分析，得到的FC大概會落在
里，原始矩陣不同歸一化矩陣不同，一頓操作得到的兩個FC和它們之間關(guān)系肯定也不同，so,大概也并不是很普適(＠_＠;) plot(abs(correctFC), abs(correctFC)/abs(odered_NMFC))
plot(common_FC, common_FC/odered_NMFC)
目前看起來，以GSE12452為例，表達(dá)矩陣被歸一化后再做差異分析，FC閾值縮為原來的一半，絕大部分FC0是FC的1.25倍，極少數(shù)在3倍以上，F(xiàn)C0和FC0/FC的關(guān)系在0.5x以上。
■ ■ ■
|
