1. #复习seurat流程
    2. rm(list = ls())
    3. library(Seurat)
    4. library(SeuratData)
    5. library(ggplot2)
    6. library(patchwork)
    7. library(dplyr)
    9. # Load the PBMC dataset
    10. # https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz
    11. pbmc.data <- Read10X(data.dir = "filtered_gene_bc_matrices/hg19/")
    12. # Initialize the Seurat object with the raw (non-normalized data).
    14. pbmc <- CreateSeuratObject(counts = pbmc.data, 
    15.                            project = "pbmc3k",
    16.                            min.cells = 3,
    17.                            min.features = 200)
    18. pbmc
    20. pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, 
    21.                                              pattern = "^MT-") #计算属于给定特性集的所有计数的百分比
    22. #将QC度量可视化为小提琴图
    23. VlnPlot(pbmc, 
    24.         features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), 
    25.         ncol = 3)
    26. #FeatureScatter通常用于可视化特征-特征关系，但也可以用于对象计算的任何内容，例如对象元数据的列，PC分数等。
    27. plot1 <- FeatureScatter(pbmc, 
    28.                         feature1 = "nCount_RNA", 
    29.                         feature2 = "percent.mt") #单细胞数据的散点图
    30. plot2 <- FeatureScatter(pbmc, 
    31.                         feature1 = "nCount_RNA", 
    32.                         feature2 = "nFeature_RNA") #单细胞数据的散点图
    33. plot1 + plot2
    35. pbmc <- subset(pbmc,
    36.                subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
    37. pbmc <- NormalizeData(pbmc,
    38.                       normalization.method = "LogNormalize",
    39.                       scale.factor = 1e4) 
    40. pbmc <- FindVariableFeatures(pbmc, 
    41.                              selection.method = 'vst', 
    42.                              nfeatures = 2000)
    43. #找出10个最易突变的基因
    44. top10 <- head(VariableFeatures(pbmc), 10)
    45. #绘制带有或不带有标签的变量特征
    46. plot1 <- VariableFeaturePlot(pbmc)
    47. plot2 <- LabelPoints(plot = plot1, 
    48.                      points = top10, 
    49.                      repel = TRUE)
    50. plot1 + plot2
    52. pbmc <- ScaleData(pbmc, 
    53.                   vars.to.regress = "percent.mt")
    54. pbmc <- RunPCA(pbmc, 
    55.                features = VariableFeatures(object = pbmc))  
    57. pbmc <- FindNeighbors(pbmc, dims = 1:10)
    58. pbmc <- FindClusters(pbmc, resolution = 0.5)
    60. #查看前5个单元格的集群id
    61. head(Idents(pbmc), 5)
    62. table(pbmc$seurat_clusters) 
    63. pbmc <- RunUMAP(pbmc, dims = 1:10)
    65. #可以设置' label = TRUE '或使用LabelClusters功能来帮助标记单个集群
    66. DimPlot(pbmc, reduction = 'umap')
    68. VlnPlot(pbmc, 
    69.         features = c("MS4A1", "CD79A")) 
    70. FeaturePlot(pbmc, 
    71.             features = c("MS4A1", "GNLY", "CD3E", 
    72.                                "CD14", "FCER1A", "FCGR3A", 
    73.                                "LYZ", "PPBP", "CD8A"))
    75. new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T",
    76.                      "B", "CD8 T", "FCGR3A+ Mono", "NK", "DC", "Platelet")
    77. names(new.cluster.ids) <- levels(pbmc)
    78. pbmc <- RenameIdents(pbmc, new.cluster.ids)
    79. DimPlot(pbmc, 
    80.         reduction = 'umap', 
    81.         label = TRUE, 
    82.         pt.size = 0.5) + NoLegend()
    84. VlnPlot(pbmc, 
    85.         features = c("MS4A1", "CD79A")) 
    86. FeaturePlot(pbmc, 
    87.             features = c("MS4A1", "CD79A"))
    88. RidgePlot(pbmc,
    89.           features = c("MS4A1", "CD79A"), 
    90.           ncol = 1)
    91. features= c('IL7R', 'CCR7','CD14', 'LYZ',  'IL7R', 'S100A4',"MS4A1", "CD8A",'FOXP3',
    92.             'FCGR3A', 'MS4A7', 'GNLY', 'NKG7',
    93.             'FCER1A', 'CST3','PPBP')
    94. DotPlot(pbmc, 
    95.         features = unique(features)) + RotatedAxis()
    96. DoHeatmap(subset(pbmc ), 
    97.           features = features,
    98.           size = 3)
    99. DoHeatmap(subset(pbmc, downsample = 100), 
    100.           features = features, 
    101.           size = 3)
    103. save(pbmc,file = 'basic.sce.pbmc.Rdata')

    image.png