盲人摸象--single cell sequencing的Seu
只有下水才会学会游泳，今天我想下水了。首先本人对单细胞测序数据的处理一无所知，但看过一些单细胞测序的文章，有一些R语言基础和跑过一些RNA-seq上游分析的经验，流程是跑过但debug是全不会。
框架
1. 单细胞样本的准备：
无论是哪种生信数据分析都非常依赖于良好的测序数据。这里的良好不仅仅是建库和测序整个过程的良好质控。事实上，测序业务已经愈发成熟，主要建库合格，数据的QC基本上都不会差。那么作为一个常年做湿实验的人，太清楚实验样本的可变性有多大了。即使是同一株细胞，百分之30融合率的状态和百分之80融合度的状态可能不一样，为了尽可能得到贴近事实的测序数据，则实验的把控则相当重要。细胞培养尚且算是容易的，动物样本的采集过程中可能会引物非常多的变量。例如：动物死亡后是否能及时取样，不同动物取样时的手法力度是否一样，取样后是否及时进行了下一步操作，操作过程中有没有无偏差，有没有混用手术器械有没有编错管子等等非常细节的东西都会影响实验的结果。尤其是大多数样本都不是同一天取的，甚至不是同时测序的，这不免还会有批次效应，在这里我也说不清楚，毕竟有些文章注意到了，而有些文章没有做，很有可能是我的浅显不能理解。
2. 单细胞测序的种类；
这里没有特意去查，在我的脑海中有：(1) smart-seq，通量小但测全长，是将一颗细胞尽可能的测完；（2）10x之类的加barcode，通量大只测带3'的，是尽可能将所有细胞都测到，因此对于单一细胞来说可能深度不够。（以下这张图是从豆豆的简书拿来的）
image.png
3. 单细胞测序数据下载；
4. 单细胞数据的格式或者类型；
5. 单细胞数据的大概流程（欢迎拍砖）：
(1) no-bias clustering;
(2) define each cluster according to top genes and GO enrichment results, calculate the proportion of all clusters;
(3) imply gene-sets for different clusters;
(4) pseudotime analysis according to specific markers.
6. 良好的统计学基础：
一直以来我都忽略了统计学的重要性，直到开始看了一些小书，理解到，原来测量值并不等于真实值，我们一切的测量、实验都是而是尽可能的描绘出真实数据所形成的函数分布。而统计学则是帮助我们用合适的实验设计、统计方法尽可能的向真实情况靠拢。统计学不仅仅是统计学，也是一门人生哲学。
1. 正文开始
为了减少学习难度，我从一个专门用于单细胞数据clustering的R包的官网获取教程：Seurat
image.png
可以看到网页上有很多按钮，about是介绍R包的基本情况，install盲猜是教如何安装，Vigneetes选项则是可以看到他们根据不同需求做了不同的教程，此外还有extension package给seurat作补充，简直不要太优秀了，还想去哪里，就在这儿学吧！这已经不是第一次发现原来我们的学习资源真的太多了，而官方有时候往往教程是做得最好的，答应我，就留在这了好吗？
2. 项目介绍
image.png
选了教程1：
A basic overview of Seurat that includes an introduction to:
(1) QC and pre-processing，QC和数据预处理
(2) Dimension reduction，降维
(3) Clustering，聚类
(4) Differential expression，差异分析
3. 准备工作：数据下载，R包安装
raw data: https://s3-us-west-2.amazonaws.com/10x.files/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz
是一个tar压缩文件，作为一个什么也不懂的人，我在电脑上解压后打开看到
image.png
文件1：barcodes
文件2：ensemble ID to gene name
文件3：真正的测序数据
由此可见，barcodes和数据数据是必不可少的，而ensemble ID to gene name，想必只要做生信这个文件则是大家都有的。
R包安装和激活
# install.packages('BiocManager') # 没有的话就安装一个
# BiocManager::install('multtest')
# install.packages('Seurat')
install.packages("patchwork")
library(dplyr)
library(Seurat)
library(patchwork)
Step1 load data and structure learning
由于数据是来自于10X genomic的cell ranger，因此我们需要使用10X专门的R function去读取它们的数据Read10X
# Load the PBMC dataset 读取数据
pbmc.data <- Read10X(data.dir =
"../data/pbmc3k/filtered_gene_bc_matrices/hg19/")
以上这些操作是教程里面都有的常规操作，那么我们作为初学者，还需要去学习解构raw data，最简单的方法就是class(pbmc.data)和dim(pbmc.data)
image.png
可见pbmc.data是一个稀疏矩阵，而稀疏矩阵对我而言是陌生的。尽可能的自我学习之后我再去问徐洲更老师
image.png
image.png
我可真是好学呢
老师的网页：http://xuzhougeng.top/archives/R-Sparse-Matrix-Note
image.png
由此可见，10x pipeline出来的结果是稀疏矩阵，只记录非零的位置，同时配备有行名和列名等尽可能详尽的信息，真是优秀呢！
Step 2 创建Seurat对象，并设置条件筛选细胞
# Initialize the Seurat object with the raw (non-normalized data).
# 根据Raw data创建Seurat对象
pbmc <- CreateSeuratObject(counts = pbmc.data,
project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc
image.png
举例查看数据（如下图）
# Lets examine a few genes in the first thirty cells
pbmc.data[c("CD3D", "TCL1A", "MS4A1"), 1:30]
image.png
符号.代表述数值0的位置，由于10X的数据有很多0值，所以10X测序的结果采用了稀疏矩阵来记录结果。
Step 3 QC and pre-processing
3.1 QC的标准
（1）能检测到某个基因的细胞数，即unique基因的分布情况，对应上面的min.cells = 3；
（2）每个细胞能检测到的基因数，对应上面的min.features = 200；
（3）线粒体gene的比例要足够小，使用PercentageFeatureSet函数计算，以MT-开头的则是线粒体gene。
image.png
3.2 添加线粒体百分比列
# 用 [[符号 可以在对象的metadata 中添加列，中间则是列名
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
image.png
一个非常好的添加列的方法，并且可以看到unique gene（nFeature）和total molecular（nCount）已经算好了，放在meta.data里面。
3.3 可视化QC矩阵
数据好不好，是骡子是马，拉出来溜溜，很简单，就是把nFeature，nCount，还有MT percentage画在同一个图上，检查数据情况。
# Visualize QC metrics as a violin plot
# 使用小提琴图可视化QC 矩阵
VlnPlot(pbmc,
features = c("nFeature_RNA",
"nCount_RNA",
"percent.mt"),
ncol = 3)
image.png
可以看到，细胞的unique gene（nFeature）和total molecular（nCount）都挺集中的，分别在1000和2000左右，而线粒体基因的比例在百分之2.5左右，合格标准大概设置为多少，需要查一下。
3.4 查看不同特征之间的关系
FeatureScatter顾名思义，就是做Feature的Scatter散点图，这个命令可以用于任何你感兴趣的feature之间的interaction探索。
plot1 <- FeatureScatter(pbmc,
feature1 = "nCount_RNA",
feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc,
feature1 = "nCount_RNA",
feature2 = "nFeature_RNA")
plot1 + plot2
image.png
图上的数字-0.13，0.95代表什么尚未厘清，但猜测是相关系数，相关性之类的意思。
3.5 设定标准，筛选数据
这里教程列出的筛选标准为：测到的total molecular（nFeature）在200-2500之间，并且线粒体基因的比例要小于5%，印证了前面QC的标准。
pbmc <- subset(pbmc,
subset = nFeature_RNA > 200 &
nFeature_RNA < 2500 &
percent.mt < 5)
4. Normalizing the data
从数据集中删除不需要的单元格后，下一步是标准化数据。
默认情况下，使用全局缩放规范化方法LogNormalize，该方法通过总表达式对每个单元格的特征表达式度量进行标准化，并将其乘以一个缩放因子(默认为10,000)，然后对结果进行log转换。标准化值存储在pbmc[["RNA"]]@data中。
pbmc <- NormalizeData(pbmc,
normalization.method = "LogNormalize",
scale.factor = 10000)
# 在默认设定的情况下，上面的的代码可简写为：
pbmc <- NormalizeData(pbmc)
5. Identification of highly variable features (feature selection)
使用FindVariableFeatures完成差异分析，选择数据集中差异较高的特征基因（默认2000）并用于下游分析。
# 差异分析，需要注意selection.method有3中，请查看帮助文档学习
pbmc <- FindVariableFeatures(pbmc,
selection.method = "vst",
nfeatures = 2000)
# 标记出top 10
top10 <- head(VariableFeatures(pbmc), 10)
# 绘图
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1,
points = top10,
repel = TRUE)
plot1 + plot2
image.png
6. Scaling the data 缩放数据
6.1 缩放的标准
**接下来，应用线性变换来缩放，这是一个标准的预处理步骤，比PCA等降维技术更重要。ScaleData函数:
（1）使每个基因在所有细胞间的表达量均值为0
（2）使每个基因在所有细胞间的表达量方差为1
缩放操作给予下游分析同等的权重，这样高表达基因就不会占主导地位，其结果存储在pbmc[["RNA"]]@scale.data中**
6.2 缩放所用到的代码，以及加速方式
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
缩放是Seurat工作流程中的一个重要步骤，但仅限于将要进行PCA的基因。因此，ScaleData中的默认值只是对前面确定的变量特性(默认值为2,000)执行缩放。要做到这一点，可以忽略前面函数调用中的features参数，即pbmc <- ScaleData(pbmc)。由于Seurat heatmap(与DoHeatmap一起产生，如下图所示)，需要对heatmap中所有基因进行缩放，以确保高表达的基因不会在heatmap中占主导地位。为了确保我们在后面的热图中不会遗漏任何基因，我们将在本教程中缩放所有基因。
6.3 删除不需要的variation source
使用ScaleData函数从单细胞数据集中删除不需要的变量。例如，我们可以“regress out”与细胞周期阶段或线粒体污染相关的异质性（去除不需要的变量，即校正协变量，校正线粒体基因比例的影响）。如下：
pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")
关于ScaleData还有一个高级版的流程sctransform，感兴趣可以点击链接https://satijalab.org/seurat/v3.0/sctransform_vignette.html
7. Perform linear dimensional reduction 线性降维
对上一步骤的变量得到的缩放数据进行PCA分析，如果想要改变变量可使用参数features
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
Seurat支持多种方式来可视化PCA 数据：VizDimReduction，DimPlot，DimHeatmap
# VizDimReduction，小范围试一下
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
# DimPlot
DimPlot(pbmc, reduction = "pca")
PCA
DimHeatmap可以探索数据集的异质性来源，帮助确定用于下游分析的主成分。细胞和features都是根据它们的PCA分数排序。
DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)
PCA
可以看到PC1-6对于数据分离比较有帮助，而越往后面帮助越小，可以通过这种方式来确定，下游分析中用于分析的主成分。
8. Determine the ‘dimensionality’ of the dataset 计算数据维度
为了克服scRNA-seq数据的technical noise，Seurat基于它们的PCA分数对单元进行分组，每个PC实质上代表一个“metafeature”，它将跨相关特征集的信息组合在一起。但是到底我们该选多少个PC呢？即，我们需要确定数据集的维度
8.1 JackStraw 流程
**随机置换一部分数据（默认为1％），然后重新 PCA，重复此过程。将包含较多低 P 值特征的主成分为“显著的”主成分。
**
# 对于大的数据集，JackStraw会耗费更多时间
pbmc <- JackStraw(pbmc, num.replicate = 100) # 重复一百次
pbmc <- ScoreJackStraw(pbmc, dims = 1:20) # 选择20个PC
# 耗时3 min 49 s
JackStrawPlot函数提供了一个可视化工具，用于用均匀分布(虚线)比较每个PC的p值分布。“显著的”PCs将显示出一个低p值(虚线以上的实线)的强富集特性。在这种情况下，在最初的10-12个pc之后，重要性似乎急剧下降。
# 可视化每个PC的P value分布
JackStrawPlot(pbmc, dims = 1:15)
image.png
8.2 Elbow plot 弯头图法？
根据每个组成部分的方差百分比对主要组成部分的方差结实率进行排序(弯头图函数)。在这个例子中，我们可以观察到PC9-10周围的一个“弯头”，这表明大部分真实信号是在前10个pc中捕获的。
ElbowPlot(pbmc)
image.png
官方补充：识别数据集的真实维度-对于用户来说是有挑战性的/不确定的。因此，建议考虑以下三种方法：
（1）more supervised，探索PCs以确定相关的异质性来源，例如可以与GSEA一起使用。
（2）基于随机null模型的统计测试，但是对于大型数据集来说非常耗时，并且可能不会返回一个明确的PC cutoff。
（3）常用的heuristic算法，可以立即计算出来。
在这个例子中，所有这三种方法都产生了相似的结果，有理由选择PC 7-12之间的任何一个作为cutoff。
在这里教程选择了10个PC，但是提醒用户考虑：
（1）树突状细胞和NK细胞爱好者可能认识到，与PCs 12和PCs 13密切相关的基因定义了罕见的免疫亚群(即MZB1是浆细胞样dc的标记)。但是，这些组非常罕见，在没有预先知识的情况下，很难从这种大小的数据集的背景噪声中区分出来。
（2）鼓励用户使用不同数量的pc(10、15甚至50!)重复下游分析。正如你将观察到的，结果往往没有显著的不同。
（3）建议用户将这个参数设置偏高。例如，仅使用5个PCs执行下游分析会显著影响结果。
9. Cluster the cells 细胞聚类
这些方法将细胞嵌入到一个图结构中——例如K-nearest neighbor (KNN)图，在具有相似特征表达模式的单元格之间绘制边缘，然后尝试将这个图划分为高度互连的“准派系”或“社区”。
（1）先在PCA空间中构造一个基于欧氏距离的KNN图，然后根据任意两个细胞在局部邻近区域的共享重叠(Jaccard相似性)来细化它们之间的边权值。此步骤使用FindNeighbors函数执行，并将之前定义的数据集维度(前10个pc)作为输入。
（2）为了对细胞进行聚类，应用模块化优化技术，如Louvain算法(默认)或SLM，以迭代方式将单元分组在一起，目标是优化标准模块化函数。FindClusters函数实现这个过程，并包含一个分辨率参数，该参数设置下游集群的“granularity”，增加的值将导致更多的集群。将该参数设置在0.4-1.2之间，对于3K左右的单细胞数据集通常会得到良好的结果。对于较大的数据集，最佳分辨率通常会增加。可以使用Idents函数找到clusters。(其实我看不懂)
pbmc <- FindNeighbors(pbmc, dims = 1:10) # 前10个PC
pbmc <- FindClusters(pbmc, resolution = 0.5) # 分辨率是0.5
# 结果
# Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
# Number of nodes: 2638
# Number of edges: 95930
# Look at cluster IDs of the first 3 cells
head(Idents(pbmc), 3)
# AAACATACAACCAC AAACATTGAGCTAC AAACATTGATCAGC
# 0 3 2
# Levels: 0 1 2 3 4 5 6 7 8
10. Run non-linear dimensional reduction (UMAP/tSNE) 非线性降维
Seurat提供了几种非线性的降维技术，如tSNE和UMAP，以可视化和探索这些数据集。这些算法的目标是学习数据的底层流形，以便在低维空间中把相似的单元放在一起。上面所确定的基于图的集群中的单元应该在这些降维图上共同定位。对于UMAP和tSNE的输入数据，我们建议使用相同的PCs作为聚类分析的输入。
# 安装UMAP
reticulate::py_install(packages ="umap-learn")
pbmc <- RunUMAP(pbmc, dims = 1:10) # 102.7 Mb
DimPlot(pbmc, reduction = "umap", pt.size=1)
# 可视化，用umap的方法
# 保存数据
saveRDS(pbmc, file = "../output/pbmc_tutorial.rds")
image.png
11. Finding differentially expressed features (cluster biomarkers) 差异基因
11.1 解释标准或参数
Seurat可以找到通过差异表达式定义集群的标记。默认情况下，它识别单个簇的阳性和阴性标记(在ident1中指定)，与所有其他细胞相比较。findallmarker为所有集群自动化这个过程，但是您也可以测试集群组之间的相互关系，或者测试所有细胞。
(1) min.pct参数要求在两组单元中至少检测到一个特征；(2) thresh.test参数要求一个特性在两组之间有一定的差异(平均)；
若两个值都设置为0，会非常耗时——因为这将测试大量不太可能具有高度差异性的features。
(3) 为了加速，可以设置max.cells.per.ident参数，这将向下采样每个标识类，使其单元格数不超过设置的单元格数。虽然通常会有功率的损失，速度的增长可能是显著的，最高度差异表达的特征可能仍然会上升到顶部。
11.2 举例
# 找cluster 1的markers
cluster1.markers <- FindMarkers(pbmc,
ident.1 = 1, # 至少检测到1个feature
min.pct = 0.25) # 设置后会加速，默认是0.1
# 耗时23 s
head(cluster1.markers, n = 5) # 显示！
image.png
# find all markers distinguishing cluster 5 from clusters 0 and 3
cluster5.markers <- FindMarkers(pbmc,
ident.1 = 5, # 找cluster5的marker
ident.2 = c(0, 3), # 和cluster 0-3对比
min.pct = 0.25)
head(cluster5.markers, n = 5)
image.png
# find markers for every cluster compared to all remaining cells,
# report only the positive ones
pbmc.markers <- FindAllMarkers(pbmc,
only.pos = TRUE,
min.pct = 0.25,
logfc.threshold = 0.25)
pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)
image.png
11.3 更加精细的寻找marker的方法
我们可以用 test.use 参数设置多种检验方式。例如，用 ROC 返回任何单个标记的“分类能力”。
# 用ROC方法
cluster1.markers <- FindMarkers(pbmc,
ident.1 = 0,
logfc.threshold = 0.25,
test.use = "roc",
only.pos = TRUE)
11.4 小提琴图可视化marker
画法1
# 小提琴图将marker可视化，随便挑上部分结果中的3个基因
# LDHB, CCR7,LEF1可以将cluster0分别于其他cluster
VlnPlot(pbmc, features = c("LDHB", "CCR7","LEF1"))
image.png
画法2
# raw counts也可以画
VlnPlot(pbmc,
features = c("LDHB", "CCR7","LEF1"),
slot = "counts",
log = TRUE)
2
画法3
# cluster 图展示
FeaturePlot(pbmc,
features = c("LDHB", "CCR7","LEF1"))
image.png
画法4
DoHeatmap为给定的细胞和特征生成一个表达式heatmap。在本例中，绘制每个集群的前20个标记(如果小于20，则绘制所有标记)。
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()
image.png
12. Assigning cell type identity to clusters
Fortunately in the case of this dataset, we can use canonical markers to easily match the unbiased clustering to known cell types:
image.png
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 = 1.2,
label.size=5) + NoLegend()
image.png
收尾
感觉收获很大，至少跑完了一个流程，在心里有了大致的印象，在自学的过程中我用的最多的求助不是问别人，而是?Rfunction查询帮助文档，理解函数有助于加深我对数据转换的理解。整个学习当然少不了很多小伙伴的帮助，徐洲更老师是给我解释了R语言的稀疏矩阵，卖萌哥和小洁，让我明白%>%原来是R语言中的管道符，跟linux中的|是一样的含义，这个看似不起眼的point，完全阻碍了我对数据处理的理解。当然我也有参考生信菜鸟团中的同一个教程，想对于那个教程的简练和精要，我的文字更加生硬和生疏，还有非常大的进步空间。
整一个流程增加了我对单细胞数据处理的的了解，我原先的理解是：QC--filtration--cluster--GO--identify marker
现在看来：
（1）QC的有很多门路，可以变着法的做feature-feature interaction，同时，线粒体基因比例则是我第一次听到的，相信我不知道的东西还有很多很多；
（2）从文献阅读的经验来看，筛选一般是对细胞测序基因数的筛选，没想到线粒体基因比例前面已经提到，但让我意想不到的是，原来droplet竟然会有dowble和mutiple的情况发生，可见对建库流程和质控的不了解；
（3）数据缩放中的校正协变量亦是针对线粒体gene；
（4）第一次认识到PCA还有这么多做法；
（5）确定数据维度，应该是最难的部分之一，对于PC数量的选择，将对后续分析具有决定性作用，而PC数目的确定也是有多种不同的方法；
（6）细胞聚类则是依据上一步PCA分析的结果，选择合适的PC数目和cutoff。
感谢（改天再检查错别字）！
