前些天我们兼职工程师交流群的小伙伴们反馈了一个“鸡同鸭讲”的数据分析售后单子，客户强烈要求必须复现出来跟张泽民大佬的单细胞原文一模一样的降维聚类分群和生物学命名结果。
这个文章很经典，发表在2019的CELL期刊：《Landscape and Dynamics of Single Immune Cells in Hepatocellular Carcinoma》，其中 单细胞表达量矩阵文件是：
ls -lh *gz|cut -d" " -f 6-
250M 8 30 11:20 GSE140228_UMI_counts_Droplet.mtx.gz
308K 8 30 11:19 GSE140228_UMI_counts_Droplet_barcodes.tsv.gz
478K 8 30 11:19 GSE140228_UMI_counts_Droplet_cellinfo.tsv.gz
1.0M 8 30 11:19 GSE140228_UMI_counts_Droplet_genes.tsv.gz
我们可以很简单的读取上面的表达量矩阵文件以及作者给出来了各个细胞的身份信息，然后做降维聚类分群，如下所示 ：
这个数据分析过程，最大的难点，就是读取作者给出来的文件，如下所示：
library(data.table)
library(Matrix)
mtx=readMM( "inputs/GSE140228_UMI_counts_Droplet.mtx.gz" )
cl=fread( "inputs/GSE140228_UMI_counts_Droplet_barcodes.tsv.gz" ,
header = F,data.table = F )
rl=fread( "inputs/GSE140228_UMI_counts_Droplet_genes.tsv.gz" ,
header = T,data.table = F )
kp=!duplicated(rl$SYMBOL)
table(kp)
dim(mtx)
mtx=mtx[kp,]
rl=rl[kp,]
rownames(mtx)=rl$SYMBOL
colnames(mtx)=cl$V1 #paste0('c',cl$V2)
mtx[1:4,1:4]
dim(mtx)
meta=fread( "inputs/GSE140228_UMI_counts_Droplet_cellinfo.tsv.gz" ,
header = T,data.table = F )
head(meta)
rownames(meta)=meta[,1]
meta=meta[,-1]
head(meta)
identical(rownames(meta),colnames(mtx))
library(Seurat)
sce.all=CreateSeuratObject(counts = mtx ,
meta.data = meta,
min.cells = 5,
min.features = 300 )
如果我们没有读取作者自己的单细胞亚群名字文件，是不可能拿到跟文章一模一样的各个单细胞亚群结果 ：
首先呢，作者的每个单细胞亚群进行了大量的细分，比如cd4和cd8还有nk都是各自有好多个细分亚群，这些信息在我们自己的UMAP里面肯定是会糅合在一起的，因为哪怕是cd4和cd8还有nk都是有互相掺杂的，如下所示：
因为cd4和cd8还有nk在第一层次降维聚类分群的UMAP本来就是不应该是泾渭分明，就算是我们强行的想把细分亚群添加在umap上面，也是看起来就“杂乱无章”：
但是客户根本就不理解各个结果复现的重要性
毕竟张泽民大佬的权威性太高，所以初次接触单细胞的小伙伴们肯定是想朝他看齐。
