自学Seurat进阶——官网教程代码手把手解读（上）

建立Seurat对象

官网中的tutorial使用的是10x公司测序、通过cell ranger初步处理后的pbmc3k的数据，下载的示例数据链接：https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz

解压后发现文件夹包括三个文件barcodes.tsv, features.tsv, matrix.mtx。

将整个文件夹用Read10X()函数读入，会得到一个UMIcount矩阵，接着使用这个矩阵创建Seurat对象，这个对象包含很多slot，即有关对象的属性信息，不仅包括count matrix，还包括后续进行PCA，聚类等分析不断积累更新的信息。

上一篇教程中自学Seurat的一点记录_weixin_54709098的博客-CSDN博客讲过，为了节约存储空间，使用的是稀疏矩阵，矩阵行名为基因名，列名为细胞名（barcode）。

由于单细胞测序的对象——每个单细胞都没有replicate，且限于样本大小，会导致很多dropout的存在，因此对象建立好后，需进行预处理过滤掉低质量的细胞和基因，留下最能反映细胞间差异的信息，这也是后续降维的核心目标。

预处理

质量控制和选择用于后续分析的细胞

我们主要通过控制每个细胞探测到的基因数和线粒体基因比值进行QC

基因数过高可能是液滴中包含了双细胞或多细胞

基因数过低可能是空液滴，或者RNA泄露的低质量细胞，通常受损伤的低质量细胞线粒体应激，其中的RNA外漏，因此线粒体基因比例偏高，这两个指标有一定的负相关性，它们的负相关性和基因数与count数之间的正相关性都可以通过feature scatter体现且会计算相关系数。

线粒体基因比例用下面代码实现

# The [[ operator can add columns to object metadata. This is a great place to stash QC stats
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

# Visualize QC metrics as a violin plot
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

# FeatureScatter is typically used to visualize feature-feature relationships, but can be used
# for anything calculated by the object, i.e. columns in object metadata, PC scores etc.

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)

 

nCount_RNA, nFeature_RNA, percent.mt这三个QC Metrics被存储在meta.data中。

用VlnPlot，Feature Scatter对这些QC Metrics可视化，用subset函数取符合标准的数据。

Normalize数据

由于不同的细胞测序深度不同。所以导致每个细胞的总reads数不同，简单来说假如某个基因在某几个细胞中高表达，我们倾向于将它们归为一类，但我们如何去定义高表达，显然是该基因表达量相对于该细胞中所有基因的表达量的比值，是一个相对值，而不是绝对值，因此我们用normalize来消除不同细胞测序深度的影响。Log Normalize将基因的相对表达值乘以10000，默认的换算系数，再进行对数化。10000是默认的平均每个细胞的总reads数，而log可以有效的降低数据的离散性，如对10和1000相差100倍的两个数取对数，则得到1和3，大大降低了离散程度。利于数据的分析，为了解决0无法取对数的问题，将每个数+1。

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

高变基因的识别

高变基因就是在一些细胞中高表达，另一些细胞中低表达，它们可以充分反映细胞间的差异，高变基因是通过标准差（Standard deviation）来体现的。标准差越大，越高变。

pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
#返回2000个基因，用于后续的降维聚类

# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(pbmc), 10)

# plot variable features with and without labels
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2

Scale 数据

scaling是一种线性转化，即将某基因在不同细胞中的分布转化成均值为0，方差为1的分布，进行z-score的转换。这使每个基因在下游分析中的具有相同的权重，高表达基因不那么显著。

pbmc <- ScaleData(pbmc)