7_Addmodule和基因加和法add 得到的细胞类型打分在空间上空转对比

7_Addmodule和基因加和法add 得到的细胞类型打分在空间上空转对比

library(ggplot2)
library(ggalluvial)
library(svglite)
library(Seurat)
library(openxlsx)
library(tibble)
library(ggpubr)



path = "G:/silicosis/sicosis/silicosis_ST/overlapped_map/contrast_between_addmodule_and addtion"
dir.create(path)
setwd(path)
getwd()

marker=read.xlsx('G:/silicosis/sicosis/yll/macrophage/no cluster2/0.3/pure_cluster01345_dotplot/cluster01345.xlsx')
#marker=read.xlsx("D:/Win10 System/Desktop/singlecell _marker.xlsx",sheet=1)
head(marker)
#names(marker) = c("B cell", "AT", "Fib", "Endo", "neutrophil", "Cycling base", "T", "Mono", "Macro")
colnames(marker)


load("G:/silicosis/需求/矽肺-数据分析结果-0308_yll/矽肺-数据分析结果-0308/NS_7_sct.rds")
load("G:/silicosis/需求/矽肺-数据分析结果-0308_yll/矽肺-数据分析结果-0308/NS_56_sct.rds")
load("G:/silicosis/需求/矽肺-数据分析结果-0308_yll/矽肺-数据分析结果-0308/sio2_7_sct.rds")
load("G:/silicosis/需求/矽肺-数据分析结果-0308_yll/矽肺-数据分析结果-0308/sio2_56_sct.rds")

load()


library(dplyr)
library(Hmisc)

#以ciliated 为例子
marker
names(marker)

marker$`cluster1_AM`
cellname="cluster1_AM"

mymarker=na.omit(unique(marker$"cluster1_AM")) %>% capitalize() %>% as.character() %>%list()
mymarker
length(mymarker) #查看一共有几个基因



length(mymarker) #查看一共有几个基因
#或者直接使用随机函数来取
number=c(1:5)  #命名用的，根据取样的大小自己设置
number
cellname  #命名用的
mymarker

library(ggpubr)
library(ggplot2)

#单纯的基因之和加起来
for (i in 1:6) {
     #循环的次数最好由总的mymarker 数目length（mymarkers）决定
  markers=sample(mymarker,size=5)   #取样大小必须小于length（mymarkers）
  #markers=mymarker[i]
  number=c(1:length(markers))
  print(paste0("第 ",i," 次:",markers,collapse=";"))
  if (1==1){
    
    #ns_7d
    #markers = rownames(NS_7_sct)[number]## input markes
    expr = colMeans(NS_7_sct@assays$SCT@scale.data[markers, ])
    NS_7_sct@assays$SCT@scale.data[markers[1], ] = expr
    p1=SpatialFeaturePlot(NS_7_sct, features = markers[1], slot = "scale.data")+ 
      ggtitle(paste(markers, collapse = "|"))
    
    #ns_56d
    #markers = rownames(NS_56_sct)[number]## input markes
    expr = colMeans(NS_56_sct@assays$SCT@scale.data[markers, ])
    NS_56_sct@assays$SCT@scale.data[markers[1], ] = expr
    p2 = SpatialFeaturePlot(NS_56_sct, features = markers[1], slot = "scale.data") +
      ggtitle(paste(markers, collapse = "|"))
    
    #sio2_7d
    #markers = rownames(sio2_7_sct)[number]## input markes
    expr = colMeans(sio2_7_sct@assays$SCT@scale.data[markers, ])
    sio2_7_sct@assays$SCT@scale.data[markers[1], ] = expr
    p3 = SpatialFeaturePlot(sio2_7_sct, features = markers[1], slot = "scale.data")+ 
      ggtitle(paste(markers, collapse = "|"))
    
    #sio2_56d
    #markers = rownames(sio2_56_sct)[number]## input markes
    expr = colMeans(sio2_56_sct@assays$SCT@scale.data[markers, ])
    sio2_56_sct@assays$SCT@scale.data[markers[1], ] = expr
    p4 = SpatialFeaturePlot(sio2_56_sct, features = markers[1], slot = "scale.data")+ 
      ggtitle(paste(markers, collapse = "|"))
    #ggarrange(p1,p2,p3,p4,ncol = 2,nrow =2)
    
    jpeg(paste0(cellname,"_","total_",length(mymarker),"_",paste0(min(number),"-",max(number)), 
                paste(markers,collapse = "_"),"_.jpeg"),
         height = 12, width = 12, units = 'in', res=300)
    p=ggarrange(p1,p2,p3,p4,ncol = 2,nrow =2)
    print(p)
    dev.off()
  }  
}



getwd()
#基因之和加起来 addmodule
DefaultAssay(NS_7_sct)
for (i in 1) {
     #循环的次数最好由总的mymarker 数目length（mymarkers）决定
  markers=sample(mymarker,size=6,replace = FALSE)   #取样大小必须小于length（mymarkers）
  #markers=mymarker[i]
  number=c(1:length(markers))
  print(paste0("第 ",i," 次:",markers,collapse=";"))
  
  
  
  if (1==1){
    
    #ns_7d
    #markers = rownames(NS_7_sct)[number]## input markes
    expr = colMeans(NS_7_sct@assays$SCT@scale.data[markers, ])
    NS_7_sct@assays$SCT@scale.data[markers[1], ] = expr
    p1=SpatialFeaturePlot(NS_7_sct, features = markers[1], slot = "scale.data")+ 
      ggtitle(paste(markers, collapse = "|"))
    
#ns_7d addmodule 
if(1==1){
    
  markerlist=list(markers)
  NS_7_sct=AddModuleScore(NS_7_sct,
                          features =markerlist,
                          name = "markerlist")
  #结果保存在这里
  colnames(NS_7_sct@meta.data)
  head(NS_7_sct@meta.data)
  colnames(NS_7_sct@meta.data)[7]<-"markerlist"
  
  p1addmodule=SpatialFeaturePlot(NS_7_sct, features = "markerlist", slot = "scale.data")+
    ggtitle(paste(unlist(markerlist), collapse = "|"))
  
  NS_7_sct@meta.data<-NS_7_sct@meta.data[,-7]
}
    
    
    
    #ns_56d
    #markers = rownames(NS_56_sct)[number]## input markes
    expr = colMeans(NS_56_sct@assays$SCT@scale.data[markers, ])
    NS_56_sct@assays$SCT@scale.data[markers[1], ] = expr
    p2 = SpatialFeaturePlot(NS_56_sct, features = markers[1], slot = "scale.data") +
      ggtitle(paste(markers, collapse = "|"))

    
#ns_56d addmodule
if(1==1){
    
  markerlist=list(markers)
  NS_56_sct=AddModuleScore(NS_56_sct,
                           features =markerlist,
                           name = "markerlist")
  #结果保存在这里
  colnames(NS_56_sct@meta.data)
  head(NS_56_sct@meta.data)
  colnames(NS_56_sct@meta.data)[7]<-"markerlist"
  
  p2addmodule=SpatialFeaturePlot(NS_56_sct, features = "markerlist", slot = "scale.data")
  NS_56_sct@meta.data<-NS_56_sct@meta.data[,-7]
  }
    
    #sio2_7d
    #markers = rownames(sio2_7_sct)[number]## input markes
    expr = colMeans(sio2_7_sct@assays$SCT@scale.data[markers, ])
    sio2_7_sct@assays$SCT@scale.data[markers[1], ] = expr
    p3 = SpatialFeaturePlot(sio2_7_sct, features = markers[1], slot = "scale.data")+ 
      ggtitle(paste(markers, collapse = "|"))
    
#sio2_7d addmodule 
if(1==1){
       
markerlist=list(markers)
sio2_7_sct=AddModuleScore(sio2_7_sct,
                          features =markerlist,
                          name = "markerlist")
#结果保存在这里
colnames(sio2_7_sct@meta.data)
head(sio2_7_sct@meta.data)
colnames(sio2_7_sct@meta.data)[7]<-"markerlist"

p3addmodule=SpatialFeaturePlot(sio2_7_sct, features = "markerlist", slot = "scale.data")+
  ggtitle(paste(cellname)) 
sio2_7_sct@meta.data<-sio2_7_sct@meta.data[,-7]

}
       
    
    
    
    #sio2_56d
    #markers = rownames(sio2_56_sct)[number]## input markes
    expr = colMeans(sio2_56_sct@assays$SCT@scale.data[markers, ])
    sio2_56_sct@assays$SCT@scale.data[markers[1], ] = expr
    p4 = SpatialFeaturePlot(sio2_56_sct, features = markers[1], slot = "scale.data")+ 
      ggtitle(paste(markers, collapse = "|"))
    #ggarrange(p1,p2,p3,p4,ncol = 2,nrow =2)
  
#sio2_56d addmodule
    if(1==1){
    
      markerlist=list(markers)
      sio2_56_sct=AddModuleScore(sio2_56_sct,
                                features =markerlist,
                                name = "markerlist")
      #结果保存在这里
      colnames(sio2_56_sct@meta.data)
      head(sio2_56_sct@meta.data)
      colnames(sio2_56_sct@meta.data)[7]<-"markerlist"
      
      p4addmodule=SpatialFeaturePlot(sio2_56_sct, features = "markerlist", slot = "scale.data")+
        ggtitle(paste(unlist(markerlist), collapse = "|")) 
      sio2_56_sct@meta.data<-sio2_56_sct@meta.data[,-7]
      }
    
    
    #addition
    jpeg(paste0(cellname,"_","total_",length(mymarker),"_",paste0(min(number),"-",max(number)), 
                paste(markers,collapse = "_"),"_.jpeg"),
         height = 12, width = 12, units = 'in', res=300)
    p=ggarrange(p1,p2,p3,p4,ncol = 2,nrow =2)
    print(p)
    dev.off()
    
    
    #addmodule
    jpeg(paste0(cellname,"_","total_",length(unlist(markerlist)),"_",paste0(min(number),"-",max(number)), 
                paste(unlist(markerlist)[1:15],collapse = "_"),"_.jpeg"), #只取前15个
         height = 12, width = 12, units = 'in', res=600)
    p=ggpubr::ggarrange(p1addmodule,p2addmodule,p3addmodule,p4addmodule,ncol = 2,nrow =2)
    print(p)
    dev.off()
    
    
  }  
}






while (dev.off()) {
    
  dev.off()
}
getwd()
#判断基因是否在空转内
for (each in mymarker) {
    
  if (each %in% c(rownames(NS_7_sct@assays$SCT@scale.data),
                  rownames(NS_56_sct@assays$SCT@scale.data),
                  rownames(sio2_7_sct@assays$SCT@scale.data),
                  rownames(sio2_56_sct@assays$SCT@scale.data)
  )==FALSE)
    print(paste0(each,"未在空转数据检测到"))
  
  if (each %in% c(rownames(NS_7_sct@assays$SCT@scale.data))==FALSE)
    print(paste0(each,"未在NS_7空转数据检测到"))
  if (each %in% c(rownames(NS_56_sct@assays$SCT@scale.data))==FALSE)
    print(paste0(each,"未在NS_56空转数据检测到"))
  if (each %in% c(rownames(sio2_7_sct@assays$SCT@scale.data))==FALSE)
    print(paste0(each,"未在sio2_7空转数据检测到"))
  if (each %in% c(rownames(sio2_56_sct@assays$SCT@scale.data))==FALSE)
    print(paste0(each,"未在sio2_56空转数据检测到"))
} #判断基因是否在空转内