单细胞.png

    1. rm(list=ls())
    3. library(dplyr)
    4. library(Seurat)
    5. library(patchwork)
    7. ## =============1.Load the PBMC dataset
    8. data_dir <- "data/pbmc3k_filtered_gene_bc_matrices/"
    9. pbmc.data <- Read10X(data.dir = data_dir)
    10. pbmc.data
    12. # Initialize the Seurat object with the raw (non-normalized data)
    13. # min.cell：每个feature至少在多少个细胞中表达
    14. # min.features：每个细胞中至少有多少个feature被检测到
    15. pbmc <- CreateSeuratObject(counts = pbmc.data, 
    16.                            project = "pbmc3k", 
    17.                            min.cells = 3,
    18.                            min.features = 200)
    19. pbmc
    23. ## =============2.探索一下数据
    24. # count matrix长什么样？
    25. # 首先看看三个基因的count值，'.'values in the matrix represent 0s (no molecules detected)
    26. pbmc.data[1:30, 1:30]
    27. pbmc.data[c("CD3D", "TCL1A", "MS4A1"), 1:30]
    28. rownames(pbmc.data)
    29. colnames(pbmc.data)
    31. # 0值矩阵大小
    32. raw_counts <- as.matrix(pbmc.data)
    33. raw_counts[1:4,1:4]
    34. dense.size <- object.size(raw_counts)
    35. dense.size
    37. # 稀疏矩阵大小
    38. sparse.size <- object.size(pbmc.data)
    39. sparse.size
    41. # 压缩比
    42. dense.size/sparse.size
    46. ## =============3.QC: 
    47. # 1.每个细胞中RNA的个数；
    48. # 2，每个细胞中线粒体基因表达的比例
    49. # The [[ operator can add columns to object metadata. This is a great place to stash QC stats
    50. grep(pattern = "^MT-", rownames(pbmc), value = T)
    51. pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
    54. # QC指标
    55. head(pbmc@meta.data, 5)
    58. # QC指标使用小提琴图可视化
    59. VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), 
    60.         ncol = 3,pt.size = 0.1)
    63. # 两个指标之间的相关性
    64. plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", 
    65.                         feature2 = "percent.mt")
    66. plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", 
    67.                         feature2 = "nFeature_RNA")
    68. plot1 + plot2
    70. # 过滤
    71. pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
    72. pbmc
    76. ## =============4.标准化
    77. # LogNormalize
    78. pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", 
    79.                       scale.factor = 10000)
    81. # 标准化后的值保存在：pbmc[["RNA"]]@data
    82. normalized.data <- pbmc[["RNA"]]@data
    83. normalized.data[1:20,1:4]
    84. dim(normalized.data)
    88. ## =============5.鉴定高变基因
    89. # 高变基因：在一些细胞中表达高，另一些细胞中表达低的基因
    90. # 变异指标： mean-variance relationship
    91. # 默认返回两千个高变基因，用于下游如PCA降维分析。
    92. pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", 
    93.                              nfeatures = 2000)
    95. # 提取前10的高变基因
    96. top10 <- head(VariableFeatures(pbmc), 10)
    97. top10
    99. # 展示高变基因
    100. plot1 <- VariableFeaturePlot(pbmc)
    101. plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
    102. plot1
    103. plot2
    107. ## =============6.Scaling the data
    108. # 归一化处理：每一个基因在所有细胞中的均值变为0，方差标为1，对于降维来说是必需步骤
    109. # 归一化后的值保存在：pbmc[["RNA"]]@scale.data
    110. pbmc <- ScaleData(pbmc)
    111. scale.data <- pbmc[["RNA"]]@scale.data
    112. dim(scale.data)
    113. scale.data[1:10,1:4]
    116. # 后面绘制热图使用的是归一化后的值，可以选择全部基因归一化，但是耗时久一点
    117. all.genes <- rownames(pbmc)
    118. pbmc <- ScaleData(pbmc, features = all.genes)
    122. ## =============7.降维
    123. # PCA降维，默认使用前面2000个高变基因，可以使用features改变用于降维的基因集
    124. pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
    126. # 可视化
    127. VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
    128. DimPlot(pbmc, reduction = "pca")
    129. DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
    130. DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)
    134. ## =============8.确定使用PC个数
    135. # # each PC essentially representing a ‘metafeature’
    136. # 方式1：根据p值判断，比较耗时
    137. # pbmc <- JackStraw(pbmc, num.replicate = 100)
    138. # pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
    139. # 
    140. # JackStrawPlot(pbmc, dims = 1:15) 
    142. # 方式2：肘部法，根据拐点判断，比较主观，尽量多选几个PC
    143. ElbowPlot(pbmc, ndims = 30) 
    147. ## =============9.对细胞聚类
    148. # 选取10个PC，首先基于PCA空间构建一个基于欧氏距离的KNN图
    149. pbmc <- FindNeighbors(pbmc, dims = 1:10)
    151. # 聚类并最优化
    152. # resolution参数：值越大，细胞分群数越多，
    153. # 0.4-1.2 typically returns good results for single-cell datasets of around 3K cells
    154. # Optimal resolution often increases for larger datasets. 
    155. pbmc <- FindClusters(pbmc, resolution = 0.5)
    157. # 查看聚类数ID
    158. head(Idents(pbmc), 5)
    160. # 查看每个类别多少个细胞
    161. head(pbmc@meta.data)
    162. table(pbmc@meta.data$seurat_clusters)
    165. # resolution对聚类的影响
    166. res.used <- seq(0.1,1,by=0.2)
    167. res.used
    168. for(i in res.used){
    169.   pbmc <- FindClusters(pbmc, verbose = T, resolution = res.used)
    170. }
    172. # 可视化
    173. library(clustree)
    174. clus.tree.out <- clustree(pbmc) +
    175.                  theme(legend.position = "bottom") + 
    176.                  scale_color_brewer(palette = "Set1") +
    177.                  scale_edge_color_continuous(low = "grey80", high = "red")
    179. clus.tree.out
    181. pbmc <- FindClusters(pbmc, verbose = T, resolution = 0.5)
    182. ## =============10.将细胞在低维空间可视化UMAP/tSNE
    183. set.seed(123)
    184. pbmc <- RunUMAP(pbmc, dims = 1:10)
    185. pbmc <- RunTSNE(pbmc, dims = 1:10)
    187. # 可视化
    188. DimPlot(pbmc, reduction = "umap", label = T, label.size = 5)
    189. DimPlot(pbmc, reduction = "tsne", label = T, label.size = 5)
    190. saveRDS(pbmc, file = "data/pbmc_tutorial.rds")
    192. #pbmc <- readRDS("data/pbmc_tutorial.rds")
    195. ## =============11.差异表达分析
    196. # 在cluster2 vs else中差异表达
    197. cluster1.markers <- FindMarkers(pbmc, ident.1 = 2, min.pct = 0.25)
    198. head(cluster1.markers, n = 5)
    201. # 指定两个类cluster 5 from clusters 0 and 3
    202. cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
    203. head(cluster5.markers, n = 5)
    206. # 所有类的差异表达基因
    207. # only.pos：只保留上调差异表达的基因
    208. pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
    209. write.csv(pbmc.markers,file = "data/pbmc.markers.csv")
    210. head(pbmc.markers)
    213. #  筛选
    214. pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)
    217. # 选择一些基因进行可视化，作者这里根据自己的知识背景选择的相关基因
    218. VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
    219. VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)
    220. FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
    221. DotPlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A")) + coord_flip() + 
    222.    theme(axis.text.x = element_text(angle=90, hjust=1, vjust=.5))
    224. top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)
    225. DoHeatmap(pbmc, features = top10$gene) + NoLegend()
    230. ## =============12.细胞类型鉴定：使用经典的markers进行注释
    231. new.cluster.ids <- c("Naive CD4 T",   # IL7R, CCR7
    232.                      "CD14+ Mono",    # CD14, LYZ
    233.                      "Memory CD4 T",  # IL7R, S100A4
    234.                      "B",             # MS4A1
    235.                      "CD8 T",         # CD8A
    236.                      "FCGR3A+ Mono",  # FCGR3A, MS4A7
    237.                      "NK",            # GNLY, NKG7
    238.                      "DC",            # FCER1A, CST3
    239.                      "Platelet")      # PPBP
    241. names(new.cluster.ids) <- levels(pbmc)
    242. new.cluster.ids
    244. # 原来的细胞聚类名称
    245. Idents(pbmc)
    247. # 更改细胞聚类的名字
    248. pbmc <- RenameIdents(pbmc, new.cluster.ids)
    249. Idents(pbmc)
    252. # 可视化
    253. DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 2,label.box = T) + NoLegend()
    254. DimPlot(pbmc, reduction = "tsne", label = TRUE, pt.size = 2,label.box = T) + NoLegend()
    256. # 保存结果
    257. pbmc@meta.data$cell_anno <- Idents(pbmc)
    258. head(pbmc@meta.data)
    259. write.csv(pbmc@meta.data,file = "data/metadata.csv")
    260. saveRDS(pbmc, file = "data/pbmc3k_final.rds")
    261. table(pbmc@meta.data$cell_anno)
    266. ## 去掉画图部分的流程总结
    267. data_dir <- "data/pbmc3k_filtered_gene_bc_matrices/"
    268. pbmc.counts <- Read10X(data.dir = data_dir)
    269. pbmc <- CreateSeuratObject(counts = pbmc.counts)
    270. pbmc <- NormalizeData(object = pbmc)
    271. pbmc <- FindVariableFeatures(object = pbmc)
    272. pbmc <- ScaleData(object = pbmc)
    273. pbmc <- RunPCA(object = pbmc)
    274. pbmc <- FindNeighbors(object = pbmc)
    275. pbmc <- FindClusters(object = pbmc)
    276. pbmc <- RunTSNE(object = pbmc)
    277. pbmc <- RunUMAP(object = pbmc)
    278. DimPlot(object = pbmc, reduction = "tsne")

    代码及图片均来自于生信技能树张娟老师