单细胞测序分析

1. R.studio先安装再调用相应的包（可参考上一文，这里就直接调用了）

library(tidyverse)
library(Seurat)
library(devtools)

2. 在R.studio的Terminal模块下载数据（注：提前设置好下载目录）：

wget http://cf.10xgenomics.com/samples/cell-exp/3.0.2/5k_pbmc_v3/5k_pbmc_v3_filtered_feature_bc_matrix.tar.gz

3.  利用Read10X函数读取数据集（其中包含3个文件：barcodes细胞/features基因/matrix基因表达矩阵）

pbmc.data <- Read10X(data.dir = "filtered_feature_bc_matrix/")

4. 创建Seurat对象（counts即读取的数据集，project自定义命名为pbmc5k，min.cells即筛选一个feature至少在3个细胞中存在，min.features即筛选细胞至少包含200个features，也就是基因个数）——展示数据 ——展示pbmc的前几行

pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc5k", min.cells = 3, min.features = 200)
pbmc
head(pbmc@meta.data)

6.  查找携带线粒体DNA（MT-）的标记（由于线粒体DNA含量过高一般认为是细胞活性差即将死亡的细胞，但是也要具体细胞具体分析） ——计算线粒体DNA所占百分比，将线粒体DNA含量百分比增加为pbmc数据的一列（后续要作为过滤条件） ——显示pbmc数据集的meta.data前几行 ——通过小提琴图来展示含量的差异

grep(pattern = "^MT-", rownames(pbmc), value = T) 
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
pbmc@meta.data %>% head()
vlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

7. 通过设定阈值来过滤数据（比较主观）——显示数据的变化

pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nCount_RNA < 5000 & percent.mt < 25)
pbmc

8. 标准化（全局缩放的归一化方法，目前seurat目前采用SCTransform）—— 筛选出差异表达的基因

pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

9. 选出前十个显著性差异基因并在图中进行基因名称标记。

top10 <- head(VariableFeatures(pbmc), 10)
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
CombinePlots(plots = list(plot1, plot2), ncol =1)

10. 进行数据缩放标准化ScaleData()（目的：将高表达基因和低表达基因赋予一致的权重，此步骤是后续PCA分析所必需的预处理）

pbmc <- ScaleData(pbmc)
scale.genes <-  VariableFeatures(pbmc)
pbmc <- ScaleData(pbmc, features = scale.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(pbmc)) 
#区分scale()和ScaleData()适用情况：下游分析中的PCA分析/umap、tsne聚类均用高变基因scale.data进行后续分析；基因可视化分析中，FeaturePlot、FeatureScatter、VlnPlot、DotPlot等函数默认slot =“data”，只有DoHeatmap()默认使用slot = “scale.data”，多个基因跨细胞比较；FindAllMarkers()找差异基因是默认slot =“data”，它是针对所有基因找差异基因，而不是高可变基因集

all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc[["RNA"]]@counts[1:6, 1:6]
pbmc[["RNA"]]@data[1:6, 1:6]
pbmc[["RNA"]]@scale.data[1:6, 1:6]

pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")

11. PCA分析“线性降维”并展示图

pbmc <- RunPCA(pbmc, features = VariableFeatures(pbmc), verbose = FALSE)
DimPlot(pbmc, reduction = "pca")

12. PCA分析后，确定PCs数。参考方法如下：（1）Heatmap （2）Elbow plot （3）JakStraw 

#Elbow plot 看每个PC对变异的贡献情况，通过elbow出现的地方x及后续曲线基本处于平缓，则认为真实信号主要来源于前x个PC
ElbowPlot(pbmc)

#JakStraw 函数是计算每个PC的各个基因的p值，用于判断哪些主要成分更具有统计学意义；并通过ScoreJackStraw（）量化主成分的显著性强度，富含低P值基因较多的主成分更有统计学意义。
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:50)
JackStrawPlot(pbmc, dims = 1:30)

#进一步量化确定每个PC对变异的贡献大小
mat <- pbmc[["RNA"]]@scale.data 
pca <- pbmc[["pca"]]
# Get the total variance:
total_variance <- sum(matrixStats::rowVars(mat))
eigValues = (pca@stdev)^2  ## EigenValues
varExplained = eigValues / total_variance
varExplained %>% enframe(name = "PC", value = "varExplained" ) %>%
ggplot(aes(x = PC, y = varExplained)) + 
geom_bar(stat = "identity") +
theme_classic() +
ggtitle("scree plot")

13. 细胞分群：

#识别细胞类群（通过计算每个细胞之间的相互距离，依据细胞之间“邻居”的overlap。共享邻居越多，则两个细胞群越相近）
pbmc <- FindNeighbors(pbmc, dims = 1:20)
#FindClusters（）参数resolution是影响分群的关键因素，通常认为0.4～1.2可以获得较好结果。
pbmc <- FindClusters(pbmc, resolution = 0.5)
head(Idents(pbmc), 5)
#非线性降维——UMAP
pbmc <- RunUMAP(pbmc, dims = 1:20)
DimPlot(pbmc, reduction = "umap"， label = TRUE)
#非线性降维——tSNE
pbmc <- RunTSNE(pbmc, dims = 1:20)
DimPlot(pbmc, reduction = "tSNE")

14. 筛选标记基因并进行细胞鉴定（FindMarkers函数可寻找指定细胞群的标记基因。原理是计算和其他细胞群的表达差异倍率最大的基因，即为该细胞群的标记基因）

#寻找细胞群cluster1的标记基因，参数min.pct是指基因必须至少在25%细胞中都有表达，才能被考虑在内（默认是0.1）。
cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)

#区分细胞群5与细胞群0和3的标记基因并定义为cluster5.markers
cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)

#一次性寻找每个群的标记基因
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)

#group_by将基因按细胞群分类，然后列出差异表达前两个基因。即为标记基因。
pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)

15. 标记基因可视化

#小提琴图
VlnPlot(pbmc, features = c("MS4A1", "CD79A"))

#细胞分群图
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
p<- FeaturePlot(pbmc, features = "CD14")
p_after<- p
p_after$data <- p_after$data[order(p_after$data$CD14),]
CombinePlots(plots = list(p, p_after))

#Heatmap图示（对scale进行赋值，以防标注的基因名乱码）
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)
pbmc@assays$RNA@scale.data <- scale(pbmc@assays$RNA@data, scale = TRUE)
DoHeatmap(pbmc, features = top10$gene, scale = scale.data) + NoLegend()