建立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)