单细胞转录组之Seurat包全流程-数据过滤、降维分群及可视化

if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("Seurat")

library(Seurat)
# devtools::install_github('satijalab/seurat-data')
library(SeuratData)
library(ggplot2)
library(patchwork)
library(dplyr)
##如果有未安装的包，运行如下命令安装：install.packages("包的名字")

pbmc.data <- Read10X(data.dir = "sample/outs/filtered_feature_bc_matrix")
##使用未经标准化的数据创建Seurat对象
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
##project：项目名；过滤条件：min.cells = 3 一个基因最少在3个细胞中检测到，min.features = 200一个细胞中最少检测到200个基因

str(pbmc)

Formal class 'Seurat' [package "SeuratObject"] with 13 slots
  ..@ assays      :List of 1
  .. ..$ RNA:Formal class 'Assay' [package "SeuratObject"] with 8 slots
  .. .. .. ..@ counts       :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. .. .. .. .. ..@ i       : int [1:2282976] 29 73 80 148 163 184 186 227 229 230 ...
  .. .. .. .. .. ..@ p       : int [1:2701] 0 779 2131 3260 4220 4741 5522 6304 7094 7626 ...
  .. .. .. .. .. ..@ Dim     : int [1:2] 13714 2700
  .. .. .. .. .. ..@ Dimnames:List of 2
  .. .. .. .. .. .. ..$ : chr [1:13714] "AL627309.1" "AP006222.2" "RP11-206L10.2" "RP11-206L10.9" ...
  .. .. .. .. .. .. ..$ : chr [1:2700] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
  .. .. .. .. .. ..@ x       : num [1:2282976] 1 1 2 1 1 1 1 41 1 1 ...
  .. .. .. .. .. ..@ factors : list()
  .. .. .. ..@ data         :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  .. .. .. .. .. ..@ i       : int [1:2282976] 29 73 80 148 163 184 186 227 229 230 ...
  .. .. .. .. .. ..@ p       : int [1:2701] 0 779 2131 3260 4220 4741 5522 6304 7094 7626 ...
  .. .. .. .. .. ..@ Dim     : int [1:2] 13714 2700
  .. .. .. .. .. ..@ Dimnames:List of 2
  .. .. .. .. .. .. ..$ : chr [1:13714] "AL627309.1" "AP006222.2" "RP11-206L10.2" "RP11-206L10.9" ...
  .. .. .. .. .. .. ..$ : chr [1:2700] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
  .. .. .. .. .. ..@ x       : num [1:2282976] 1 1 2 1 1 1 1 41 1 1 ...
  .. .. .. .. .. ..@ factors : list()
  .. .. .. ..@ scale.data   : num[0 , 0 ] 
  .. .. .. ..@ key          : chr "rna_"
  .. .. .. ..@ assay.orig   : NULL
  .. .. .. ..@ var.features : logi(0) 
  .. .. .. ..@ meta.features:'data.frame':    13714 obs. of  0 variables
  .. .. .. ..@ misc         : list()
  ..@ meta.data   :'data.frame':    2700 obs. of  4 variables:
  .. ..$ orig.ident  : Factor w/ 1 level "pbmc3k": 1 1 1 1 1 1 1 1 1 1 ...
  .. ..$ nCount_RNA  : num [1:2700] 2419 4903 3147 2639 980 ...
  .. ..$ nFeature_RNA: int [1:2700] 779 1352 1129 960 521 781 782 790 532 550 ...
  .. ..$ percent.mt  : num [1:2700] 3.02 3.79 0.89 1.74 1.22 ...
  ..@ active.assay: chr "RNA"
  ..@ active.ident: Factor w/ 1 level "pbmc3k": 1 1 1 1 1 1 1 1 1 1 ...
  .. ..- attr(*, "names")= chr [1:2700] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ...
  ..@ graphs      : list()
  ..@ neighbors   : list()
  ..@ reductions  : list()
  ..@ images      : list()
  ..@ project.name: chr "pbmc3k"
  ..@ misc        : list()
  ..@ version     :Classes 'package_version', 'numeric_version'  hidden list of 1
  .. ..$ : int [1:3] 4 1 0
  ..@ commands    : list()
  ..@ tools       : list()

pbmc@assays

$RNA
Assay data with 13714 features for 2700 cells
First 10 features:
 AL627309.1, AP006222.2, RP11-206L10.2, RP11-206L10.9, LINC00115,
NOC2L, KLHL17, PLEKHN1, RP11-54O7.17, HES4 
##pbmc@assays$RNA中包含着2700个细胞的13714个基因的信息，包括原始counts(pbmc@assays$RNA@counts)

head(pbmc@meta.data)
                 orig.ident nCount_RNA nFeature_RNA percent.mt
AAACATACAACCAC-1     pbmc3k       2419          779  3.0177759
AAACATTGAGCTAC-1     pbmc3k       4903         1352  3.7935958
AAACATTGATCAGC-1     pbmc3k       3147         1129  0.8897363
AAACCGTGCTTCCG-1     pbmc3k       2639          960  1.7430845
AAACCGTGTATGCG-1     pbmc3k        980          521  1.2244898
AAACGCACTGGTAC-1     pbmc3k       2163          781  1.6643551

#pbmc[[]] 可进入meta.data 相当于pbmc@meta.data，，对meta.data下级数据调用方式跟数据框取列类似，pbmc@meta.data$percent.mt，另外pbmc[["percent.mt"]]也可以表示选择percent.mt这一列

#计算线粒体基因占总基因数目的百分比，pattern后可以是正则表达式。
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

#小提琴图展示meta.data中的数据，ncol 表示显示多个图时的列数
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt") 
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA") 
plot1 + plot2

pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 1e4) 

#normalization.method有LogNormalize，CLR，RC三种。
#LogNormalize:每个细胞的基因数数除以该细胞的总基因数，再乘以scale.factor。然后使用log（x+1）进行自然对数转换。
#CLR：应用居中的对数比率变换
#RC: 相对计数。每个细胞的功能计数除以该细胞的总计数，再乘以scale.factor。没有应用log转换。对于每百万计数(CPM)设置规模。系数= 1 e6

#选取2000个差异基因
pbmc <- FindVariableFeatures(pbmc, selection.method = 'vst', nfeatures = 2000)

#selection.method：vst，mean.var.plot (mvp)，dispersion (disp)。
#nfeatures:选择差异基因的数量;仅在selection.method设置为“dispersion”或“vst”时使用。

# 选择前10的差异基因
top10 <- head(VariableFeatures(pbmc), 10)
# 带标签与否的展示差异基因
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2

#消除线粒体基因的影响
pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")

pbmc <- RunPCA(pbmc, features = VariableFeatures( pbmc))  

#dim输入降维的维度，resolution分辨率，判断近邻距离大小，值越低，聚类越少。

pbmc <- FindNeighbors(pbmc, dims = 1:10)

?FindNeighbors
Description:
    Constructs a Shared Nearest Neighbor (SNN) Graph for a given
     dataset. We first determine the k-nearest neighbors of each cell.
     We use this knn graph to construct the SNN graph by calculating
     the neighborhood overlap (Jaccard index) between every cell and
     its k.param nearest neighbors.
这个地方说明，这个函数首先是计算每个细胞的KNN，也就是计算每个细胞之间的相互距离，依据细胞之间距离的graph来构建snn graph（依据细胞之间“邻居”的overlop）
这里有三个问题：1、knn是什么，2、Jaccard index又是什么  3、邻居的判定
我们来看看参数（主要参数）：
distance.matrix: Boolean value of whether the provided matrix is a
          distance matrix; note, for objects of class ‘dist’, this
          parameter will be set automatically
      这个参数我们通常不会设置，但是默认是TRUE。
k.param: Defines k for the k-nearest neighbor algorithm
      这个参数就是用来定义最相近的几个细胞作为邻居，默认是20
compute.SNN: also compute the shared nearest neighbor graph
     计算共享邻居的数量，一般不设置
prune.SNN: Sets the cutoff for acceptable Jaccard index when computing
          the neighborhood overlap for the SNN construction. Any edges
          with values less than or equal to this will be set to 0 and
          removed from the SNN graph. Essentially sets the strigency of
          pruning (0 - no pruning, 1 - prune everything).
     在计算SNN构造的邻域重叠时，为可接受的Jaccard index设置截止值。 值小于或等于此值的任何边将被设置为0并从SNN图中删除。 本质上设置修剪的严格程度（0-不修剪，1-修剪所有内容）。
nn.method: Method for nearest neighbor finding. Options include: rann,
          annoy
      这个参数提供了如何判断邻居的方法，提供的可选是rann和annoy，**这个我们在后面讨论**。
      annoy.metric: Distance metric for annoy. Options include: euclidean,
          cosine, manhattan, and hamming
      （annoy距离的方式）
force.recalc: Force recalculation of SNN
      SNN强制重新计算，一般不设置

pbmc <- FindClusters(pbmc, resolution = 0.5)

?FindClusters
Description:

     Identify clusters of cells by a shared nearest neighbor (SNN)
     modularity optimization based clustering algorithm. First
     calculate k-nearest neighbors and construct the SNN graph. Then
     optimize the modularity function to determine clusters. For a full
     description of the algorithms, see Waltman and van Eck (2013) _The
     European Physical Journal B_. Thanks to Nigel Delaney
     (evolvedmicrobe@github) for the rewrite of the Java modularity
     optimizer code in Rcpp!
这个说明，依据SNN来识别类群，当然算法很复杂，我们可以参考给的网址。
我们来看看主要参数
resolution: Value of the resolution parameter, use a value above
          (below) 1.0 if you want to obtain a larger (smaller) number
          of communities.
这个参数可以理解为清晰度，值越低，可以容纳更少的共享邻居数量，聚类数也会变少。
modularity.fxn: 计算模块系数函数，1为标准函数；2为备选函数，这里没有具体说明是什么函数，我认为1是上面提到的Kronecker delta函数。
method: Method for running leiden (defaults to matrix which is fast
          for small datasets). Enable method = "igraph" to avoid
          casting large data to a dense matrix.
      这个参数表示leiden算法的计算方式，（我对算法是小白~，求大神告知）
algorithm: 模块系数优化算法，1使用原始Louvain算法；2使用Louvain algorithm with multilevel refinement；3使用SLM算法；4使用Leiden算法（注：4需要额外安装插件）
n.start: 随机开始的数量，默认是10
random.seed: 随机数种子，默认是0

head(Idents(pbmc), 5)

AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 
               0                3                2 
AAACCGTGCTTCCG-1 AAACCGTGTATGCG-1 
               1                6 
Levels: 0 1 2 3 4 5 6 7 8

table(pbmc$seurat_clusters) 

  0   1   2   3   4   5   6   7   8 
709 480 429 342 316 162 154  32  14 

pbmc <- RunUMAP(pbmc, dims = 1:10)

#note that you can set `label = TRUE` or use the LabelClusters function to help label individual clusters

DimPlot(pbmc, reduction = 'umap',label = TRUE)

#FindMarkers查找亚群的标记基因
# 发现聚类1的所有biomarkers
cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)
head(cluster1.markers, n = 5)

# 查找将聚类1与聚类2和3区分的所有标记基因
cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(2, 3), min.pct = 0.25)
head(cluster5.markers, n = 5)

# 与所有其他亚群相比，找到每个亚群的标记，仅报告阳性细胞
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
top10markers<-pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)

cluster1.markers <- FindMarkers(pbmc, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)

DoHeatmap(pbmc , features = top10markers$gene, size = 3)

FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", 
"LYZ", "PPBP", "CD8A"))

#新id与之前id(0-8)一一对应
new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T","B", "CD8 T", "FCGR3A+ Mono", "NK", "DC", "Platelet")

names(new.cluster.ids) <- levels(pbmc)
#更改seuratd对象中的Idents
pbmc <- RenameIdents(pbmc, new.cluster.ids)

DimPlot(pbmc, reduction = 'umap', 
        label = TRUE, pt.size = 0.5) + NoLegend()

##以下features可换成任意基因
##小提琴图
VlnPlot(pbmc, features = c("MS4A1", "CD79A")) 

##坐标映射图
FeaturePlot(pbmc, features = c("MS4A1", "CD79A"))

##峰峦图
RidgePlot(pbmc, features = c("MS4A1", "CD79A"), ncol = 1)


features= c('IL7R', 'CCR7','CD14', 'LYZ',  'IL7R', 'S100A4',"MS4A1", "CD8A",'FOXP3',
            'FCGR3A', 'MS4A7', 'GNLY', 'NKG7',
            'FCER1A', 'CST3','PPBP')

##气泡图
DotPlot(pbmc, features = unique(features)) + RotatedAxis()

#downsample 在每个细胞亚群中均抽样100个细胞做热图
DoHeatmap(subset(pbmc, downsample = 100), 
          features = features, size = 3)