电竞比分网-中国电竞赛事及体育赛事平台

搜索
分享

【原】如果表達(dá)矩陣被歸一化會發(fā)生什么

 健明 2021-07-14
事情起源于，某個(gè)吃了很多漢堡王一起學(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í)這個(gè)表達(dá)矩陣是已經(jīng)歸一化過的。

為了探究事物的本質(zhì)和宇宙奧義，拿一個(gè)優(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大神把擬合過程包在了一個(gè)函數(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è)計(jì)矩陣、比較矩陣和前面的是一樣的

差異分析結(jié)果

NM_DEGs <- deg(NM_expmat,NMdesign,NMcontrast.matrix)

比較兩個(gè)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 

神奇的是，兩個(gè)FC中正負(fù)值的數(shù)量是一樣的，難道沒變嗎？

再用負(fù)負(fù)得正檢驗(yàn)一下

test <- (common_FC * weird_FC)
table(test > 0)
# FALSE  TRUE 
# 21382 22767 

emmmm這說明有21382個(gè)FC值的正負(fù)變了….嗎?

再確認(rèn)一下：

乍一看沒什么問題，順序是一致的，所以直接畫圖，乘除沒問題。

但，問題就是，順序居然是一致的？

再檢查一下發(fā)現(xiàn)自己太年輕，順序怎么可能一致?。。?！

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)變換這么恐怖的事….

手動(dòng)根據(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大概會落在

里，原始矩陣不同歸一化矩陣不同，一頓操作得到的兩個(gè)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以上。

A

E

作者

假裝

挖金礦的獅毛

編輯

小潔

忘了怎么分身

■   ■   ■ 

    轉(zhuǎn)藏 分享 獻(xiàn)花（0）

    0條評論

    發(fā)表

    請遵守用戶 評論公約

    類似文章 更多