用SOM package画一个还行(我认为)的图
SOM package 和kohonen package都可以做SOM cluster分析。kohonen的图形展示要漂亮很多,还有ggsom包进一步优化。我没找到SOM包的图形优化包,决定自己画一个。
SOM package
require(som)
SOMdata=data.frame(read.table("xxx.xls", header = TRUE))
data_som=SOMdata[,2:7] #这一步为了去掉rownames
#data_som.f <- filtering(as.matrix(data_som)) 这步
scaledata<-scale(data_som)
normaldata<-normalize(as.matrix(data_som),byrow = T)
data_som=som(normaldata,xdim=5, ydim=6, init="linear", alpha=NULL,
alphaType="inverse", neigh="gaussian", topol="rect")
plot(data_som)
plot(data_som)这一步会得到一张中心数据点趋势线图,但是我不想跟大家都一样,所以决定用ggplot2把数据重新画一遍。中心数据在xx$code中
Make datasheet for ggplot2
得到的数据list(data_som)不能直接用ggplot2画图,必须转化为dataframe
#make a dataframe for ggplot2
data_som_frame<- as.data.frame(data_som$code)
colnames(data_som_frame)<-c("A", "B", "C", "D", "E", "F")
data_som_frame$sub<-c(1:30) #30是是SOM subclusters的数目= xdim × ydim
data_som_frame$nobs<-data_som$code.sum$nobs #nobs是每个subcluster中gene的个数
#reshape datasheet into long-sheet that ggplot2 can handle
require(reshape2)
data_long<- melt(data_som_frame, id = c("sub","nobs"), variable.name = "condition", value.name = "code")
然后就可以用ggplot2想画什么样画什么样啦
p<-ggplot(data_long, aes(x=condition, y=code, group= sub))+geom_point(pch = 1,
size = 2)+ scale_shape(solid=FALSE)+geom_line()
p<-p+facet_wrap(~sub,labeller = "label_both")
labeller
其实到这一步就可以了,会优化的继续优化。我不太会,我想把subcluster(1~30)标在图上,同时显示每个分组中gene的数目,意味着我需要在每一个facet的label上引用两个单元格的内容,另外也不想要原始的灰框,我采用了下面的方法:
p<-p+facet_wrap(~sub,labeller = "label_both")+theme(strip.background = element_blank())
p<-p+geom_text(aes(label = nobs), x = Inf, y = Inf,
hjust = 1.5, vjust = 3)
p
另外有一篇帖子 Put multi-variable facet_wrap labels on one line 这篇看起来非常简洁,也可以。
p2<-ggplot(data_long, aes(x=condition, y=code, group=sub)) +
geom_line()+ geom_point(size = 2) + facet_wrap(~paste(sub,nobs, sep="-"))
p2
最后用ggthemer 美化一下
require(ggthemes)
ggthemr('pale', layout = 'clean', spacing = 0.8, type = 'inner')
p
p2
最后是关于ggplot2和labeller()的使用,有机会应该仔细研读啊……画个图愁死了
更新:
用log2FoldChange后的数据画图
除了采用中心数据(xx$code),文献中很多图使用log2FoldChange来表示gene的变化,我也想用这种方法(而非code)来表示每个基因在每个条件下的变化,用一条平均数线来表示subcluster整体的变化,然后用hline()画出零点水平直线。
#SOM-graph using original data (log2foldchange)
#将分组以indices的方式赋回原数据组,data222.log是原值经过log2变换后的矩阵表
#data_som$visual$x[i]是SOM算法聚类后相应基因归属的组(x,y)
data222.log$x = c()
for(i in 1:2336) {data222.log$x[i] <-c(data_som$visual$x[i])}
data222.log$y = c()
for(i in 1:2336) {data222.log$y[i] <-c(data_som$visual$y[i])}
#reshape datasheet into long-sheet that ggplot2 can handle
data222.log.long<- melt(data222.log, id = c("x","y","gene"),
variable.name = "condition", value.name = "code")
#ggplot2+hline()
p<-ggplot(data222.log.long, aes(x=condition, y=code, group=gene)) +
geom_line(alpha = 3/5, colour = "#f3c57b")+geom_hline(aes(yintercept=0), colour="black")
#加入平均值,加粗
p+ stat_summary(aes(group=x+y), fun.y=mean, geom="line", colour="#db735c",size=1)
+ facet_wrap(~paste(x,y, sep="-"),scales = "free")
p