每个单细胞亚群取子集后继续降维聚类分群标准操作（以b细胞为例）

> print(cell_groups)
          Group Count               Genes
1    Epithelial  1743 EPCAM, KRT19, CLDN4
2       Stromal  1288 PECAM1, CLO1A2, VWF
3 Proliferative  1089  MKI67, STMN1, PCNA
4             T 24448     CD3D, CD3E, CD2
5             B  7708 CD79A, IGHG1, MS4A1
6            NK  1173  KLRD1, GNLY, KLRF1
7       Myeloid  5519  CSF1R, CSF3R, CD68

rm(list=ls())
options(stringsAsFactors = F) 
source('../../scRNA_scripts/lib.R')
set.seed(123456789) 
###### step1:导入数据 ######    
sce.all.int = readRDS('../../2-harmony/sce.all_int.rds')
colnames(sce.all.int@meta.data) 
load('../../phe.Rdata')
table(phe$celltype) 
# 因为前面的b细胞和浆细胞被我拆开命名，所以这个时候细分子集需要把两个都提取到
choose_cellIds=rownames(phe)[phe$celltype %in% c('Bcells','plasma') ]
sce.all=sce.all.int[,choose_cellIds]
sce.all=CreateSeuratObject(
  counts = sce.all@assays$RNA$counts,
  meta.data = sce.all@meta.data
) 

###### step2:QC质控 ######
dir.create("./1-QC")
setwd("./1-QC")
# 如果过滤的太狠，就需要去修改这个过滤代码
source('../../../scRNA_scripts/qc.R')
sce.all.filt = basic_qc(sce.all)
print(dim(sce.all))
print(dim(sce.all.filt))
setwd('../')
getwd()

sp='human' 
###### step3: harmony整合多个单细胞样品 ######

if(T){
  dir.create("2-harmony")
  getwd()
  setwd("2-harmony")
  source('../../../scRNA_scripts/harmony.R')
  # 默认 ScaleData 没有添加"nCount_RNA", "nFeature_RNA"
  # 默认的
  sce.all.int = run_harmony(sce.all.filt)
  setwd('../')
  
}
 
###### step4:  看标记基因库 ######
# 原则上分辨率是需要自己肉眼判断，取决于个人经验
# 为了省力，我们直接看  0.1和0.8即可

table(Idents(sce.all.int))
table(sce.all.int$seurat_clusters)
table(sce.all.int$RNA_snn_res.0.1) 
table(sce.all.int$RNA_snn_res.0.8) 

getwd()
dir.create('check-by-0.1')
setwd('check-by-0.1')
sel.clust = "RNA_snn_res.0.1"
sce.all.int <- SetIdent(sce.all.int, value = sel.clust)
table(sce.all.int@active.ident) 

source('../../../scRNA_scripts/check-all-markers.R')
setwd('../') 
getwd()

dir.create('check-by-0.5')
setwd('check-by-0.5')
sel.clust = "RNA_snn_res.0.5"
sce.all.int <- SetIdent(sce.all.int, value = sel.clust)
table(sce.all.int@active.ident) 
source('../../../scRNA_scripts/check-all-markers.R')
setwd('../') 
getwd()

dir.create('check-by-0.8')
setwd('check-by-0.8')
sel.clust = "RNA_snn_res.0.8"
sce.all.int <- SetIdent(sce.all.int, value = sel.clust)
table(sce.all.int@active.ident) 
source('../../../scRNA_scripts/check-all-markers.R')
setwd('../') 
getwd() 

  celltype[celltype$ClusterID %in% c( 1 ),2]='navie'
  celltype[celltype$ClusterID %in% c( 0 ),2]='memory'
  celltype[celltype$ClusterID %in% c(2),2]='plasma'  
  celltype[celltype$ClusterID %in% c(4),2]='GC'   
  celltype[celltype$ClusterID %in% c(3),2]='Tcells'      
  celltype[celltype$ClusterID %in% c(5),2]='myeloids'  

clustermemory:TNFRSF13B,LINC01781,ZBTB20,TNFSF9,IER5,WSB1
 clustermyeloids:FCGR3B,S100A9,S100A8,CXCL8,CMTM2,G0S2
 clusterplasma:MZB1,FKBP11,DERL3,PRDX4,SLAMF7,FNDC3B
 clusterTcells:IL32,CD3E,CD3D,CD2,CD3G,IL7R
 clusternavie:YBX3,IGHD,TCL1A,FCER2,ZBTB16,CLEC2B
 clusterGC:LMO2,SUGCT,SERPINA9,RAPGEF5,MYBL1,AC023590.1 

load('../check-by-celltype/qc-_marker_cosg.Rdata') 
head(marker_cosg)
 symbols_list <-  as.list(as.data.frame(apply(marker_cosg$names,2,head,100)))
 symbols_list
 source('../../com_go_kegg_ReactomePA_human.R')
#source('../com_go_kegg_ReactomePA_mice.R')
com_go_kegg_ReactomePA_human(symbols_list, pro='b' )
setwd('../')