########################################################
#-------------------------------------------------------
# Topic:小鼠单细胞数据跑SOM流程
# Author:Wang Haiquan
# Date:Fri Jun 12 09:54:57 2020
# Mail:mg1835020@smail.nju.edu.cn
#-------------------------------------------------------
########################################################
library(kohonen)
library(Seurat)
library(stringr)
library(ggplot2)
library(pheatmap)
sample_mouse<-readRDS("../20200604human_mouse/data/tang_mouse.rds")
dim(sample_mouse)
sample_mouse_group<-str_split(colnames(sample_mouse),"_",simplify = T)[,2]
names(sample_mouse_group)<-colnames(sample_mouse)
sample_mouse_group<-sample_mouse_group[grep("^RS",sample_mouse_group)]
table(sample_mouse_group)
sample_mouse<-sample_mouse[,names(sample_mouse_group)]
sample_mouse<-CreateSeuratObject(sample_mouse)
sample_mouse<-NormalizeData(sample_mouse)%>%ScaleData()%>%FindVariableFeatures(nfeatures=5000)
sample_mouse<-sample_mouse@assays$RNA@scale.data[VariableFeatures(sample_mouse),]
dim(sample_mouse)
#-------------------------------------------------------
#Function:进行PCA及相关性检查
#-------------------------------------------------------
sample_mouse_pca<-prcomp(t(sample_mouse),scale. = F,center = F)
plot(sample_mouse_pca)
sample_mouse_pca<-sample_mouse_pca$x
sample_mouse_pca
ggplot(as.data.frame(sample_mouse_pca),aes(x=PC1,y=PC2,color=sample_mouse_group))+geom_point()
sample_anno_col<-data.frame(row.names = colnames(sample_mouse),cell_type=sample_mouse_group)
pheatmap(cor((sample_mouse)),show_rownames = F,show_colnames = F,annotation_col = sample_anno_col)
#-------------------------------------------------------
#Function:使用均值作为SOM的输入
#-------------------------------------------------------
sample_mouse<-readRDS("../20200604human_mouse/data/tang_mouse.rds")
dim(sample_mouse)
sample_mouse_group<-str_split(colnames(sample_mouse),"_",simplify = T)[,2]
names(sample_mouse_group)<-colnames(sample_mouse)
sample_mouse_group<-sample_mouse_group[grep("^RS",sample_mouse_group)]
table(sample_mouse_group)
sample_mouse<-sample_mouse[,names(sample_mouse_group)]
sample_mouse<-CreateSeuratObject(sample_mouse)
sample_mouse<-NormalizeData(sample_mouse)%>%ScaleData()%>%FindVariableFeatures(nfeatures=5000)
sample_mouse@active.ident<-as.factor(sample_mouse_group)
sample_mouse<-AverageExpression(sample_mouse,features = VariableFeatures(sample_mouse))
sample_mouse<-sample_mouse$RNA
dim(sample_mouse)
sample_mouse<-t(scale(t(sample_mouse)))
#做SOM
mouse_som_grid<-somgrid(xdim = 20, ydim=20, topo="hexagonal")
mouse_som<-supersom(sample_mouse,mouse_som_grid,rlen = 500)
par(mfcol=c(2,2))
plot(mouse_som,"changes")
plot(mouse_som,"counts",palette.name = coolBlueHotRed)
plot(mouse_som,"codes",palette.name = coolBlueHotRed)
plot(mouse_som,"quality",palette.name = coolBlueHotRed)
par(mfcol=c(1,2))
plot(mouse_som,"mapping",palette.name = coolBlueHotRed)
plot(mouse_som,"dist.neighbours",palette.name = coolBlueHotRed)
par(mfcol=c(2,2))
plot(mouse_som,"property",property = getCodes(mouse_som)[,1],palette.name = coolBlueHotRed,main=colnames(getCodes(mouse_som))[1])
plot(mouse_som,"property",property = getCodes(mouse_som)[,2],palette.name = coolBlueHotRed,main=colnames(getCodes(mouse_som))[2])
plot(mouse_som,"property",property = getCodes(mouse_som)[,3],palette.name = coolBlueHotRed,main=colnames(getCodes(mouse_som))[3])
plot(mouse_som,"property",property = getCodes(mouse_som)[,4],palette.name = coolBlueHotRed,main=colnames(getCodes(mouse_som))[4])
#进行簇的划分
#定义1为阈值，其它cluster为1
mouse_code<-getCodes(mouse_som)
mouse_code_cluster<-t(apply(mouse_code,1,function(x){ifelse(x>1,x,0)}))
mouse_code_cluster<-apply(mouse_code_cluster,1,function(x){ifelse(sum(x)>1,colnames(mouse_code)[which(x==max(x))],5)})
mouse_code_cluster
par(mfcol=c(2,2))
plot(mouse_som,"property",property = getCodes(mouse_som)[,1],palette.name = coolBlueHotRed,main=colnames(getCodes(mouse_som))[1])
add.cluster.boundaries(mouse_som,mouse_code_cluster)
plot(mouse_som,"property",property = getCodes(mouse_som)[,2],palette.name = coolBlueHotRed,main=colnames(getCodes(mouse_som))[2])
add.cluster.boundaries(mouse_som,mouse_code_cluster)
plot(mouse_som,"property",property = getCodes(mouse_som)[,3],palette.name = coolBlueHotRed,main=colnames(getCodes(mouse_som))[3])
add.cluster.boundaries(mouse_som,mouse_code_cluster)
plot(mouse_som,"property",property = getCodes(mouse_som)[,4],palette.name = coolBlueHotRed,main=colnames(getCodes(mouse_som))[4])
add.cluster.boundaries(mouse_som,mouse_code_cluster)
#做热图看差异基因的结果
mouse_som_genelist<-data.frame(gene=rownames(sample_mouse),cluster1=mouse_som$unit.classif)
mouse_code_df<-data.frame(cluster1=str_split(names(mouse_code_cluster),"V",simplify = T)[,2],cluster2=mouse_code_cluster)
mouse_som_genelist<-merge(mouse_som_genelist,mouse_code_df,by="cluster1")
mouse_som_genelist<-mouse_som_genelist[order(factor(mouse_som_genelist$cluster2)),]
rownames(mouse_som_genelist)<-mouse_som_genelist$gene
head(mouse_som_genelist)
pheatmap(sample_mouse[as.character(mouse_som_genelist$gene),],
cluster_rows = F,cluster_cols = F,show_rownames = F,show_colnames = F,
annotation_row =mouse_som_genelist[,-2])
sample_mouse_scale<-readRDS("../20200604human_mouse/data/tang_mouse.rds")
sample_mouse_scale<-sample_mouse_scale[rownames(mouse_som_genelist),names(sample_mouse_group)]
sample_mouse_scale<-CreateSeuratObject(sample_mouse_scale)
sample_mouse_scale<-NormalizeData(sample_mouse_scale)%>%ScaleData()
sample_mouse_scale<-sample_mouse_scale@assays$RNA@scale.data
head(sample_mouse_scale)
pheatmap(sample_mouse_scale[,order(factor(sample_mouse_group))],
show_rownames = F,show_colnames = F,
cluster_rows = F,cluster_cols = F,
annotation_row = mouse_som_genelist[,-2],
annotation_col = sample_anno_col,
color = colorRampPalette(c("purple","black","yellow"))(100),
breaks = seq(-4,4,8/100))
#与findmark的比较
sample_mouse_scale<-readRDS("../20200604human_mouse/data/tang_mouse.rds")
sample_mouse_scale<-sample_mouse_scale[rownames(mouse_som_genelist),names(sample_mouse_group)]
sample_mouse_scale<-CreateSeuratObject(sample_mouse_scale)
sample_mouse_scale@active.ident<-as.factor(sample_mouse_group)
sample_mouse_scale<-NormalizeData(sample_mouse_scale)%>%ScaleData()
sample_mouse_find_marker<-FindAllMarkers(sample_mouse_scale,only.pos = T,test.use = "t")
sample_mouse_find_marker[1:6,]
library(UpSetR)
sample_seurat_marker_res<-sample_mouse_find_marker[sample_mouse_find_marker$p_val<0.05,]
sample_seurat_marker_res<-tapply(sample_seurat_marker_res$gene,sample_seurat_marker_res$cluster,print)
sample_som_marker_res<-tapply(as.character(mouse_som_genelist$gene),mouse_som_genelist$cluster2,print)
names(sample_seurat_marker_res)
names(sample_som_marker_res)
a=4
sample_mark_compare<-fromList(list(seurat_1=sample_seurat_marker_res[[a]],
som_1=sample_som_marker_res[[a+1]]))
upset(sample_mark_compare)
sample_mouse_scale<-readRDS("../20200604human_mouse/data/tang_mouse.rds")
sample_mouse_scale<-sample_mouse_scale[rownames(mouse_som_genelist),names(sample_mouse_group)]
sample_mouse_scale<-CreateSeuratObject(sample_mouse_scale)%>%NormalizeData()%>%FindVariableFeatures%>%ScaleData()%>%
RunPCA()%>%RunUMAP(dim=1:10)
sample_mouse_scale@active.ident<-as.factor(sample_mouse_group)
VlnPlot(sample_mouse_scale,sample_seurat_marker_res$RS1o2[1:10])
#########################################################################
#结论，感觉效果挺差的，如果用来处理组内均一性很好的样本的话应该效果会不错但是用来处理单细胞，感觉emmmm
#########################################################################
修改一下上面错误的结论
通过如下代码的更细致的观察
发现在处理单细胞的数据中，使用SOM方法找到的模式表达基因更倾向于特异性，缺点是会受一些低表达基因异常值的影响，而错将低表达基因也纳入特异性表达基因列表，因此使用前需要对低表达基因做细致筛查。（简单看了一下基本上som单独筛选的都是低表达的基因，但是在热图展示上由于是根据行scale的，因此很难看出来，所以我们需要做之前细致筛查！）
普通的差异分析方法优点是能够过滤低表达基因，但是在模式的特异性上可能比较差的。涉及到至少有一组显著就显著还是唯一显著这样的问题，寻找差异基因时一般是满足前者即可，而我们在寻找阶段性特异表达基因的时候都希望是唯一显著！
#比较两种方法筛到基因的异同
sample_mouse_scale<-readRDS("../20200604human_mouse/data/tang_mouse.rds")
sample_mouse_scale<-sample_mouse_scale[rownames(mouse_som_genelist),names(sample_mouse_group)]
sample_mouse_scale<-CreateSeuratObject(sample_mouse_scale)%>%NormalizeData()%>%FindVariableFeatures%>%ScaleData()%>%
RunPCA()%>%RunUMAP(dim=1:10)
sample_mouse_scale@active.ident<-as.factor(sample_mouse_group)
#我们以二者结果相差最大的RS3o4组作为比较
#二者筛到的基因
VlnPlot(sample_mouse_scale,intersect(sample_som_marker_res$RS3o4,sample_seurat_marker_res$RS3o4)[sample(1:50,5)])
#som筛到 findmarker没筛到
VlnPlot(sample_mouse_scale,setdiff(sample_som_marker_res$RS3o4,sample_seurat_marker_res$RS3o4)[sample(1:50,5)])
#som没有筛到，findmaker筛到了
VlnPlot(sample_mouse_scale,setdiff(sample_seurat_marker_res$RS3o4，sample_som_marker_res$RS3o4)[sample(1:50,5)])
两种方法筛到基因的差别
二者都筛到的基因
仅som筛到
仅findmarke筛到
寻找阶段特异性基因，可以通过设定codebook vector阈值的方式查找，从下图我们可以看到，在第30个神经元中，四类细胞的值分别为1.4840752 -0.3065276 -0.6562483 -0.5212993，因为该神经元下的基因的表达情况应尽量与该权重向量靠近，因此我们有理由推断该神经元下的基因在第一类细胞中表达很高，在其它三类细胞中表达很低，如下图所示
Rplot05.png
Rplot06.png
