用SOM package畫一個還行(我認為)的圖
SOM package 和kohonen package都可以做SOM cluster分析。kohonen的圖形展示要漂亮很多,還有ggsom包進一步優(yōu)化。我沒找到SOM包的圖形優(yōu)化包,決定自己畫一個。
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)這一步會得到一張中心數(shù)據(jù)點趨勢線圖,但是我不想跟大家都一樣,所以決定用ggplot2把數(shù)據(jù)重新畫一遍。中心數(shù)據(jù)在xx$code中
Make datasheet for ggplot2
得到的數(shù)據(jù)list(data_som)不能直接用ggplot2畫圖,必須轉(zhuǎn)化為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的數(shù)目= xdim × ydim
data_som_frame$nobs<-data_som$code.sum$nobs #nobs是每個subcluster中gene的個數(shù)
#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
其實到這一步就可以了,會優(yōu)化的繼續(xù)優(yōu)化。我不太會,我想把subcluster(1~30)標在圖上,同時顯示每個分組中gene的數(shù)目,意味著我需要在每一個facet的label上引用兩個單元格的內(nèi)容,另外也不想要原始的灰框,我采用了下面的方法:
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后的數(shù)據(jù)畫圖
除了采用中心數(shù)據(jù)(xx$code),文獻中很多圖使用log2FoldChange來表示gene的變化,我也想用這種方法(而非code)來表示每個基因在每個條件下的變化,用一條平均數(shù)線來表示subcluster整體的變化,然后用hline()畫出零點水平直線。
#SOM-graph using original data (log2foldchange)
#將分組以indices的方式賦回原數(shù)據(jù)組,data222.log是原值經(jīng)過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