单细胞分析实录(7): 差别表达分析/细胞类型注释

时间 2021-01-08

标签 ide 函数学习 code orm 对象 blog element get 栏目 Java 繁體版

原文原文链接

前面已经讲解了：
单细胞分析实录(1): 认识Cell Hashing
单细胞分析实录(2): 使用Cell Ranger获得表达矩阵
 单细胞分析实录(3): Cell Hashing数据拆分
 单细胞分析实录(4): doublet检测
 单细胞分析实录(5): Seurat标准流程
 单细胞分析实录(6): 去除批次效应/整合数据ide

这一节咱们进入细胞类型注释的学习，总的来讲，两条路线，手动注释和软件注释。我在实际处理的过程当中，主要仍是手动注释，软件的注释结果只做为一个参考，来辅助证明手动注释的结果是准确无误的。
相信看过几篇单细胞文献的小伙伴基本都会知道几种常见细胞大类的marker，咱们在注释的时候用这些marker就能够基本肯定细胞大类了。另外，以SingleR为例，对于细胞大类的注释还算准确，固然也没有到很准确的程度，我试过更换SingleR里面不一样的参考集，最后获得的大类注释结果一致性不到80%。对于细胞小类的注释，软件就更加没法胜任了，我几乎没有看太高分的文献会用SingleR的小类注释结果。固然不排除随着单细胞研究的普及和深刻，之后能有更准确的软件出现。
接下来，我以上一节的数据为例，走一遍个人分析流程。函数

以前的数据呈现出了16个cluster，至于resolution参数的选择，个人原则是能在降维图上分开的细胞亚群，能有它本身的cluster label，而且成团较好，较紧致的一群细胞没有必要再强行分群（好比上图的第4群）。
这16个cluster不必定都会用到，由于doublet、细胞数太少等缘由，可能还得去掉。
我推荐用常见的marker先区分一下，大概能知道cluster对应什么类型。这里用到的marker都是在文献里面常常出现的，我本身也整理了一份更全一点的细胞类型marker清单，有须要的能够在公众号后台说明。学习

celltype_marker=c(
  "EPCAM",#上皮细胞 epithelial
  "PECAM1",#内皮细胞 endothelial
  "COL3A1",#成纤维细胞 fibroblasts
  "CD163","AIF1",#髓系细胞 myeloid
  "CD79A",#B细胞
  "JCHAIN",#浆细胞 plasma cell
  "CD3D","CD8A","CD4",#T细胞
  "GNLY","NKG7",#NK细胞
  "PTPRC"#免疫细胞
)

VlnPlot(test.seu,features = celltype_marker,pt.size = 0,ncol = 2)
ggsave(filename = "marker.png",device = "png",width = 44,height = 33,units = "cm")

结合对应的关系，初步确认细胞类型以下，须要说明的是，NK细胞和T细胞不是很能区分开，一些文献也是直接把这两个当作一个大群来作的，此处我根据GNLY基因肯定第5个cluster为NK细胞，是否正确咱们后面再看。code

#免疫细胞
0: B_cell
4: Plasma_cell
1 2: T_cell
5: NK_cell
13: Unknown

#非免疫细胞
3 6 7 8 10 11 12: Epithelial
14: Endothelial
9: Fibroblasts

#doublet
15: Doublet (Myeloid+CD4)

以上根据经典的marker初步肯定了细胞大类，接着咱们找差别基因，看看找出来的每个cluster的差别基因是否是和前面鉴定的类型一致。orm

1. 找差别基因

markers <- FindAllMarkers(test.seu, logfc.threshold = 0.25, min.pct = 0.1, 
                          only.pos = TRUE, test.use = "wilcox")
markers_df = markers %>% group_by(cluster) %>% top_n(n = 500, wt = avg_logFC)
write.table(markers_df,file="test.seu_0.5_logfc0.25_markers500.txt",quote=F,sep="\t",row.names=F,col.names=T)

FindAllMarkers()这个函数会将cluster中的某一类和剩下的那些类比较，来找差别基因。与其对应的有个FindMarkers()函数，能够指定哪两个cluster对比。logfc.threshold表示logfc的阈值，这里有两个地方须要注意：一是Seurat里面的logfc计算公式很特别，并非咱们日常在bulk里面那样算均值，相除，求log，但其实也不要纠结怎么算的，只需知道这是表示倍数的一个指标便可；二是若是想画火山图，这个阈值能够设为0，否则最后画出来的火山图中间会缺一段，咱们彻底能够把所有基因先拿出来，画火山图，再根据p value, logFC这些阈值本身过滤。对象

min.pct表示基因在多少细胞中表达的阈值，only.pos = TRUE表示只求高表达的基因，test.use表示检测差别基因所用的方法。
第2行代码，根据avg_logFC这一列把前500个差别基因取出来，小于500个时，有多少保留多少。
事实上，没必要在乎咱们保留多少个差别基因，实际用到的，也就前面十几（几十）个基因。
最终获得的差别基因list以下：blog

> head(markers_df)
# A tibble: 6 x 7
# Groups:   cluster [1]
  p_val avg_logFC pct.1 pct.2 p_val_adj cluster gene    
  <dbl>     <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>   
1     0      3.58 0.999 0.114         0 0       HLA-DRA 
2     0      2.72 1     0.331         0 0       CD74    
3     0      2.59 0.983 0.169         0 0       HLA-DPA1
4     0      2.50 0.985 0.148         0 0       HLA-DPB1

而后，咱们能够根据每个cluster排在前面的基因验证前面鉴定的结果是否正确。我把刚刚获得的的表格看一下，基本符合前面的鉴定。element

2. SingleR注释

我刚开始用SingleR是在19年的年末，如今使用方法可能有一些变化了。由于Single在注释细胞的过程当中，会用到一些参考数据集，咱们能够把这些数据集保存下来，下一次使用就能够直接加载，省去了从新下载参考集的时间。get

library(celldex)
hpca.se <- HumanPrimaryCellAtlasData()
hpca.se
## class: SummarizedExperiment 
## dim: 19363 713 
## metadata(0):
## assays(1): logcounts
## rownames(19363): A1BG A1BG-AS1 ... ZZEF1 ZZZ3
## rowData names(0):
## colnames(713): GSM112490 GSM112491 ... GSM92233 GSM92234
## colData names(3): label.main label.fine label.ont

save(hpca.se,file="/ref/singler/HPCA/hpca.se.RData")

注意参考集里面是logcounts矩阵，后面对于单细胞数据集也要作相似的处理。
导入UMI count矩阵it

library(SingleR)
library(scater)
library(SummarizedExperiment)
test.count=as.data.frame(test.seu[["RNA"]]@counts)

导入参考集，保证两个数据集的基因相同，而后log处理

load(file="hpca.se.RData")
common_hpca <- intersect(rownames(test.count), rownames(hpca.se))
hpca.se <- hpca.se[common_hpca,]
test.count_forhpca <- test.count[common_hpca,]
test.count_forhpca.se <- SummarizedExperiment(assays=list(counts=test.count_forhpca))
test.count_forhpca.se <- logNormCounts(test.count_forhpca.se)

接下来是注释步骤，在这一步里，我只用了main大类注释，还有一个fine小类注释，我就不演示了，由于我以为小类注释不太准。

###main
pred.main.hpca <- SingleR(test = test.count_forhpca.se, ref = hpca.se, labels = hpca.se$label.main)
result_main_hpca <- as.data.frame(pred.main.hpca$labels)
result_main_hpca$CB <- rownames(pred.main.hpca)
colnames(result_main_hpca) <- c('HPCA_Main', 'CB')
write.table(result_main_hpca, file = "HPCA_Main.txt", sep = '\t', row.names = FALSE, col.names = TRUE, quote = FALSE) #保存下来，方便之后调用

获得的结果是这样的，每一个CB都有一个label

> head(result_main_hpca)
         HPCA_Main                 CB
1           B_cell A_AAACCCAAGGGTCACA
2          T_cells A_AAACCCAAGTATAACG
3 Epithelial_cells A_AAACCCAGTCTCTCAC
4           B_cell A_AAACCCAGTGAGTCAG
5          T_cells A_AAACCCAGTGGCACTC
6          T_cells A_AAACGAAAGCCAGAGT

咱们接下来要把这个结果整合到test.seu@meta.data中，而后画tsne/umap展现一下

test.seu@meta.data$CB=rownames(test.seu@meta.data)
test.seu@meta.data=merge(test.seu@meta.data,result_main_hpca,by="CB")
rownames(test.seu@meta.data)=test.seu@meta.data$CB

p5 <- DimPlot(test.seu, reduction = "tsne", group.by = "HPCA_Main", pt.size=0.5)+theme(
  axis.line = element_blank(),
  axis.ticks = element_blank(),axis.text = element_blank()
)
p6 <- DimPlot(test.seu, reduction = "tsne", group.by = "ident",   pt.size=0.5, label = TRUE,repel = TRUE)+theme(
  axis.line = element_blank(),
  axis.ticks = element_blank(),axis.text = element_blank()
)
fig_tsne <- plot_grid(p6, p5, labels = c('ident','HPCA_Main'),rel_widths = c(2,3))
ggsave(filename = "tsne4.pdf", plot = fig_tsne, device = 'pdf', width = 36, height = 12, units = 'cm')

总体看上去也是符合最初的鉴定的

到这儿，算是把大类鉴定了，最后咱们把注释信息添加上去

B_cell=c(0)
Plasma_cell=c(4)
T_cell=c(1,2)
NK_cell=c(5)
Unknown=c(13)
Epithelial=c(3, 6, 7, 8, 10, 11, 12)
Endothelial=c(14)
Fibroblasts=c(9)
Doublet=c(15)

current.cluster.ids <- c(B_cell,
                         Plasma_cell,
                         T_cell,
                         NK_cell,
                         Unknown,
                         Epithelial,
                         Endothelial,
                         Fibroblasts,
                         Doublet)

new.cluster.ids <- c(rep("B_cell",length(B_cell)),
                     rep("Plasma_cell",length(Plasma_cell)),
                     rep("T_cell",length(T_cell)),
                     rep("NK_cell",length(NK_cell)),
                     rep("Unknown",length(Unknown)),
                     rep("Epithelial",length(Epithelial)),
                     rep("Endothelial",length(Endothelial)),
                     rep("Fibroblasts",length(Fibroblasts)),
                     rep("Doublet",length(Doublet))
)

test.seu@meta.data$celltype <- plyr::mapvalues(x = as.integer(as.character(test.seu@meta.data$seurat_clusters)), from = current.cluster.ids, to = new.cluster.ids)

新的test.seu@meta.data是这样的：

统计一下各类细胞的数目

> table(test.seu@meta.data$celltype)

     B_cell     Doublet Endothelial  Epithelial Fibroblasts     NK_cell Plasma_cell 
       1188          22          27        2023         205         441         638 
     T_cell     Unknown 
       2167          35

咱们把Doublet和Unknown去掉，画最后的tsne图

plotCB=as.data.frame(test.seu@meta.data%>%filter(seurat_clusters!="13" &
                                                seurat_clusters!="15"))[,"CB"]
DimPlot(test.seu, reduction = "tsne", group.by = "celltype", pt.size=0.5,cells = plotCB)
ggsave(filename = "tsne5.pdf", width = 15, height = 12, units = 'cm')
saveRDS(test.seu,file = "test.seu.rds") #保存test.seu对象，下次能够直接调用，不用重复以上步骤

最后分享一个找差别基因的小技巧，有时候，咱们但愿对本身定义的类群找差别基因，能够用下面这种方法

Idents(test.seu)="celltype"
markers2 <- FindAllMarkers(test.seu, logfc.threshold = 0.25, min.pct = 0.1, 
                          only.pos = TRUE)

因水平有限，有错误的地方，欢迎批评指正！