mia kaiti

9个cluster

#request 2
.libPaths(c( "/home/data/t040413/R/x86_64-pc-linux-gnu-library/4.2",
             "/home/data/t040413/R/yll/usr/local/lib/R/site-library",  
             "/home/data/refdir/Rlib/", "/usr/local/lib/R/library"))




setwd('/home/data/t040413/silicosis/spatial/')




library(Seurat)
library(dplyr)
#d.all=readRDS("~/silicosis/spatial/integrated_slides/integrated_slides.rds")
load("/home/data/t040413/silicosis/spatial_transcriptomics/silicosis_ST_harmony_SCT_r0.5.rds")
dim(d.all)
DefaultAssay(d.all)="Spatial"
#visium_slides=SplitObject(object = d.all,split.by = "stim")

names(d.all);dim(d.all)
d.all@meta.data %>%head()
head(colnames(d.all))
#给d.all 添加meta信息------
adata_obs=read.csv("~/silicosis/spatial/adata_obs.csv")
head(adata_obs)
mymeta=  paste0(d.all@meta.data$orig.ident,"_",colnames(d.all))  %>% gsub("-.*","",.)  # %>%  head()
head(mymeta)
tail(mymeta)

#掉-及其之后内容
adata_obs$col= adata_obs$spot_id %>% gsub("-.*","",.)    # %>%  head()
head(adata_obs)

rownames(adata_obs)=adata_obs$col
adata_obs=adata_obs[mymeta,]
head(adata_obs)
identical(mymeta,adata_obs$col)
d.all=AddMetaData(d.all,metadata = adata_obs)

DefaultAssay(d.all)='Spatial'
DimPlot(d.all,group.by = 'stim')
DimPlot(d.all,label = TRUE)
dev.off()


####st-----
cluster_gene_stat=FindAllMarkers(d.all,only.pos = TRUE,logfc.threshold = 0.5)
head(cluster_gene_stat)
table(cluster_gene_stat$cluster)

# Initialize a dataframe for us to store values in:
st.clusts=paste0("cluster",seq(0,8,1))


head(cluster_gene_stat)
st.marker.list=split(cluster_gene_stat,f = cluster_gene_stat$cluster)
st.marker.list = lapply(st.marker.list,rownames)
head(st.marker.list)
#st.clusts=lapply(st.marker.list,rownames)

st.clusts
st.marker.list


###sc----
cellType_gene_stat =readRDS("~/silicosis/spatial2/marker_list_silicosis_nogrem1.rds")
head(cellType_gene_stat)
table(cellType_gene_stat$cluster)

cellType_marker = list()
for(i in unique(cellType_gene_stat$cluster)){ #按照标准 获得20个细胞类型的marker基因 列表
  cellType_marker[[i]] = cellType_gene_stat[cellType_gene_stat$cluster==i & cellType_gene_stat$p_val_adj<0.1 &cellType_gene_stat$avg_log2FC>0.7, "gene"]
}
cellType_marker
head(cellType_marker)
names(cellType_marker)
sc.marker.list=cellType_marker



sc.clusts=names(cellType_marker)

sc.clusts
cellType_marker



st.clusts
st.marker.list


N <- length(st.clusts) 
M <- length(sc.clusts)
MIA.results <- matrix(0,nrow = M, ncol = N)
row.names(MIA.results) <- sc.clusts
colnames(MIA.results) <- st.clusts
head(MIA.results)




# Gene universe
gene.universe <- length(rownames(cluster_gene_stat))



# Loop over ST clusters
for (i in 1:N) {
  # Then loop over SC clusters
  # i=1
  for (j in 1:M) {
    genes1 <- st.marker.list[[st.clusts[i]]]
    genes2 <- sc.marker.list[[sc.clusts[j]]]
    
    # Hypergeometric    
    A <- length(intersect(genes1,genes2))
    B <- length(genes1)
    C <- length(genes2)
    enr <- -log10(phyper(A, B, gene.universe-B, C, lower.tail = FALSE))
    dep <- -log10(1-phyper(A, B, gene.universe-B, C, lower.tail = FALSE))
    if (enr < dep) {
      MIA.results[j,i] = -dep
    } else {
      MIA.results[j,i] = enr
    }
  }
}


# Some results were -Inf...check why this is the case...
MIA.results[is.infinite(MIA.results)] <- 0
MIA.results


pheatmap::pheatmap(MIA.results,scale ='row' )


MIA.results2=MIA.results %>%as.matrix()

MIA.results2[ MIA.results >5 ] =5


pheatmap::pheatmap(MIA.results2,scale ='row' )


# Visualize as heatmap
library(reshape2)
library(ggplot2)
library(dplyr)
library(Seurat)
heatmap_df <- data.frame('Cell types' = melt(MIA.results)[,1],
                         'Tissue regions' = melt(MIA.results)[,2],
                         enrichment = melt(MIA.results)[,3])
range(heatmap_df$enrichment)

head(heatmap_df)

boxplot(heatmap_df$enrichment)


heatmap_df[abs(heatmap_df$enrichment) >7,]$enrichment=7


ggplot(data = heatmap_df, aes(x = Tissue.regions, y = Cell.types, fill = enrichment)) +
  geom_tile() + 
  scale_fill_gradient2(low = "navyblue", high = "red", mid = "white", 
                       midpoint = 0, limit = c(-7,7), space = "Lab", 
                       name="Enrichment \n -log10(p)") +
  ylim(heatmap_df$Cell.types %>% levels() %>% sort() %>% rev())+
  theme_minimal()+RotatedAxis() 


library(reshape2)
library(ggplot2)
library(dplyr)
library(Seurat)
library(scales)
heatmap_df <- data.frame('Cell types' = melt(MIA.results)[,1],
                         'Tissue regions' = melt(MIA.results)[,2],
                         enrichment = melt(MIA.results)[,3])
head(heatmap_df)

ggplot(data = heatmap_df, aes(x = Tissue.regions, y = Cell.types, fill = enrichment)) +
  geom_tile() + 
  scale_fill_gradient2(low = "blue", high = "red",mid = 'white', 
                       midpoint = 0, limit = c(-1,1), space = "Lab",oob=squish,
                       name="Enrichment \n -log10(p)") +
  ylim(heatmap_df$Cell.types %>% levels() %>% sort() %>% rev())+
  theme_minimal()+RotatedAxis() 

四个区

#request 2
.libPaths(c( "/home/data/t040413/R/x86_64-pc-linux-gnu-library/4.2",
             "/home/data/t040413/R/yll/usr/local/lib/R/site-library",  
             "/home/data/refdir/Rlib/", "/usr/local/lib/R/library"))




setwd('/home/data/t040413/silicosis/spatial/')




library(Seurat)
library(dplyr)
#d.all=readRDS("~/silicosis/spatial/integrated_slides/integrated_slides.rds")
load("/home/data/t040413/silicosis/spatial_transcriptomics/silicosis_ST_harmony_SCT_r0.5.rds")
dim(d.all)
DefaultAssay(d.all)="Spatial"
#visium_slides=SplitObject(object = d.all,split.by = "stim")

names(d.all);dim(d.all)
d.all@meta.data %>%head()
head(colnames(d.all))
#给d.all 添加meta信息------
adata_obs=read.csv("~/silicosis/spatial/adata_obs.csv")
head(adata_obs)
mymeta=  paste0(d.all@meta.data$orig.ident,"_",colnames(d.all))  %>% gsub("-.*","",.)  # %>%  head()
head(mymeta)
tail(mymeta)

#掉-及其之后内容
adata_obs$col= adata_obs$spot_id %>% gsub("-.*","",.)    # %>%  head()
head(adata_obs)

rownames(adata_obs)=adata_obs$col
adata_obs=adata_obs[mymeta,]
head(adata_obs)
identical(mymeta,adata_obs$col)
d.all=AddMetaData(d.all,metadata = adata_obs)

DefaultAssay(d.all)='Spatial'
DimPlot(d.all,group.by = 'stim')
DimPlot(d.all,label = TRUE)
dev.off()

Idents(d.all)=d.all$clusters

####st-----
cluster_gene_stat=FindAllMarkers(d.all,only.pos = TRUE,logfc.threshold = 0.7)
head(cluster_gene_stat)
table(cluster_gene_stat$cluster)



head(cluster_gene_stat)
st.marker.list=split(cluster_gene_stat,f = cluster_gene_stat$cluster)
st.marker.list = lapply(st.marker.list,rownames)
head(st.marker.list)

# Initialize a dataframe for us to store values in:
st.clusts=names(st.marker.list)



st.clusts
st.marker.list


###sc----
cellType_gene_stat =readRDS("~/silicosis/spatial2/marker_list_silicosis_nogrem1.rds")
head(cellType_gene_stat)
table(cellType_gene_stat$cluster)

cellType_marker = list()
for(i in unique(cellType_gene_stat$cluster)){ #按照标准 获得20个细胞类型的marker基因 列表
  cellType_marker[[i]] = cellType_gene_stat[cellType_gene_stat$cluster==i & cellType_gene_stat$p_val_adj<0.1 &cellType_gene_stat$avg_log2FC>0.7, "gene"]
}
cellType_marker
head(cellType_marker)
names(cellType_marker)
sc.marker.list=cellType_marker



sc.clusts=names(cellType_marker)

sc.clusts
cellType_marker



st.clusts
st.marker.list


N <- length(st.clusts) 
M <- length(sc.clusts)
MIA.results <- matrix(0,nrow = M, ncol = N)
row.names(MIA.results) <- sc.clusts
colnames(MIA.results) <- st.clusts
head(MIA.results)




# Gene universe
gene.universe <- length(rownames(cluster_gene_stat))



# Loop over ST clusters
for (i in 1:N) {
  # Then loop over SC clusters
  # i=1
  for (j in 1:M) {
    genes1 <- st.marker.list[[st.clusts[i]]]
    genes2 <- sc.marker.list[[sc.clusts[j]]]
    
    # Hypergeometric    
    A <- length(intersect(genes1,genes2))
    B <- length(genes1)
    C <- length(genes2)
    enr <- -log10(phyper(A, B, gene.universe-B, C, lower.tail = FALSE))
    dep <- -log10(1-phyper(A, B, gene.universe-B, C, lower.tail = FALSE))
    if (enr < dep) {
      MIA.results[j,i] = -dep
    } else {
      MIA.results[j,i] = enr
    }
  }
}


# Some results were -Inf...check why this is the case...
MIA.results[is.infinite(MIA.results)] <- 0
MIA.results


pheatmap::pheatmap(MIA.results,scale ='row' )


MIA.results2=MIA.results %>%as.matrix()

MIA.results2[ MIA.results >5 ] =5


pheatmap::pheatmap(MIA.results2,scale ='row' )


# Visualize as heatmap
library(reshape2)
library(ggplot2)
library(dplyr)
library(Seurat)
heatmap_df <- data.frame('Cell types' = melt(MIA.results)[,1],
                         'Tissue regions' = melt(MIA.results)[,2],
                         enrichment = melt(MIA.results)[,3])
range(heatmap_df$enrichment)

head(heatmap_df)

boxplot(heatmap_df$enrichment)


heatmap_df[abs(heatmap_df$enrichment) >7,]$enrichment=7


ggplot(data = heatmap_df, aes(x = Tissue.regions, y = Cell.types, fill = enrichment)) +
  geom_tile() + 
  scale_fill_gradient2(low = "navyblue", high = "red", mid = "white", 
                       midpoint = 0, limit = c(-7,7), space = "Lab", 
                       name="Enrichment \n -log10(p)") +
  ylim(heatmap_df$Cell.types %>% levels() %>% sort() %>% rev())+
  theme_minimal()+RotatedAxis() 


library(reshape2)
library(ggplot2)
library(dplyr)
library(Seurat)
library(scales)
heatmap_df <- data.frame('Cell types' = melt(MIA.results)[,1],
                         'Tissue regions' = melt(MIA.results)[,2],
                         enrichment = melt(MIA.results)[,3])
head(heatmap_df)

ggplot(data = heatmap_df, aes(x = Tissue.regions, y = Cell.types, fill = enrichment)) +
  geom_tile() + 
  scale_fill_gradient2(low = "blue", high = "red",mid = 'white', 
                       midpoint = 0, limit = c(-1,1), space = "Lab",oob=squish,
                       name="Enrichment \n -log10(p)") +
  ylim(heatmap_df$Cell.types %>% levels() %>% sort() %>% rev())+
  theme_minimal()+RotatedAxis()