0.背景
神奇的曾老板好久沒出現(xiàn)在我的推文里了检柬，今日份出沒：多次差異分析難道就需要多個(gè)火山圖嗎
就是呢婿斥，當(dāng)分組非常多的時(shí)候要展示汛兜，可以以這個(gè)barplot的形式去展示巴粪，省版面費(fèi)，也非常直觀。
然后他就找了一個(gè)數(shù)據(jù)GSE116439肛根，作為學(xué)徒作業(yè)了辫塌。既然是學(xué)徒作業(yè)，為啥是我這個(gè)講師來做呢派哲。臼氨。。說來話長芭届，那是一個(gè)月黑風(fēng)高的晚上储矩，我看到這張圖發(fā)表了評(píng)論。褂乍。持隧。
那。逃片。屡拨。所以我就做了。
1.數(shù)據(jù)下載和分組信息獲取
這個(gè)表達(dá)矩陣有50多M褥实，用AnnoProbe::geoChina()下載更快更絲滑~
#devtools::install_github("jmzeng1314/AnnoProbe")
library(tidyverse)
library(GEOquery)
library(AnnoProbe)
g = geoChina("GSE116439",destdir = ".")[[1]]
exp = exprs(g)
gr = pData(g)[1]
gr = separate(gr,title,sep = "_",letters[1:4])
head(gr)
## a b c d
## GSM3232610 A498 cisplatin 0nM 24h
## GSM3232611 A498 cisplatin 0nM 2h
## GSM3232612 A498 cisplatin 0nM 6h
## GSM3232613 A498 cisplatin 15000nM 24h
## GSM3232614 A498 cisplatin 15000nM 2h
## GSM3232615 A498 cisplatin 15000nM 6h
分組已經(jīng)很清楚的寫在了title列呀狼，這是把他們拆分開成4列的結(jié)果∷鹄耄可以看到我們要的是濃度和時(shí)間的信息哥艇，在3-4列。這里我把3-4列放在一起咯草冈。
Group = paste(gr$c,gr$d,sep = "_")
Group = factor(Group,levels = c("0nM_2h", "0nM_6h", "0nM_24h", "15000nM_2h", "15000nM_6h",
"15000nM_24h", "3000nM_2h", "3000nM_6h", "3000nM_24h"))
table(Group)
## Group
## 0nM_2h 0nM_6h 0nM_24h 15000nM_2h 15000nM_6h 15000nM_24h
## 55 56 56 55 56 55
## 3000nM_2h 3000nM_6h 3000nM_24h
## 54 56 56
上面是分組信息她奥，下面是design矩陣，第一列是截距怎棱，無需理會(huì)哩俭。第二列往后，對(duì)應(yīng)著1的行拳恋，說明表達(dá)矩陣?yán)飳?duì)應(yīng)的列（樣本）屬于這個(gè)分組凡资。比如第三個(gè)樣本就屬于Group0nM_6h。
design = model.matrix(~Group)
head(design)
## (Intercept) Group0nM_6h Group0nM_24h Group15000nM_2h Group15000nM_6h
## 1 1 0 1 0 0
## 2 1 0 0 0 0
## 3 1 1 0 0 0
## 4 1 0 0 0 0
## 5 1 0 0 1 0
## 6 1 0 0 0 1
## Group15000nM_24h Group3000nM_2h Group3000nM_6h Group3000nM_24h
## 1 0 0 0 0
## 2 0 0 0 0
## 3 0 0 0 0
## 4 1 0 0 0
## 5 0 0 0 0
## 6 0 0 0 0
這是design矩陣谬运，是所有其他分組和0nM_2h對(duì)比隙赁。后文畫圖就按照這種分組方式來。
還可以有別的比較方式梆暖，比如同時(shí)考慮這兩個(gè)因素進(jìn)行差異分析伞访，就是說在做濃度比較時(shí)，是按照時(shí)間配對(duì)去比較的轰驳，同樣厚掷，在做時(shí)間比較的時(shí)候弟灼，也是按照濃度配對(duì)去比較。兩種方式任選咯冒黑。
nm = factor(pd$c,levels = c("0nM", "15000nM", "3000nM"))
tm = factor(pd$d,levels = c("2h", "6h", "24h"))
design2 = model.matrix(~nm+tm)
head(design2)
# (Intercept) nm15000nM nm3000nM tm6h tm24h
# 1 1 0 0 0 1
# 2 1 0 0 0 0
# 3 1 0 0 1 0
# 4 1 1 0 0 1
# 5 1 1 0 0 0
# 6 1 1 0 1 0
直接pia~復(fù)制過來的代碼完成差異分析和畫圖(顏色我掰回來了田绑，強(qiáng)迫癥不可能允許紅色表示下調(diào))
library(limma)
glm = lmFit(exp , design = design )
glm = eBayes(glm)
testResults <- decideTests(glm, method="hierarchical",adjust.method="BH", p.value=0.05)[,-1]
significantGenes <- sapply(1:ncol(testResults), function(j){
c <- glm$coefficients[testResults[,j]!=0,j+1]
table(cut(c, breaks=c(-5,seq(-1.5,1.5,l=7),5)))
})
colnames(significantGenes) <- colnames(testResults)
# Barplot
library(RColorBrewer)
library(Hmisc)
# author mg14 , https://rdrr.io/github/mg14/mg14/
rotatedLabel <- function(x0 = seq_along(labels), y0 = rep(par("usr")[3], length(labels)), labels, pos = 1, cex=1, srt=45, ...) {
w <- strwidth(labels, units="user", cex=cex)
h <- strheight(labels, units="user",cex=cex)
u <- par('usr')
p <- par('plt')
f <- par("fin")
xpd <- par("xpd")
par(xpd=NA)
text(x=x0 + ifelse(pos==1, -1,1) * w/2*cos(srt/360*2*base::pi), y = y0 + ifelse(pos==1, -1,1) * w/2 *sin(srt/360*2*base::pi) * (u[4]-u[3])/(u[2]-u[1]) / (p[4]-p[3]) * (p[2]-p[1])* f[1]/f[2] , labels, las=2, cex=cex, pos=pos, adj=1, srt=srt,...)
par(xpd=xpd)
}
par(bty="n", mgp = c(2.5,.33,0), mar=c(3,3.3,2,0)+.1, las=2, tcl=-.25)
b <- barplot(significantGenes, las=2, ylab = "Differentially expressed genes",
col=rev(brewer.pal(8,"RdYlBu")),
legend.text=FALSE , border=0, xaxt="n")#, col = set1[simple.annot[names(n)]], border=NA)
print(b)
## [1] 0.7 1.9 3.1 4.3 5.5 6.7 7.9 9.1
rotatedLabel(x0=b, y0=rep(10, ncol(significantGenes)),
labels=colnames(significantGenes), cex=.7, srt=45,
font=ifelse(grepl("[[:lower:]]", colnames(design))[-1], 1,3) )
clip(0,30,0,1000)
#text(b+0.2, colSums(n)+50, colSums(n), pos=3, cex=.7, srt=90)
x0 <- 21.5
image(x=x0+c(0,0.8), y=par("usr")[4]+seq(-100,100,l=9), z=matrix(1:8, ncol=8), col=brewer.pal(8,"RdYlBu"), add=TRUE)
text(x=x0+1.5, y=par("usr")[4]+seq(-50,50,l=3), format(seq(-1,1,l=3),2), cex=0.66)
lines(x=rep(x0+.8,2), y=par("usr")[4]+c(-75,75))
segments(x0+.8,par("usr")[4]+seq(-75,75,l=7),x0+.9,par("usr")[4]+seq(-75,75,l=7))
text(x0+.8, par("usr")[4]+125, "log2 FC", cex=.66)
rotatedLabel(b-0.1, colSums(significantGenes), colSums(significantGenes), pos=3, cex=, srt=45)
ggplot2 畫同款圖
哈，基礎(chǔ)包代碼是現(xiàn)成的抡爹，沒有玩盡興掩驱，所以用ggplot2代碼再玩一玩。
library(tidyverse)
library(ggplot2)
df = significantGenes %>%
as.data.frame() %>%
rownames_to_column(var = "logFC") %>%
pivot_longer(cols = starts_with("Group"),
names_to = "Group",
values_to = "count")
my_color = rev(brewer.pal(8,"RdYlBu"))
names(my_color) = unique(df$logFC)
df$Group = str_remove(df$Group,"Group")
df$logFC = factor(df$logFC,levels = rev(unique(df$logFC)))
df2 = group_by(df,Group) %>%
summarise(count = sum(count))
ggplot(df)+
geom_bar(aes(x = Group,y = count,fill = logFC),stat="identity")+
geom_text(data = df2,aes(x = Group,y = count,label = count),vjust = -1,angle = 45)+
theme_bw()+
scale_fill_manual(values = my_color)+
theme(axis.text.x = element_text(angle=50,vjust = 0.5))+
guides(fill = guide_legend(reverse=T))
ggplot2這些代碼解決的問題：
1.矩陣畫圖冬竟，得改成數(shù)據(jù)框欧穴，并且寬變長。
2.堆疊式直方圖的疊放次序诱咏，是用因子水平規(guī)定的苔可，不規(guī)定就會(huì)自動(dòng)缴挖。
3.直方圖頭頂加數(shù)字袋狞，用summarise配group_by實(shí)現(xiàn)計(jì)算，冷知識(shí)：ggplot2不同圖層可以使用不同的數(shù)據(jù)畫映屋。
4.最后一句代碼實(shí)現(xiàn)顏色圖例順序逆轉(zhuǎn)
寫點(diǎn)閑話
這個(gè)周是我們線上直播課的間隙苟鸯，又是豆豆辦港澳通行證等待的一個(gè)星期，等同于兩人一起放假棚点。我們已經(jīng)來了珠海兩年早处，長隆海洋王國就在幾公里外，愣是到現(xiàn)在才想起來去瘫析，可能就是傳說中的家門口的景點(diǎn)永遠(yuǎn)不想去吧砌梆。天實(shí)在是太熱了，戶外活動(dòng)基本告別贬循，豆豆好不容易放假又不想浪費(fèi)咸包，所以直接就去了。長隆的單日門票400杖虾，年卡800烂瘫，所以我們選了年卡哈哈。分享兩張有意思的照片
祝大家笑口常開奇适。另外下周一（8月2號(hào)）開始又是一輪新的數(shù)據(jù)挖掘和生信入門線上課坟比，如果你的暑假還沒過完，不妨來跟我學(xué)習(xí)咯嚷往。
