0.背景
神奇的曾老板好久没出现在我的推文里了，今日份出没：多次差异分析难道就需要多个火山图吗
就是呢，当分组非常多的时候要展示，可以以这个barplot的形式去展示，省版面费，也非常直观。
然后他就找了一个数据GSE116439，作为学徒作业了。既然是学徒作业，为啥是我这个讲师来做呢。。。说来话长，那是一个月黑风高的晚上，我看到这张图发表了评论。。。
那。。。所以我就做了。
1.数据下载和分组信息获取
这个表达矩阵有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
分组已经很清楚的写在了title列，这是把他们拆分开成4列的结果。可以看到我们要的是浓度和时间的信息，在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矩阵，第一列是截距，无需理会。第二列往后，对应着1的行，说明表达矩阵里对应的列（样本）属于这个分组。比如第三个样本就属于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对比。后文画图就按照这种分组方式来。
还可以有别的比较方式，比如同时考虑这两个因素进行差异分析，就是说在做浓度比较时，是按照时间配对去比较的，同样，在做时间比较的时候，也是按照浓度配对去比较。两种方式任选咯。
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~复制过来的代码完成差异分析和画图(颜色我掰回来了，强迫症不可能允许红色表示下调)
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 画同款图
哈，基础包代码是现成的，没有玩尽兴，所以用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.矩阵画图，得改成数据框，并且宽变长。
2.堆叠式直方图的叠放次序，是用因子水平规定的，不规定就会自动。
3.直方图头顶加数字，用summarise配group_by实现计算，冷知识：ggplot2不同图层可以使用不同的数据画。
4.最后一句代码实现颜色图例顺序逆转
写点闲话
这个周是我们线上直播课的间隙，又是豆豆办港澳通行证等待的一个星期，等同于两人一起放假。我们已经来了珠海两年，长隆海洋王国就在几公里外，愣是到现在才想起来去，可能就是传说中的家门口的景点永远不想去吧。天实在是太热了，户外活动基本告别，豆豆好不容易放假又不想浪费，所以直接就去了。长隆的单日门票400，年卡800，所以我们选了年卡哈哈。分享两张有意思的照片
愤怒鱼
笑脸鱼
祝大家笑口常开。另外下周一（8月2号）开始又是一轮新的数据挖掘和生信入门线上课，如果你的暑假还没过完，不妨来跟我学习咯。
网友评论