一、数据前处理
#保持好习惯，先清空一下环境，设置保存路径
rm(list = ls())
setwd("")1.数据读取
载入包
#载入程辑包
#multtest包安装tips,
#install.packages("BioManager")
#BiocManager::install("multtest"),
##如果出现显示R版本不符的error不要着急换R版本，直接打开R软件选择就近的镜像安装。
library(multtest)
library(Seurat)
library(dplyr)
library(patchwork)
library(R.utils)导入数据
CreateSeuratObject用于创建Seurat对象
pbmc.data <- Read10X(data.dir = "hg19")
#Load the PBMC dataset,创建Seurat对象，counts为读取的源文件，project为Seurat对象想保存的文件名，
#可以加上限定条件：min.cells为组织中分离的最少细胞数，min.features为一个细胞中测出的最少的基因数量
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
#查看pbmc中细胞数量的3个方法
pbmc
## An object of class Seurat
## 13714 features across 2700 samples within 1 assay ,细胞数量是2700个
## Active assay: RNA (13714 features, 0 variable features)
ncol(pbmc)
## [1] 2700
ncol(pbmc.data)2.进行质控
我们提到单细胞数据质控的时候，⼀般是指细胞的过滤，其实是从⼀个barcode X gene矩阵中过滤掉⼀部分不是细胞的barcode，如细胞碎⽚，双细胞，死细胞等。这三类barcode的特征可以通过其对应的基因表达情况来描述：nCount(总基因表达数)、nFeature(总基因数)、
percent.HB(红细胞基因表达⽐例)、percent.MT(线粒体基因表达⽐例)。nCount和nFeature过⾼可能是双细胞，过低可能是细胞碎⽚。percent.HB刻画红细胞⽐例，percent.MT刻画细胞状态，值过⾼可能是濒临死亡的细胞
PercentageFeatureSet 函数是根据counts总数相除算的打分：该基因集的counts总和/所有基因的counts总和。
基因集占的百分比 = 分子 / 分母 * 100；
分子： 指定基因的 counts种总和，在这次分析中是指细胞的线粒体基因转录本数
分母： 从 meta.data 获取 nCount_RNA 列，就是每个cell中所有基因的 counts总和
这个函数的意思是，每个细胞的线粒体的基因数/细胞总基因数的百分比，将作为一列数据加到pbmc的metadata中，这一列的列名为“percent.mt”
#进行质控，计算每个细胞的线粒体基因转录本数的百分比（%），使用[[ ]] 操作符存放到metadata中；
#分析的时候要确认好物种，如果是小鼠的，要用mt
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
#在创建对象 CreateSeuratObject() 的过程中，会自动计算细胞中独特基因与总基因数目，可以在目标对象中找到
head(pbmc@meta.data,5)
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.22448983.质控结果的可视化
目的:通过可视化质控结果，进一步去掉低质量的数据
（1）小提琴图，将pbmc中的meta@data中的3个指标成图
#(1)将QC结果展示为小提琴图，features中的名称加引号，ncol=3表示图形分三列展示
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
（2）散点图
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
（3）去掉极值
，数据质控的原则是：You can’t stop it if you can’t see it. 即，如果不能确定是数据中的垃圾，就保留下来，因为你去掉的可能是⼀个重要的基因或重要因素或重要稀有细胞亚群，所以去掉极值时要慎重，类似于跑流式的时候画门
#删除极值
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
#删除极值后的细胞数
ncol(as.data.frame(pbmc[["RNA"]]@counts))
##26384.数据标准化
方法一 normalization函数，默认LogNormalize的方法
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)方法二
SCTransform，是一种三合一的方法，可以将质控，归一化和去识别高变基因合为一体
pbmc <- SCTransform(pbmc, vars.to.regress = "percent.mt", verbose = FALSE)5.特征性基因的选择
计算数据集中表现出高细胞间差异的基因子集(即，它们在一些细胞中高表达，而在另一些细胞中低表达)。在下游分析中关注这些基因有助于突出单细胞数据集中的生物信号。默认情况下，每个数据集选择2000个高变异基因用于下游分析。
FindVariableFeatures 函数有 3 种选择高表达变异基因的方法，可以通过 selection.method参数来选择，它们分别是： vst（默认值）， mean.var.plot 和 dispersion。 nfeatures 参数的默认值是 2000，可以改变。如果 selection.method 参数选择的是 mean.var.plot，就不需要人为规定高表达变异基因的数目，算法会自动选择合适的数目。
#选择高变基因，使用“vst”方法选择2000个高变基因
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
#储存前10位高变基因
top10 <- head(VariableFeatures(pbmc), 10)
plot3 <- VariableFeaturePlot(pbmc)#高变基因散点图
plot4 <- LabelPoints(plot=plot3,
points = top10,
repel=TRUE)#top10加上基因名标签
CombinePlots(plots = list(plot3, plot4), ncol =2)#结合到一张图中
二、线性降维-PCA分析
1.scale data
小Q：同样是对数据进行标准化处理，单细胞分析中的NormalizeData()与ScaleData()区别在哪儿？
NormalizeData()实际上就是在消除不同细胞测序深度的影响，测序深度简而言之就是测序得到的碱基数与基因组碱基数之间的比值。一般来说，对一个细胞而言，测序深度越大，每个基因检测到的reads会越多。这实际上就会造成一个问题，在处理单细胞数据时，如果一个基因在一个细胞当中的reads数很少，这究竟是这个基因本身就表达很少，还是因为这个细胞本身测序得到的reads数就很少呢？单细胞数据当中有很多基因的reads数很多，甚至上千，但是有很多基因却是个位数甚至0，那这种数据离散程度也是很大的，但是我们会发现当我们对1000取以10为底的对数时，就变成了3，对10取以10为底的对数时，就变成了1，这样就实现了降低数据离散程度的目的。但是这也存在一定的问题，如果一个基因的reads数为0，那岂不是不能取对数？我们对所有的值都加上1不就能够解决这个问题？
所以总结起来，NormalizeData()这个函数是首先对基因的reads数进行了同一文库大小的校正，然后再对校正后的数值进行对数化操作。
ScaleData()是对基因表达量的数值进行了z-score的转换，使其服从正态分布，同时也为后面的pca分析做了铺垫，因为pca分析默认数据是服从正态分布的。
#Scale data
#将均质为0标准差为1的数据记为标准化数据,用于下一步进行PCA分析
##改变每个基因的表达，使细胞间的平均表达为0；
##缩放每个基因的表达，使细胞间的差异为1（这一步在下游分析中给予相同的权重，使得高表达基因就不会影响后续分析）；
##缩放的数据储存在pbmc[[“RNA”]]@scale.data / pbmc@assays$RNA@scale.data中。
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc,features = all.genes)
#对缩放后的数据进行降维，默认降至50维。
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))2、降维结果可视化
目的是看哪个维度需要保留或排除，不要只把它们看成一堆字母。现在我们知道如果每个PC的基因大多是细胞类型的makergene，那么该PC就是细胞类型特异的，肯定是要保留的了。
#选择前5个维度进行查看
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
PC_ 1
Positive: CST3, TYROBP, LST1, AIF1, FTL
Negative: MALAT1, LTB, IL32, IL7R, CD2
PC_ 2
Positive: CD79A, MS4A1, TCL1A, HLA-DQA1, HLA-DQB1
Negative: NKG7, PRF1, CST7, GZMB, GZMA
PC_ 3
Positive: HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPB1
Negative: PPBP, PF4, SDPR, SPARC, GNG11
PC_ 4
Positive: HLA-DQA1, CD79B, CD79A, MS4A1, HLA-DQB1
Negative: VIM, IL7R, S100A6, IL32, S100A8
PC_ 5
Positive: GZMB, NKG7, S100A8, FGFBP2, GNLY
Negative: LTB, IL7R, CKB, VIM, MS4A7
#使用DimPlot函数，提取两个主成分(默认前两个，可以修改dims选项)绘制散点图。
#这个图中每个点就是一个样本，它根据PCA的结果坐标进行画图，这个坐标保存在了cell.embeddings中
DimPlot(pbmc, reduction = "pca")
head(pbmc[['pca']]@cell.embeddings)[1:2,1:2]
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
#使用DimHeatmap可视化函数查看样品中前500个细胞在前6个PCA中的热图：
DimHeatmap(pbmc, dims = 1:6, cells = 500, balanced = TRUE)
3.确定数据集的维度
根据前边的降维结果的可视化，结合JackStraw（）或ElbowPlot()函数判断后续分析应该保留多少个维度合适。
在JackStraw()函数中, 使用基于零分布的置换检验方法。随机抽取一部分基因(默认1%)然后进行pca分析得到pca分数，将这部分基因的pca分数与先前计算的pca分数进行比较得到显著性p-Value。根据主成分（pc）所包含基因的p-value进行判断选择主成分。最终的结果是每个基因与每个主成分的关联的p-Value。保留下来的主成分是那些富集小的p-Value基因的主成分。 该 JackStrawPlot() 函数提供了一个可视化工具，用于将每个 PC 的 p 值分布与均匀分布（虚线）进行比较。“ 显着 ”PC 将显示出大量具有低 p 值的特征（虚线上方的实线）。在这种情况下，在前 10-12 个 PC之后，重要性似乎急剧下降。
ElbowPlot在图中展示了每个主成分对数据方差的解释情况（百分比表示），并进行排序。根据自己需要选择主成分，图中发现第9个主成分是一个拐点，后续的主成分(PC)变化都不大了。
JackStraw()函数
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
JackStrawPlot(pbmc, dims = 1:15)
或ElbowPlot()函数
#ElbowPlot
ElbowPlot(pbmc, ndims = 20)
三、细胞聚类
1.聚类
#细胞聚类
#计算最邻近距离
pbmc <- FindNeighbors(pbmc, dims = 1:10)
#聚类，包含设置下游聚类的“间隔尺度”的分辨率参数resolution ，增加值会导致更多的聚类。
pbmc <- FindClusters(pbmc, resolution = 0.5)
#可以使用idents函数找到聚类情况：
#查看前5个细胞的聚类id
head(Idents(pbmc), 5)2.tsne和umap可视化
注意:umap的使用需要先在电脑上搭建python环境，建议下载3.6以上版本。在b站上选一个教程跟着做就好。
安装好后，直接在Rstudio中打开Terminal，输入：sudo /usr/local/bin/pip install --upgrade virtualenv
再运行以下代码，检查是否下载成功标志就是：library(umap)能否成功，这个是Seurat调用的基础。最后特别注意，要重启一下Rstudio再运行RunUMAP
reticulate::py_install(packages ='umap-learn')# install UMAP： reticulate::py_install(packages ='umap-learn')
pbmc <- RunUMAP(pbmc, dims = 1:10)
pbmc <- RunTSNE(pbmc, dims = 1:10)
#画图
DimPlot(pbmc, reduction = "umap")
DimPlot(pbmc, reduction = "tsne")
#添加细胞类群xiba的标签
DimPlot(pbmc, reduction = "umap",label = TRUE)
LabelClusters(DimPlot(pbmc, reduction = "umap"),id = 'ident')
#可以在这里保存为以下形式，就不必再执行以上步骤
saveRDS(pbmc, file = "../output/pbmc_tutorial.rds")
四、差异分析及可视化
1.差异分析
FindAllMarkers()这个函数会将cluster中的某一类和剩下的那些类比较，来找差异基因。与其对应的有个FindMarkers()函数，可以指定哪两个cluster对比。logfc.threshold表示logfc的阈值。
min.pct参数:设定在两个细胞群中任何一个被检测到的百分比，通过此设定不检测很少表达基因来缩短程序运行时间，默认0.1
#差异分析
#寻找cluster1的差异基因
cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)
# 找到cluster5和cluster0、3之间的markgene)
cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
#找到每个cluster相比于其余cluster的markgene，只报道阳性的markgene
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
## 所有基因先分组，再根据avg_log2FC进行排序
pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)
2.添加细胞注释信息
关于怎么样通过marker确定细胞类型，参考网址：秀儿！10+生信分析最大的难点在这里！30多种方法怎么选？今天帮你解决！_pbmc (sohu.com)
#添加细胞注释信息,如果知道每一个cluster代表的细胞类群的话
new.cluster.ids <- c("Naive CD4 T", "Memory CD4 T", "CD14+ Mono", "B", "CD8 T", "FCGR3A+ Mono",
"NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
3.marker gene可视化
（1）vlnplot
#VlnPlot: 基于细胞类群的基因表达概率分布
VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)
(2)查看特定基因在聚类图中的表达量分布情况
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP",
"CD8A"))
(3) 查看特定基因在cluster中的表达山脊图
features <- c("LYZ", "CCL5")
RidgePlot(pbmc, features = features, ncol = 2)
(4)查看特定基因在cluster中的分布点图
#查看特定基因在cluster中的分布点图，点的大小代表表达该基因的细胞比例，颜色代表平均表达水平
DotPlot(pbmc, features = features) + RotatedAxis()
(5)DoHeatmap为指定的细胞和基因画表达热图
#DoHeatmap为指定的细胞和基因画表达热图。每个类群默认展示top 20标记基因
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()
文章出处登录后可见！
