1. #单细胞亚群的生物学合理命名讨论
    2. rm(list = ls())
    3. library(Seurat)
    4. library(SeuratData)
    5. library(ggplot2)
    6. library(patchwork)
    7. library(dplyr)
    8. load(file = 'basic.sce.pbmc.Rdata')
    9. levels(Idents(pbmc))
    11. # 首先提取T细胞子集
    12. DimPlot(pbmc, 
    13.         reduction = 'umap', 
    14.         label = TRUE, pt.size = 0.5) + NoLegend()
    15. FeaturePlot(pbmc, features = c("CD3D","CD3E")) #T细胞常用CD3来分
    16. sce=pbmc
    17. table(Idents(sce))
    18. #取子集
    19. t_sce = sce[, Idents(sce) %in% c('CD8 T')]
    21. #然后进行标准的降维聚类分群
    22. #代码不要变动
    23. #降维聚类分群
    24. sce=t_sce
    26. #归一化
    27. sce <- NormalizeData(sce, 
    28.                      normalization.method = "LogNormalize",
    29.                      scale.factor = 1e4) 
    30. #Find variable features
    31. sce <- FindVariableFeatures(sce, 
    32.                             selection.method = 'vst',
    33.                             nfeatures = 2000)
    34. #Scale and center the data
    35. sce <- ScaleData(sce, 
    36.                  vars.to.regress = "percent.mt")
    37. #pca
    38. sce <- RunPCA(sce, 
    39.               features = VariableFeatures(object = sce)) 
    40. #计算给定数据集的k.param近邻
    41. sce <- FindNeighbors(sce, 
    42.                      dims = 1:10)
    43. #聚类
    44. sce <- FindClusters(sce, 
    45.                     resolution = 1 )
    47. # Look at cluster IDs of the first 5 cells
    48. head(Idents(sce), 5)
    49. table(sce$seurat_clusters) 
    50. sce <- RunUMAP(sce, dims = 1:10)
    51. DimPlot(sce,
    52.         reduction = 'umap')
    54. #为数据集中的每个标识类查找markers(差异表达基因)
    55. sce.markers <- FindAllMarkers(object = sce,
    56.                               only.pos = TRUE, 
    57.                               min.pct = 0.25, 
    58.                               thresh.use = 0.25)
    59. DT::datatable(sce.markers)
    60. # write.csv(sce.markers,file=paste0(pro,'_sce.markers.csv'))
    61. save(sce.markers,
    62.      file = 'sce.markers-for-cd8-subsets.Rdata')
    64. library(dplyr) 
    65. #在seurat4里面,是 avg_log2FC,是V3版本是avg_logFC
    66. top10 <- sce.markers %>% group_by(cluster) %>% top_n(5, avg_log2FC)
    67. DoHeatmap(sce,
    68.           top10$gene,
    69.           size=3)  
    70. p <- DotPlot(sce, 
    71.              features = unique(top10$gene),
    72.              assay='RNA' )  + coord_flip()
    74. p
    75. ggsave('check-top5-for-cd8-subsets.pdf' )
    78. # cytotoxicity (GZMB, PRF1),  0 
    79. # naive (LEF1, SELL, TCF7), 1 
    80. # LTB T cells , 2
    81. # RPl/s 
    83. # 参考： https://mp.weixin.qq.com/s/1YjoX8lTIB0e2vV0Eqey2Q 
    85. marker_genes= c("MKI67","TOP2A",'TNFRSF9','MX1',
    86.                 "SELL","IL7R","CD40LG","ANXA1","FOS")
    89. p <- DotPlot(sce, 
    90.              features = marker_genes,
    91.              assay='RNA'  )  + coord_flip() + theme(axis.text.x = element_text(angle = 90))
    93. p
    94. ggsave('check_g1_markers_by_Tcell-SubType.pdf')
    96. marker_genes= c("MKI67","TOP2A",'TNFRSF9','MX1')
    98. FeaturePlot(pbmc, 
    99.             features = marker_genes)