数据分析能力确实是跟能发表的期刊质量有关

学员在交流群问到了如何比较不同单细胞数据挖掘文章，其中一个新鲜出炉的单细胞层面的泛癌文章，是英国的Kevin Litchfield课题组，发表在期刊《cancer cell 》，另外一个是张泽民老师发表在CELL杂志的同样的泛癌b细胞单细胞转录组数据探索性课题。

文献标题是：《A pan-cancer single-cell RNA-seq atlas of intratumoral B cells》，文献里面 纳入了15 studies spanning 7 cancer types ，最后是 126,101 TIL-Bs across 10 clusters ，但是图可以说是丑爆了：

虽然说确实是很丑，但是还是蛮符合规则的，因为里面的plasma和其它b细胞确实是还算是泾渭分明，而且增殖的细胞会成为两个完全不同的b细胞大亚群的桥梁，另外就是gc的b细胞也是可以跟增殖的区分开。但是跟张泽民老师发表在CELL杂志的数据分析比起来就有“云泥之别”，如下所示：

值得注意是，张泽民老师其实分析的公共数据集更多一点：

• Our newly generated pan-cancer scRNA-seq dataset was assembled with 54 additional published studies to elucidate the pan-cancer characteristics of TIBs
• containing 511,847 B cells derived from 948 samples of 649 patients across 19 major can- cer types, primarily covering tumor tissues, adjacent non-tumor tissues (ANTs), and peripheral blood

使用singleR对两个单细胞转录组数据集进行细胞亚群注释信息的迁移

首先读取张泽民老师的数据挖掘里面的b细胞部分信息进行降维聚类分群（略）后，取子集 ：

sce.all = readRDS('../2-harmony/sce.all_int.rds')
sp='human' 
load('../phe.Rdata') 
identical(rownames(phe) , colnames(sce.all))
sce.all@meta.data = phe
sel.clust = "celltype"
sce.all <- SetIdent(sce.all, value = sel.clust)
table(sce.all@active.ident) 
DimPlot(sce.all) 
colnames(sce.all@meta.data) 
sel.clust = "celltype"
sce.all <- SetIdent(sce.all, value = sel.clust)  
sce.all
table(Idents(sce.all))
sce.main.100 = subset(sce.all,downsample=100)
sce.main.100
save(sce.main.100,file = 'GSE233236-zzm-sub.sce.Rdata')

然后读取英国的Kevin Litchfield课题组的数据挖掘里面的b细胞部分信息进行降维聚类分群后（略），使用singleR对两个单细胞转录组数据集进行细胞亚群注释信息的迁移：

ource('scRNA_scripts/lib.R')
sce.all = readRDS('2-harmony/sce.all_int.rds')
sce.all  
sce.all.int= sce.all
 
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(dplyr)
library(stringr)
library(ggsci)  
library(patchwork)
load('./GSE233236-zzm-sub.sce.Rdata')  
table(Idents(sce.main.100) )
##NOTE：以前是AverageExpression
av <-AggregateExpression(sce.main.100 ,      # group.by = "celltype",
                         assays = "RNA") 
Ref = av[[1]]
head(Ref)
ref_sce=SingleCellExperiment::SingleCellExperiment(assays=list(counts=Ref))
library(scater)
ref_sce=scater::logNormCounts(ref_sce)
library(SingleCellExperiment)
logcounts(ref_sce)[1:4,1:4]
colData(ref_sce)$Type=colnames(Ref)
table(colnames(Ref))
ref_sce

testdata <- GetAssayData(sce.all, slot="data")
testdata[1:4,1:4]
dim(testdata)
testdata <-sce.all@assays$RNA$counts 
library(SingleR)
pred <- SingleR(test=testdata, ref=ref_sce, 
                labels=ref_sce$Type
)
as.data.frame(table(pred$labels))
head(pred) 
labels=pred$labels
table(labels)  
save(labels,file = 'SingleR_celltype.Rdata') 

这样的话，b细胞的细分亚群的生物学注释信息勉强看得过去，起码比文章里面的UMAP要好一点：

不过，这样的话，就很难复现出来文献里面的分群啦， 比如文章里面的  MT1X-high Plasma/Plasmablasts   ，就是 金属离子酶亚群！

其实 不同单细胞转录组数据集的降维聚类分群其实都会有 热激蛋白的亚群，细胞增殖亚群，干扰素亚群，金属离子酶亚群，线粒体或者核糖体亚群，或者低质量亚群...  我们之前就分享过：【巨噬细胞新分类体系（放弃传统M1和M2）】，一切都是数据结果合理的挑选和解释而已。大家在自己的单细转录组数据降维聚类分群的时候也可以看到CXCL9 and SPP1的排他性，跟前面提到的TREM2联合SPP1去和FOLR2基因的排他性类似的。另外，我们通常是并不会选择提高分辨率这个手段来获取精细的单细胞亚群，而是取巨噬细胞子集后，继续进行降维聚类分群后再命名的策略。参考:  取单细胞亚群子集细分的时候一定会出现干扰亚群（所以不要惊慌）， 这样的话你就可以看到细胞亚群里面的混杂因素，而且可以手动删除到干扰因素。

                               cycle    GC memory naive plasma
  Naive B cells                    11   167   4846 10969      7
  Activated B cells                48   141  15443  6037    562
  Resting Memory B cells          138   392  20651  3901   1781
  IGKC-high Plasma/Plasmablasts   113    16    334   129   6593
  Conventional Plasma cells       438    38    346    69  19588
  Stressed Plasma cells            88    27   3733   447   6103
  Atypical Memory B cells         109    65   6707   876    278
  Proliferative B cells          5463   132    363   131     34
  GC B cells                      316  2975    522   132      1
  MT1X-high Plasma/Plasmablasts    41     7    692   179   1587

写在文末