随着单细胞技术的发展和普及，越来越多的课题组开始使用单细胞技术进行自己的课题研究。然而，在单细胞的注释上，商业公司提供的流程化注释往往难以满足实际的需求，甚至得到牛头不对马嘴的结果。这期讲带大家实现全方位的单细胞注释指南。
1. 引言
单细胞RNA测序（scRNA-seq）是一种先进的分子生物学技术，可以在单个细胞水平上评估基因表达。在过去的几年里，scRNA-seq在细胞注释中扮演了重要角色，它已经成为揭示细胞异质性和发展过程中转录组动态的有力工具。目前，常见的scRNA-seq注释方法包括基于表达模式的细胞标记、参考基因表达和聚类特征等。
然而，尽管scRNA-seq可以在单细胞的水平上揭示基因的表达，但是目前在细胞注释上依然存在不小的挑战。首先，技术本身的特点如扩增偏差、噪音和低覆盖率等，可能导致结果的偏差和不确定性。其次，数据的处理和分析也是一个复杂而关键的步骤，如何准确地鉴定和分类不同的细胞类型仍然是一个挑战。此外，由于细胞的异质性和动态性，以及在特定生理状态下转录组的变化，细胞状态鉴定和动态分析也是进一步研究的难点。
2.常见注释策略
细胞注释按照是否有人工参与可以分为自动注释和人工注释。
2.1 自动注释
最为常见的就是使用SingleR进行，其根据参考数据集的表达矩阵，和输入数据表达矩阵的相似性来识别细胞类型。其使用非常的方便，是绝大多数商业公司提供的全自动流程方案。然而，缺点是显而易见的。首先，SingleR包含的参考数据集有限，难以满足实际样本的需求，参考的有限导致注释出现牛头不对马脚的结果。
另一类策略则是根据Marker基因进行，这里常用的软件有cellassign和scCATCH等，他们的思路是类似的，首先需要我们输入参考的marker基因，marker基因可以从cell marker2.0、sctype或已经发表文献等获得。这里推荐大家使用scCATCH，其自带了cell marker的数据库对新手较为友好。通过marker基因进行细胞分群的鉴定可以特异性的根据物种、组织的特点，注释到更为合理的结果。然而，这也并非“万全之策”，由于样本之间、组织之间甚至是批次之间，分析策略之间等一系列因素的差异，适用于别人的marker不见得就一定适用于自己的样本。
2.2 人工注释
人工注释，顾名思义，就是自己去注释。如何自己去注释呢？手动注释在很大程度上比较类似基于marker基因注释的方法。基于参考文献和以往经验得到的marker基因的表达丰度进行。
3.实战
3.1 SingleR自动注释
使用的数据是经典的pbmc3k （https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz）
安装好 R并安装R包
Seurat、SingleR、tidyverse等。
首先我们导入数据，pbmc解压出来的数据是包括三个文件barcodes.tsv genes.tsv matrix.mtx。首先是进行导入数据 -> 质控处理 -> 标准化 -> 降纬 -> 分群
library(Seurat)
library(tidyverse)
library(SingleR)
pbmc <- Read10X('~/pbmc/filtered_gene_bc_matrices/hg19') #这里填写解压后的目录
pbmc <- CreateSeuratObject(pbmc) #转化为Seurat对象
Pattern = '^MT-' #线粒体基因的名字，根据实际去匹配，人是^MT-开头，小鼠^mt-，大鼠^Mt-
pbmc[["percent.mt"]] <- PercentageFeatureSet(object = pbmc, pattern = Pattern)
pbmc <- subset(x = pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 8000 & percent.mt < 15) #进行细胞过滤，这里是要求每个细胞里面至少有200个的基因，但不大于8000个基因，线粒体基因含量不大于15%。这里的过滤就根据实际来，没有绝对的标准。
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000) #进行log归一化
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000) #寻找高变基因
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes) #根据基因数进行数据的缩放
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc)) #进行PCA降纬，特征是我们前面找到的高变基因
pbmc <- FindNeighbors(pbmc, dims = 1:30) #寻找临近的特征点
pbmc <- FindClusters(object = pbmc, verbose = T, resolution = 0.7) #分群，关键在于resolution调整分辨率，在0～1之间，约大就分越多分群
这样我们就得到就将细胞进行了初步的分群，随后开始导入SingleR的参数数据集。
也有几种方法，一种我们手动下载，一种是安装``celldex```包导入， 如 celldex::HumanPrimaryCellAtlasData，或者导入别人已经有注释的rds,方法也非常简单。只需要
ref <- celldex::HumanPrimaryCellAtlasData() #用celldex包导入，第一次需要一点时间下载，该数包还有其他更多的参考数据集，可以详细看手册
ref <- readRDS('pbmc.rds') #用别人注释好的rds，需要先准备，没有准备的话用方法一
pbmc_anno <- SingleR(test = pbmc@assays$RNA@data, ref = ref, labels = ref$label.main, cluster = pbmc$seurat_clusters) #这里我们选择按照分群去注释
这个时候我们可以看到，已经注释好了。主要有T细胞和B细胞，单核细胞和血小板等我们进行一下可视化展示
image.png
pbmc <- RunUMAP(object = pbmc, dims = 1:30) #首先进行UMAP降纬
pbmc$CellType <- pbmc_anno$labels[match(pbmc$seurat_clusters, row.names(pbmc_anno))] #加上SingleR的注释结果
p1 <- DimPlot(pbmc, group.by = "seurat_clusters",reduction = "umap", label = T)
p2 <- DimPlot(pbmc, group.by = "CellType",reduction = "umap", label = T)
p1+p2
ggsave('CellAnnot_umap.png', height = 4.5, width = 12, dpi = 300)
CellAnnot_umap.png
我们对比一下官方给出的手动注释结果，可以看到已经是基本接近了。
image.png
3.2 scCATCH注释
这是我个人更推荐的一种方法，首先主要安装 scCATCH包，我们接着上面的代码。
library(scCATCH)
obj <- createscCATCH(data = pbmc@assays$RNA@data, cluster = as.character(pbmc$seurat_clusters) ) #创建CATCH对象
obj <- findmarkergene(obj, marker = cellmatch, if_use_custom_marker = F ,use_method = "2",logfc = 0.4, cell_min_pct = 0.25, pvalue = 0.05, cancer = 'Normal', species = 'Human', tissue = 'Blood') #这里需要注意的关键参数是marker我们使用了自带的cellmatch， 指定好癌症的类型，物种和组织，我们可以通过unique(cellmatch$tissue)去查看自带的组织等。
obj <- findcelltype(obj)
obj@celltype$cell_type
我们看看效果
marker <- tibble( cluster = parse_number(obj@celltype$cluster), cell_tpye = obj@celltype$cell_type)
pbmc$CellType <- marker$cell_tpye[match(pbmc$seurat_clusters,marker$cluster)]
p1 <- DimPlot(pbmc, group.by = "seurat_clusters",reduction = "umap", label = T)
p2 <- DimPlot(pbmc, group.by = "CellType",reduction = "umap", label = T,label.size=2)
p1+p2
ggsave('CellAnnot_umap2.png', height = 4.5, width = 12, dpi = 300)
CellAnnot_umap2.png
效果也是很不错。那么如果我们需要参考别人的结果如何进行呢？这里还是只需要借助scCATCH包即可。我们需要准备以下样式的表格，包含以下的行头，比如我们使用pbmc3k提供的marker
celltype gene subtype1 subtype2 subtype3 pmid
| Cluster ID | Markers | Cell Type |
|---|---|---|
| 0 | IL7R, CCR7 | Naive CD4+ T |
| 1 | CD14, LYZ | CD14+ Mono |
| 2 | IL7R, S100A4 | Memory CD4+ |
| 3 | MS4A1 | B |
| 4 | CD8A | CD8+ T |
| 5 | FCGR3A, MS4A7 | FCGR3A+ Mono |
| 6 | GNLY, NKG7 | NK |
| 7 | FCER1A, CST3 | DC |
| 8 | PPBP | Platelet |
对应代码如下：
cell_marker <- tribble(
~celltype, ~gene, ~subtype1,~subtype2,~subtype3, ~pmid,
"T cell", "IL7R", "Naive CD4+", NA, NA, NA,
"T cell", "CCR7", "Naive CD4+", NA, NA, NA,
"Mono", "CD14", "CD14+", NA, NA, NA,
"Mono", "LYZ", "CD14+", NA, NA, NA,
"T cell", "IL7R", "Memory CD4+", NA, NA, NA,
"T cell", "S100A4", "Memory CD4+", NA, NA, NA,
"B cell", "MS4A1", NA, NA, NA, NA,
"T cell", "CD8A", "CD8+", NA, NA, NA,
"Mono", "FCGR3A", "FCGR3A+", NA, NA, NA,
"Mono", "MS4A7", "FCGR3A+", NA, NA, NA,
"NK", "GNLY", NA, NA, NA, NA,
"NK", "NKG7", NA, NA, NA, NA,
"DC", "FCER1A", NA, NA, NA, NA,
"DC", "CST3", NA, NA, NA, NA,
"Platelet", "PPBP", NA, NA, NA, NA,
)
obj <- findmarkergene(obj, marker = cell_marker, if_use_custom_marker = T , use_method = "1",logfc = 0.4, cell_min_pct = 0.25, pvalue = 0.05)
obj <- findcelltype(obj)
marker <- tibble( cluster = parse_number(obj@celltype$cluster), cell_tpye = obj@celltype$cell_type)
pbmc$CellType <- marker$cell_tpye[match(pbmc$seurat_clusters,marker$cluster)]
p1 <- DimPlot(pbmc, group.by = "seurat_clusters",reduction = "umap", label = T)
p2 <- DimPlot(pbmc, group.by = "CellType",reduction = "umap", label = T,label.size=2)
p1+p2
ggsave('CellAnnot_umap3.png', height = 4.5, width = 12, dpi = 300)
CellAnnot_umap3.png
这样我们就得到了一个和官方给到的注释完美一致的了。而我们如果需要手动在修改只需要修改marker中的数据就实现了所谓的手动注释了。
网友评论