本文首发于“生信补给站”公众号 https://mp.weixin.qq.com/s/hoFgGZRXCp0NnumALISVSw
学习单细胞转录组肯定先来一遍Seurat官网的标准流程。
数据来源于Peripheral Blood Mononuclear Cells (PBMC)，共2700个单细胞， Illumina NextSeq 500平台。下载链接在这：https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz
解压到目录下，有以下三个文件，也是10X的标准文件
barcodes.tsv ,genes.tsv, matrix.mtx
一 加载R包 数据集
Seurat可以直接用Read10X函数读取cellranger的结果数据，使用pbmc数据初始化Seurat对象
library(dplyr) library(Seurat) library(patchwork) # Load the PBMC dataset pbmc.data <- Read10X(data.dir = "./data/hg19/") #解压缩后的路径 # Initialize the Seurat object with the raw (non-normalized data). pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200) pbmc
查看count matrix矩阵
#查看数据 dim(pbmc.data) # 32738 2700 # Lets examine a few genes in the first thirty cells pbmc.data[c("CD3D", "TCL1A", "MS4A1"), 1:30] 3 x 30 sparse Matrix of class "dgCMatrix" [[ suppressing 30 column names ‘AAACATACAACCAC-1’, ‘AAACATTGAGCTAC-1’, ‘AAACATTGATCAGC-1’ ... ]] CD3D 4 . 10 . . 1 2 3 1 . . 2 7 1 . . 1 3 . 2 3 . . . . . 3 4 1 5 TCL1A . . . . . . . . 1 . . . . . . . . . . . . 1 . . . . . . . . MS4A1 . 6 . . . . . . 1 1 1 . . . . . . . . . 36 1 2 . . 2 . . . .
其中的.表示0，即no molecules detected。seurat使用一个稀疏矩阵来保存count matrix，节省存储空间。
二 数据预处理
2.1 QC
一般使用以下三个标准，也可以参考commonly used
- 每个细胞中检测到的唯一基因数
- 低质量的细胞或者空的droplet液滴通常含有很少的基因
- Cell doublets 或 multiplets 可能表现出异常高的基因计数
- 每个细胞中检测到的分子总数
- 线粒体基因含量比例
- 低质量或者死亡细胞含有很高的线粒体基因
- 使用PercentageFeatureSet()计算线粒体QC比例
- MT-开头的基因是线粒体基因
# The [[ operator can add columns to object metadata. This is a great place to stash QC stats pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-") # Show QC metrics for the first 5 cells 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.2244898
注意PercentageFeatureSet函数可以计算每个CELL中的指定基因子集的计数百分比。
小提琴图：查看基因数目, UMI数目, 线粒体基因占比
# Visualize QC metrics as a violin plot
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
基因数目, 线粒体基因占比与UMI数目相关性
# 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
过滤
- 过滤具有UMI计数超过 2500 或小于 200的细胞
- 过滤具有>5%线粒体的细胞
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
2.2 数据标准化
默认标准化方法为LogNormalize，标化后的数据存在pbmc[["RNA"]]@data中。
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000) #pbmc <- NormalizeData(pbmc)
如果我记不住这些[[ 和 @ 怎么办？
使用最“简单，通用”的方式，str(pbmc) 查看一下数据结构，然后根据结构对应使用@ 或者 $ 。
str(pbmc)
#截取部分展示用
Formal class 'Seurat' [package "SeuratObject"] with 13 slots ..@ assays :List of 1 .. ..$ RNA:Formal class 'Assay' [package "SeuratObject"] with 8 slots .. .. .. ..@ counts :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots .. .. .. .. .. ..@ i : int [1:2238732] 29 73 80 148 163 184 186 227 229 230 ... .. .. .. .. .. ..@ p : int [1:2639] 0 779 2131 3260 4220 4741 5522 6304 7094 7626 ... .. .. .. .. .. ..@ Dim : int [1:2] 13714 2638 .. .. .. .. .. ..@ Dimnames:List of 2 .. .. .. .. .. .. ..$ : chr [1:13714] "AL627309.1" "AP006222.2" "RP11-206L10.2" "RP11-206L10.9" ... .. .. .. .. .. .. ..$ : chr [1:2638] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ... .. .. .. .. .. ..@ x : num [1:2238732] 1 1 2 1 1 1 1 41 1 1 ... .. .. .. .. .. ..@ factors : list() .. .. .. ..@ data :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots .. .. .. .. .. ..@ i : int [1:2238732] 29 73 80 148 163 184 186 227 229 230 ... .. .. .. .. .. ..@ p : int [1:2639] 0 779 2131 3260 4220 4741 5522 6304 7094 7626 ... .. .. .. .. .. ..@ Dim : int [1:2] 13714 2638 .. .. .. .. .. ..@ Dimnames:List of 2 .. .. .. .. .. .. ..$ : chr [1:13714] "AL627309.1" "AP006222.2" "RP11-206L10.2" "RP11-206L10.9" ... .. .. .. .. .. .. ..$ : chr [1:2638] "AAACATACAACCAC-1" "AAACATTGAGCTAC-1" "AAACATTGATCAGC-1" "AAACCGTGCTTCCG-1" ... .. .. .. .. .. ..@ x : num [1:2238732] 1.64 1.64 2.23 1.64 1.64 ... .. .. .. .. .. ..@ factors : list() .. .. .. ..@ scale.data : num[0 , 0 ] .. .. .. ..@ key : chr "rna_" .. .. .. ..@ assay.orig : NULL .. .. .. ..@ var.features : chr(0) .. .. .. ..@ meta.features:'data.frame': 13714 obs. of 0 variables .. .. .. ..@ misc : list() ..@ meta.data :'data.frame': 2638 obs. of 5 variables: .. ..$ orig.ident : Factor w/ 1 level "pbmc3k": 1 1 1 1 1 1 1 1 1 1 ... .. ..$ nCount_RNA : num [1:2638] 2419 4903 3147 2639 980 ... .. ..$ nFeature_RNA: int [1:2638] 779 1352 1129 960 521 781 782 790 532 550 ... .. ..$ percent.mt : num [1:2638] 3.02 3.79 0.89 1.74 1.22 ... .. ..$ percent.HB : num [1:2638] 0 0 0 0 0 0 0 0 0 0 ... 根据层级结构找到 data 即可， pbmc@assays$RNA@data （标准化后的数据）
另：pbmc@assays$RNA@counts为原始count数据
简单看一下标准化前后的数据
par(mfrow = c(1,2)) hist(colSums(pbmc$RNA@counts),breaks = 50) hist(colSums(pbmc$RNA@data),breaks = 50)
2.3 高变基因（特征选择）
选择数据集中高变异的特征子集（在某些细胞中高表达，在其他细胞中低表达）。通过Seurat内置的FindVariableFeatures()函数，计算每一个基因的均值和方差，默认选择高变的2000个基因用于下游分析。
标示Top10的基因 以及 标示目标基因
#标示感兴趣的基因
plot3 <- LabelPoints(plot = plot1, points = c(top10,"HLA-DPB1","CCL4"), repel = TRUE) plot2 + plot3
2.4 归一化
使用ScaleData()函数线性转换（"归一化"）数据，使得每一个基因在所有cell中的表达均值为0，方差为1.
结果在pbmc[["RNA"]]@scale.data #同样可以通过str的方式根据数据结构获取。
all.genes <- rownames(pbmc) pbmc <- ScaleData(pbmc, features = all.genes)
有需要的话也可以移出不必要的来源数据？
pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")
三 PCA 降维
3.1 使用scale后的数据进行PCA分析，默认表达变化大的基因。可以使用features参数重新定义数据集。
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc)) 可视化方式包含以下几种，VizDimReduction(), DimPlot() 和 DimHeatmap() # Examine and visualize PCA results a few different ways 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, GZMA, GZMB PC_ 3 Positive: HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPA1 Negative: PPBP, PF4, SDPR, SPARC, GNG11 PC_ 4 Positive: HLA-DQA1, CD79B, CD79A, MS4A1, HLA-DQB1 Negative: VIM, IL7R, S100A6, S100A8, IL32 PC_ 5 Positive: GZMB, FGFBP2, S100A8, NKG7, GNLY Negative: LTB, IL7R, CKB, MS4A7, RP11-290F20.3
3.2 可视化展示
1）dims指定展示几个PCs ， nfeatures指定每个PC展示多少基因
VizDimLoadings(pbmc, dims = 1:2, nfeatures = 20, reduction = "pca")
2） PCA点图
DimPlot(pbmc, reduction = "pca")
3） PCA热图
DimHeatmap(pbmc, dims = 1:4, #几个PCs cells = 600, #多少cell ncol = 2, #图形展示几列 balanced = TRUE)
