TME髓系亚群细分“一本通”

# Step1.读入数据 --------------------------------------------------------------
### 1.读入数据
seurat.data = readRDS("./GSE266919_Myeloid.rds")
seurat.data

# 计算 logcounts（对数归一化）
logcounts(seurat.data) <- log1p(counts(seurat.data))  # 使用 log(x + 1) 转换

seurat.obj <- as.Seurat(seurat.data,counts = "counts",data = "logcounts")
seurat.obj

# Step2.三步走：标准化、特征选择和归一化分析 ------------------------------------------------
### 2.三步走：标准化、特征选择和归一化分析
seurat.obj <- seurat.obj %>% NormalizeData(verbose = F) %>%
  FindVariableFeatures(selection.method = "vst", nfeatures = 2000, verbose = F) %>% 
  ScaleData(verbose = F)

seurat.obj


# Step3.降维和聚类，检查批次 --------------------------------------------------------
### 3.1 降维和聚类
seurat.obj = seurat.obj %>% 
  RunPCA(npcs = 30, verbose = F) %>% 
  #RunTSNE(reduction = "pca", dims = 1:30, verbose = F) %>% 
  RunUMAP(reduction = "pca", dims = 1:30, verbose = F)

seurat.obj

### 3.2 检查批次
options(repr.plot.width = 10, repr.plot.height = 4.5)
p1.compare=wrap_plots(ncol = 2,
                      DimPlot(seurat.obj, reduction = "pca", group.by = "Sample")+NoAxes()+ggtitle("Before_PCA"),
                      DimPlot(seurat.obj, reduction = "umap", group.by = "Sample")+NoAxes()+ggtitle("Before_UMAP"),
                      guides = "collect"
)
p1.compare

### 3.3 若无批次，继续运行聚类
seurat.obj <- FindNeighbors(seurat.obj, dims = 1:30)
for (res in c(0.05,0.1,0.3,0.5,0.8,1,1.2,1.4,1.5,2)){
  print(res)
  seurat.obj <- FindClusters(seurat.obj,resolution = res, algorithm = 1)
}

#umap可视化
cluster_umap <- wrap_plots(ncol = 5,
                           DimPlot(seurat.obj, reduction = "umap", group.by = "originalexp_snn_res.0.05", label = T) & NoAxes(),  
                           DimPlot(seurat.obj, reduction = "umap", group.by = "originalexp_snn_res.0.1", label = T) & NoAxes(),
                           DimPlot(seurat.obj, reduction = "umap", group.by = "originalexp_snn_res.0.3", label = T)& NoAxes(),
                           DimPlot(seurat.obj, reduction = "umap", group.by = "originalexp_snn_res.0.5", label = T) & NoAxes(),
                           DimPlot(seurat.obj, reduction = "umap", group.by = "originalexp_snn_res.0.8", label = T) & NoAxes(), 
                           DimPlot(seurat.obj, reduction = "umap", group.by = "originalexp_snn_res.1", label = T) & NoAxes(),
                           DimPlot(seurat.obj, reduction = "umap", group.by = "originalexp_snn_res.1.2", label = T) & NoAxes(),
                           DimPlot(seurat.obj, reduction = "umap", group.by = "originalexp_snn_res.1.4", label = T)& NoAxes(),
                           DimPlot(seurat.obj, reduction = "umap", group.by = "originalexp_snn_res.1.5", label = T)& NoAxes(),
                           DimPlot(seurat.obj, reduction = "umap", group.by = "originalexp_snn_res.2", label = T)& NoAxes()
)
cluster_umap
# ggsave(cluster_umap,filename = "Outplot/Step3.Before_inter.cluster_umap.pdf",
#        width = 25, height = 9)

# 选择一个合适的分辨率
Idents(object = seurat.obj) <- "originalexp_snn_res.0.3"
options(repr.plot.width = 6, repr.plot.height = 5)
p1 <- DimPlot(seurat.obj, reduction = "umap", group.by = "originalexp_snn_res.0.3", label = T)& NoAxes()

# Step4.亚群注释 --------------------------------------------------------------
options(repr.plot.width = 7.5, repr.plot.height = 10)
check_genes = c("S100A8","S100A9",'FCN1','CD52', # Mono-S100A8 
                'SOD2',"CXCL2",'C5AR1', # Mono-IL1B 
                'IFI30',"NUPR1","GPNMB",'CTSD', # Macro-NUPR1 
                "ACP5","HLA-DQA1",'PLD3',# Macro-CXCL9 
                "CCL2","HMOX1","IFI6","TNF",# Macro-CCL2 
                "IFI27","CCL3L1",'FABP5', # Macro-CCL3L1 
                "C3","FCGR2A", "ALOX5AP",'PLD4',# Macro-FOLR2 
                "FTH1","FTL","B2M","HSP90AA1", "TMSB4X", # Macro-FTH1 
                "SELENOP","LGMN","CD5L", "MUCL1","RSAD2", # Macro-SELENOP 
                "G0S2", "IL1B","RPS29", # Macro-IL1B 
                "FN1","CCL7","MT2A",'TGFBI','HSP90AB1', # Macro-FN1 
                "RNASE1","KLF6","CXCL3",# Macro-RNASE1 
                "CXCL10","IRF1","ISG15", # Macro-ISG15 
                "RAB42","PLIN2","MT1X","CFD", # Macro-MT1X 
                "TPSAB1","TPSB2",'CPA3', # MAST-TPSAB1
                "CPVL","CST3","LYZ","SNX3", # cDC1-CLEC9A 
                "RGS2","CXCR4","INSIG1","GPR183", # cDC2-CLEC10A
                "MARCKSL1","FSCN1","TXN","CCR7", # mDC-LAMP3
                "S100B","MMP9","LMNA","CD1A", # cDC2-CD1A
                "GZMB","JCHAIN","IRF7","ITM2C" # pDC-LILRA4
)

## 运行以下命令检查 check_genes 中是否有重复的基因名称
duplicated_genes <- check_genes[duplicated(check_genes)]
duplicated_genes

p1 <- DotPlot(object = seurat.obj, features = check_genes,assay = "originalexp",scale = T) + 
  coord_flip()