单细胞的细胞类型注释包CAMML
大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~
就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~
这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我!
起初找到CAMML包，是因为它用了GSE72056这个数据。此数据出自2016年的一篇science文章，Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq.是smartseq数据，黑色素瘤的，4000多个细胞。
我查找资料时，发现CAMML这个包用它做了示例数据，这是一个根据细胞所对应的基因和每个基因的权重，注释细胞类型的包。出于好奇，跑了跑示例代码。
CAMML函数需要2个参数，一个是seurat对象，一个是基因集数据框。
seurat: a Seurat Object of the scRNA-seq data, gene.set.df: a data frame with a row for each gene and the following required columns: “ensemb.id” and “cell.type” and optional columns of “gene.weight” and “gene.symbol”.
这里所说的gene.weight其实就是差异分析的logFC。
1.获取Seurat对象
1.1整理数据
library(CAMML)
library(Seurat)
library(edgeR)
library(org.Hs.eg.db)
library(msigdbr)
#access data1
ground_truth <- data.table::fread("GSE72056_melanoma_single_cell_revised_v2.txt",data.table = F)
#前三行是这个数据细胞的一些信息
cells <- ground_truth[1:3,]
cells[,1:4]
## Cell
## 1 tumor
## 2 malignant(1=no,2=yes,0=unresolved)
## 3 non-malignant cell type (1=T,2=B,3=Macro.4=Endo.,5=CAF;6=NK)
## Cy72_CD45_H02_S758_comb CY58_1_CD45_B02_S974_comb Cy71_CD45_D08_S524_comb
## 1 72 58 71
## 2 1 1 2
## 3 2 1 0
cells <- unlist(matrix(cells[3,-1]));table(cells)
## cells
## 0 1 2 3 4 5 6
## 1758 2068 515 126 65 61 52
#remove cell classifications from count matrix
gt <- ground_truth[-(1:3),]
#remove duplications
gt <- gt[!duplicated(gt$Cell),cells!=0]
#reset rownames
rownames(gt) <- gt[,1]
gt <- gt[,-1]
boxplot(gt[,1:10])
注意看箱线图的取值范围？GEO页面上说这个数据是log2(tmp/10+1),不是原始count，但是帮助文档里是没做什么处理，直接用的。
为啥除以10，这么解释的： TPM values were divided by 10 since we estimate the complexity of our single cell libraries to be on the order of 100,000 transcripts and would like to avoid counting each transcript ~10 times, as would be the case with TPM, which may inflate the difference between the expression level of a gene in cells in which the gene is detected and those in which it is not detected.
Seurat github上对TMP数据的说明是：
If you have TPM data, you can simply manually log transform the gene expression matrix in the object@data slot before scaling the data. You have to replace your object@data slot with the desired gene expression matrix as follows: pbmc@data = log(x = norm + 1))
对于logTPM数据，是不做normalize这一步的，参考：https://mp.weixin.qq.com/s/LuM9EimBTGK0eDGdhv-LxA。
因此我把这个文档里的normalize步骤去掉了，后面的结果也和原文档里不一样了，但它本来就是log的数据，在normalize的时候再log一次，可就不对了，所以不一样就不一样吧。
1.2 seurat流程
#Seurat pipeline
gse72056 <- CreateSeuratObject(gt, project = "gse72056",
min.cells=100, min.features=500,
num.var.features=2000)
#data filtering and normalization
gse72056[["percent.mt"]] <- PercentageFeatureSet(gse72056, pattern = "^MT-")
gse72056 <- subset(gse72056, subset = percent.mt < 5)
#gse72056 <- NormalizeData(gse72056)
gse72056 <- FindVariableFeatures(gse72056, selection.method = "vst", nfeatures = 2000)
# Data clustering
gse72056 <- ScaleData(gse72056)
gse72056 <- RunPCA(gse72056)
gse72056 <- FindNeighbors(gse72056, dims = 1:30)
gse72056 <- FindClusters(gse72056, resolution = 0.25)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
##
## Number of nodes: 2887
## Number of edges: 99233
##
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.9330
## Number of communities: 11
## Elapsed time: 0 seconds
gse72056 <- RunUMAP(gse72056, dims = 1:30)
UMAPPlot(gse72056)
气氛到这了，我觉得得弄个singleR泡泡。
再跑一下singleR
# 注释
library(celldex)
library(SingleR)
#ref <- celldex::HumanPrimaryCellAtlasData()
#因为下载太慢，使用了提前存好的本地数据
ref <- get(load("../single_ref/ref_Human_all.RData"))
library(BiocParallel)
sce.all = gse72056
pred.scRNA <- SingleR(test = sce.all@assays$RNA@data,
ref = ref,
labels = ref$label.main,
clusters = sce.all@active.ident)
pred.scRNA$pruned.labels
## [1] "T_cells" "T_cells" "B_cell"
## [4] "Epithelial_cells" "Neurons" "Tissue_stem_cells"
## [7] "Monocyte" "T_cells" "Endothelial_cells"
## [10] "Neurons" "B_cell"
plotScoreHeatmap(pred.scRNA, clusters=pred.scRNA@rownames, fontsize.row = 9,show_colnames = T)
new.cluster.ids <- pred.scRNA$pruned.labels
names(new.cluster.ids) <- levels(sce.all)
levels(sce.all)
## [1] "0" "1" "2" "3" "4" "5" "6" "7" "8" "9" "10"
sce.all <- RenameIdents(sce.all,new.cluster.ids)
levels(sce.all)
## [1] "T_cells" "B_cell" "Epithelial_cells"
## [4] "Neurons" "Tissue_stem_cells" "Monocyte"
## [7] "Endothelial_cells"
p = UMAPPlot(object = sce.all, pt.size = 0.5, label = TRUE)
p
2.构建基因集数据
reference <- get(load("../single_ref/ref_Human_all.RData"))
#set labels
labs <- unique(reference@colData$label.main)
labs <- sort(labs)
#isolate labels for this dataset
labs <- labs[c(2,9,12,35,28,20)]
labs <- sort(labs)
labs
## [1] "B_cell" "Endothelial_cells" "Fibroblasts"
## [4] "Macrophage" "NK_cell" "T_cells"
#define relevant columns and counts
colcount <- reference@assays@data$logcounts
counts <- reference@colData$label.main
colcount <- colcount[,counts %in% labs]
dim(colcount)
## [1] 19363 263
colcount[1:4,1:4]
## GSM1209554_HH1763_UI33plus2_201004 GSM1209555_HH1778_u133plus2_211004
## A1BG 7.8836 6.6865
## A1BG-AS1 6.2969 6.8955
## A1CF 4.4967 4.1635
## A2M 5.8868 5.1582
## GSM1209556_HH1786_U133plus2_091104 GSM1209557_HH1791_u133plus2_251104
## A1BG 6.7626 7.2963
## A1BG-AS1 7.3084 6.4962
## A1CF 4.5236 4.2356
## A2M 4.8328 5.8057
counts <- counts[counts %in% labs]
table(counts)
## counts
## B_cell Endothelial_cells Fibroblasts Macrophage
## 26 64 10 90
## NK_cell T_cells
## 5 68
2.1 edgeR差异分析获取logFC
是用上面的colcount做差异分析，不是用GSE72056这个数据做哦。
f = "deg.rdata"
if(!file.exists(f)){
v <- data.frame()
for (i in 1:length(labs)){
#i = 1
d <- DGEList(counts=exp(colcount), group = ifelse(counts == labs[i],1,0))
d <- calcNormFactors(d)
d1 <- estimateCommonDisp(d, verbose=T)
d1 <- estimateTagwiseDisp(d1)
et12 <- exactTest(d1, pair = c(1,2))
gp <- et12$table
gp <- gp[order(gp$logFC, decreasing = T),]
#save gene symbols
r <- rownames(gp[gp$logFC>5,])
#save log fc
gw <- (gp[gp$logFC>5,1])
v <- rbind(v,cbind(rep(labs[i], length(r)), r,gw))
}
save(v,file = f)
}
## Disp = 0.68524 , BCV = 0.8278
## Disp = 0.60415 , BCV = 0.7773
## Disp = 0.70597 , BCV = 0.8402
## Disp = 0.57566 , BCV = 0.7587
## Disp = 0.74498 , BCV = 0.8631
## Disp = 0.61096 , BCV = 0.7816
load(f)
head(v)
## V1 r gw
## 1 B_cell RGS13 11.1408222922782
## 2 B_cell IGHV3-21 10.4089458471889
## 3 B_cell TNFRSF17 10.2031186351758
## 4 B_cell IGLV3-25 9.77139307126979
## 5 B_cell IGHV3-23 9.75375179084989
## 6 B_cell TCL1A 9.69362790838796
V1是细胞类型，r是每种细胞对应的基因，gw是logFC，这里只取了logFC>5的基因。
2.2 从msigdbr中找每种细胞的marker基因
不得不吐槽一下这段代码，太。。。了，不仅每句都多写个没必要的which，还都加上个c，示例代码也免不了画蛇添足。我懒得改啦。
#incorporate and intersect with C8
x <- c()
library(msigdbr)
m <- msigdbr(category = "C8")
#B-cell genes
r <- c(m$gene_symbol[which(m$gs_name == "HAY_BONE_MARROW_FOLLICULAR_B_CELL")])
r <- intersect(r, v[which(v[,1] == "B_cell"),2])
x <- rbind(x,cbind(rep("B_cell", length(r)), r))
#T-cell genes
r <- c(m$gene_symbol[which(m$gs_name == "HAY_BONE_MARROW_CD8_T_CELL")])
r <- intersect(r, v[which(v[,1] == "T_cells"),2])
x <- rbind(x,cbind(rep("T_cells", length(r)), r))
r <- c(m$gene_symbol[which(m$gs_name == "HAY_BONE_MARROW_NAIVE_T_CELL")])
r <- intersect(r, v[which(v[,1] == "T_cells"),2])
x <- rbind(x,cbind(rep("T_cells", length(r)), r))
#NK cells
r <- c(m$gene_symbol[which(m$gs_name == "HAY_BONE_MARROW_NK_CELLS")])
r <- intersect(r, v[which(v[,1] == "NK_cell"),2])
x <- rbind(x,cbind(rep("NK_cell", length(r)), r))
#Macrophages
r <- c(m$gene_symbol[which(m$gs_name == "HAY_BONE_MARROW_MONOCYTE")])
r <- intersect(r, v[which(v[,1] == "Macrophage"),2])
x <- rbind(x,cbind(rep("Macrophage", length(r)), r))
#Fibroblasts
r <- c(m$gene_symbol[which(m$gs_name == "CUI_DEVELOPING_HEART_C3_FIBROBLAST_LIKE_CELL")])
r <- intersect(r, v[which(v[,1] == "Fibroblasts"),2])
x <- rbind(x,cbind(rep("Fibroblasts", length(r)), r))
#Endothelial
r <- c(m$gene_symbol[which(m$gs_name == "CUI_DEVELOPING_HEART_C4_ENDOTHELIAL_CELL")])
r <- intersect(r, v[which(v[,1] == "Endothelial_cells"),2])
x <- rbind(x,cbind(rep("Endothelial_cells", length(r)), r))
#merge C8 and DE data
v <- data.frame(v)
df <- data.frame(x)
df <- merge(df, v, by = c("r","V1"),all.x = T)
2.3 ID转换
要把symbol转换为ensembl id，因为df参数必须有ensemblid列。原文档里的id转换十分啰嗦，我直接改用了bitr。
#convert gene symbols to Ensembl IDs
library(clusterProfiler)
library(org.Hs.eg.db)
s2e <- bitr(df$r,
fromType = "SYMBOL",
toType = "ENSEMBL",
OrgDb = org.Hs.eg.db)
df <- dplyr::inner_join(df,s2e,by=c("r"="SYMBOL"))
colnames(df) = c("gene.symbol","cell.type","gene.weight","ensembl.id")
head(df)
## gene.symbol cell.type gene.weight ensembl.id
## 1 A2M-AS1 T_cells 5.14903833976827 ENSG00000245105
## 2 ABCB4 B_cell 5.16248789428773 ENSG00000005471
## 3 ADAM28 B_cell 5.18531356835885 ENSG00000042980
## 4 ADAMTS2 Fibroblasts 5.23313468195672 ENSG00000087116
## 5 ADAMTS2 Fibroblasts 5.23313468195672 ENSG00000283802
## 6 ADAMTS9 Endothelial_cells 5.62309307256654 ENSG00000163638
终于准备好了这个df数据框！
3.CAMML 做细胞评分
前面的一大通都是在准备数据，核心代码只有下面这一点点：
gse72056 <- CAMML(gse72056,df)
#visualize results
gse72056@assays$CAMML@data[1:5,1:5]
## 5 x 5 sparse Matrix of class "dgCMatrix"
## Cy72_CD45_H02_S758_comb Cy81_Bulk_CD45_B10_S118_comb
## B-cell 0.983638495 0.53591861
## Endothelial-cells 0.079756170 0.99994960
## Fibroblasts 0.003692128 0.99927062
## Macrophage 0.061343977 0.05977083
## NK-cell . 0.35625812
## Cy72_CD45_D09_S717_comb Cy80_II_CD45_C09_S897_comb
## B-cell 0.01172034 0.2859924
## Endothelial-cells 0.05025042 0.9880788
## Fibroblasts . 0.9732322
## Macrophage 0.01625040 0.1853758
## NK-cell 0.01773781 .
## Cy74_CD45_F09_S453_comb
## B-cell 1.634602e-05
## Endothelial-cells 2.243491e-02
## Fibroblasts .
## Macrophage 8.471969e-03
## NK-cell 9.931903e-01
CAMML.results <- GetCAMMLLabels(gse72056,labels = "top10p")
#look at the number of cell types called for those > 90% of the max score
CAMML.results[1:4]
## [[1]]
## dif.inde...
## B-cell 0.9836385
##
## [[2]]
## dif.inde...
## Endothelial-cells 0.9999496
## Fibroblasts 0.9992706
##
## [[3]]
## dif.inde...
## T-cells 0.553577
##
## [[4]]
## dif.inde...
## Endothelial-cells 0.9880788
## Fibroblasts 0.9732322
对于GetCAMMLLabels这个函数的labels参数的解释：
One of the following: “top1”, “top2”, “top10p”, or “top2xmean”. “top1” will return the single-label for the top-scoring cell type for each cell. “top2” will return the labels for the two top-scoring cell types for each cell. “top10p” will return the top scoring cell type and any other cell types with scores within 10% of the top score for each cell. “top2xmean” will return any cell types with scores two times the average of all cell type scores for each cell.
sizetest <- c()
for (i in 1:length(CAMML.results)){
sizetest[i] <- nrow(CAMML.results[[i]])
}
table(sizetest)
## sizetest
## 1 2 3 4
## 1675 1035 167 10
table的结果，表示超过最高分90%以上的细胞类型有多少种。
#visualize how cell number relates to transitioning states in the data
gse72056$cellnum <- sizetest
UMAPPlot(gse72056, group.by = "cellnum")
As can be seen, cells in transition areas between cell types tend to have more cell type classifications
如果把labels参数写成top1，可以只看最高分的细胞。
CAMML.results2 <- GetCAMMLLabels(gse72056,labels = "top1")
table(unlist(CAMML.results2))
##
## B-cell Endothelial-cells Fibroblasts Macrophage
## 504 334 380 222
## NK-cell T-cells
## 502 945
table(cells)
## cells
## 0 1 2 3 4 5 6
## 1758 2068 515 126 65 61 52
gse72056$camml_label = unlist(CAMML.results2)
p1 = UMAPPlot(gse72056, group.by = "camml_label")
#这里是基于数据前三行自带的细胞类型来替换的。用于展示比较的结果
cells[cells == 1] <- "T"
cells[cells == 2] <- "B"
cells[cells == 3] <- "Macro"
cells[cells == 4] <- "Endo"
cells[cells == 5] <- "Fibro"
cells[cells == 6] <- "NK"
gse72056$cell_type <- cells[cells!=0]
p2 = UMAPPlot(gse72056, group.by = "cell_type")
p1 + p2
把他们和singleR的注释结果放在一起看看
p1+p2+p
4.关于CAMML包的一点理解
基因集数据是每次换了数据要自行构建的，也就是需要预先知道是哪几类细胞，再去找它们对应的基因，以及用于差异分析的对应的样本count数据。
示例数据里面提供了原作者给的注释结果，如果是没给，那应该是要从原文里找，或者是singleR注释、手动注释等，先得出细胞类型。
这个文档的核心结论是：处于过渡位置的那些细胞通常会有多种可能的细胞类型。
如果因为代码看不懂，而跟不上正文的节奏，可以来找我，系统学习。以下课程都是循环开课。下一期的时间，点进去咨询微信咯 生信零基础入门学习小组 生信入门班（四周线上直播课） 数据挖掘班（医生/医学生首选）