参考连接：
https://www.r-bloggers.com/self-organising-maps-for-customer-segmentation-using-r/
https://blog.csdn.net/xbinworld/article/details/50818803
https://blog.csdn.net/xbinworld/article/details/50826892
https://blog.csdn.net/xbinworld/article/details/50890900
https://blog.csdn.net/qq_40793975/article/details/82689955
http://cdmd.cnki.com.cn/Article/CDMD-10161-1016197235.htm
https://www.cnblogs.com/sylvanas2012/p/5117056.html
R包的使用
https://www.shanelynn.ie/self-organising-maps-for-customer-segmentation-using-r/
https://clarkdatalabs.github.io/soms/SOM_NBA
https://www.jianshu.com/p/dd9af17f189f
https://www.cnblogs.com/yaos/p/12739887.html
R包的使用说明书：https://cran.r-project.org/web/packages/kohonen/kohonen.pdf
文章题目：Single cell transcriptome analysis of human, marmoset and mouse embryos reveals common and divergent features of preimplantation development
########################################################
#-------------------------------------------------------
# Topic:重复一篇文章中的SOM结果图
# Author:Wang Haiquan
# Date:Tue Jun 09 23:26:11 2020
# Mail:mg1835020@smail.nju.edu.cn
#-------------------------------------------------------
########################################################
library(kohonen)
sample<-read.table("clipboard",sep = "\t",header = T)#将table4的表达矩阵连带基因名称复制到剪切板上
#write.csv(sample,"sample.csv")
dim(sample)
sample<-sample[!duplicated(sample$Gene),]
rownames(sample)<-sample$Gene
sample<-sample[,-1]
#标准化
sample<-as.matrix(t(scale(t(sample))))
sample[1:6,1:6]
names(sample)<-names(sample)
sample<-sample[!is.nan(rowSums(sample)),]
sample[1:6,1:6]
dim(sample)
# 定义网络的大小和形状
som_grid <- somgrid(xdim = 30, ydim=30, topo="hexagonal")#文中是30*30
# Train the SOM model!
system.time(som_model <- supersom(sample, grid=som_grid, keep.data = TRUE))#计算一下运行时间
plot(som_model, type = "changes")
plot(som_model,type = "codes",codeRendering = "segments")
som_cluster <- cutree(hclust(dist(as.matrix(as.data.frame(som_model$codes)))), 6)#使用层次聚类进行分类
add.cluster.boundaries(som_model,clustering =som_cluster )
#注意！！！code的顺序是从左到右，从下向上，所以要从左下角开始看
#验证1，选择code.vector>1.5为依据进行分类，在与作者的结果比较(作者文中的描述的阈值是1.5)
som_res_code<-som_model$codes[[1]][apply(som_model$codes[[1]],1,function(x){max(x)>1.5}),]
som_res_code<-som_res_code>1.5
som_res_code<-apply(som_res_code,1,function(x){colnames(som_res_code)[x]})
som_res_code
rownames(som_model$codes[[1]])
som_res_code_fina<-rep(1,length(rownames(som_model$codes[[1]])))
names(som_res_code_fina)<-rownames(som_model$codes[[1]])
som_res_code_fina[names(som_res_code_fina)%in%names(som_res_code)]<-som_res_code
som_res_code_fina
plot(som_model,type = "codes")
add.cluster.boundaries(som_model,som_res_code_fina)
par(mfcol=c(2,3))
plot(som_model,type="property",property = getCodes(som_model)[,1],palette.name = coolBlueHotRed)
add.cluster.boundaries(som_model,som_res_code_fina)
plot(som_model,type="property",property = getCodes(som_model)[,2],palette.name = coolBlueHotRed)
add.cluster.boundaries(som_model,som_res_code_fina)
plot(som_model,type="property",property = getCodes(som_model)[,3],palette.name = coolBlueHotRed)
add.cluster.boundaries(som_model,som_res_code_fina)
plot(som_model,type="property",property = getCodes(som_model)[,4],palette.name = coolBlueHotRed)
add.cluster.boundaries(som_model,som_res_code_fina)
plot(som_model,type="property",property = getCodes(som_model)[,5],palette.name = coolBlueHotRed)
add.cluster.boundaries(som_model,som_res_code_fina)
plot(som_model,type="property",property = getCodes(som_model)[,6],palette.name = coolBlueHotRed)
add.cluster.boundaries(som_model,som_res_code_fina)
#验证2，我这里选择Compacted.morula中的基因，将它的名字与原文作者的结果比较
library(stringr)
author_res<-read.csv("clipboard",header = F,sep = "\t")#将table4中作者对于Compacted.morula的阶段特异基因的预测结果复制到剪切板上
stage="Compacted.morula"
late.icm_cluster<-som_res_code_fina[som_res_code_fina==stage]
late.icm_cluster<-str_split(names(late.icm_cluster),"V",simplify = T)[,2]
late.icm_cluster<-as.numeric(late.icm_cluster)
late.icm_cluster_gene<-rownames(sample)[which(som_model$unit.classif%in%late.icm_cluster)]
author_res1<-as.character(author_res$V1[author_res$V2=="comp. MOR"])
library(UpSetR)
x<-fromList(list(author=author_res1,mine=late.icm_cluster_gene))
upset(x,text.scale = 2)
