上皮细胞里面混入了淋巴系和髓系免疫细胞呢

生信技能树

发布于 2024-11-21 09:44:39

1790

针对这个这个胃癌单细胞数据集GSE163558，我做了解读，详见：单细胞转录组降维聚类分群过滤基因和过滤细胞的区别 。而且前面已经是完成了降维聚类分群，在学习单细胞亚群命名的层次结构 演示了一个降维聚类分群结果，就有了 2-harmony/sce.all_int.rds 文件，以及对应的 phe.Rdata 注释信息。

在文献里面的单细胞亚群以及其对应的基因和细胞数量分别是：

> 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

也就是说，我把髓系免疫细胞区分成为巨噬细胞和中性粒细胞，还有mast细胞，而且把b细胞是可以区分成为浆细胞和普通b细胞。跟文章不一样的哦，所以我得到如下所示的亚群：

可以看到每个单细胞亚群都是泾渭分明的，而且有一个独立的细胞增殖状态的亚群它是多种不同单细胞亚群的混合物。

提取上皮细胞进行细分亚群

我在第一层次降维聚类分群的工作文件夹里面新建了 sub-cluster/sub-epi-inferCNV 这样的文件夹结构，然后重新开始新的r项目，所以需要如下所示的代码提取上皮细胞：

sce.all.int = readRDS('../../2-harmony/sce.all_int.rds')
colnames(sce.all.int@meta.data) 
load('../../phe.Rdata')
table(phe$celltype)  
sce.all=sce.all.int[,phe$celltype %in% c('epi' ) ]
sce.all=CreateSeuratObject(
  counts = sce.all@assays$RNA$counts,
  meta.data = sce.all@meta.data
)

就读取上面的两层的文件夹里面的 2-harmony/sce.all_int.rds 文件，以及对应的 phe.Rdata 注释信息。然后仍然是默认的降维聚类分群而已：

sp='human'
###### 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()
 
###### 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('../')
  
}

sce.all.tmp <- RunUMAP(sce.all.int, dims = 1:15, 
                       reduction = "pca")
DimPlot(sce.all.int,group.by = 'orig.ident') + 
  DimPlot(sce.all.tmp,group.by = 'orig.ident')


###### step4:  看标记基因库 ######
# 原则上分辨率是需要自己肉眼判断，取决于个人经验
# 为了省力，我们直接看  0.1 
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()

这个时候很明显的可以看到上皮细胞的降维聚类分群居然是出现了淋巴系和髓系免疫细胞的干扰哦，如下所示：

可以看到，cluster4是t淋巴细胞，而cluster6应该是髓系免疫细胞（可能是巨噬细胞或者单核细胞，中性粒细胞）。

值得注意是，这个是上皮细胞的子集，我仍然是使用了harmony算法去除个体差异，如果是要考虑肿瘤病人间的异质性就不推荐使用harmony哦！

进一步确定这些是淋巴系和髓系免疫细胞而不是上皮细胞

如果是从上面的标记基因来看，确实是全部的细分亚群都是有上皮细胞特异性的基因的，但是如果看每个亚群的top基因和功能，就可以确定这些是淋巴系和髓系免疫细胞而不是上皮细胞啦。

首先看看每个细分亚群的top10的基因，如下所示：

cluster0:GABRP,RARRES3,TRIM29,CD55,ITGB4,AL121761.2
 cluster1:TOP2A,CCNB1,ROMO1,OLFM4,GSTP1,TPX2
 cluster2:HIST1H2BG,AL137077.2,CSKMT,PCK1,LINC00668,SOCS1
 cluster3:MT1F,MT1M,METTL7A,CHGA,CPE,PTPRN2
 cluster4:CD2,CD52,CD3D,CCL5,HCST,CD3E
 cluster5:NPPB,TNF,CYP2E1,C6orf222,AC092069.1,LY6K
 cluster6:FCGR3B,PLEK,G0S2,OSM,CSF3R,BCL2A1

可以看到，cluster4是t淋巴细胞，而cluster6应该是中性粒细胞这样的髓系免疫细胞。基本上就可以手动给出来了这些单细胞亚群名字啦；

# 如果是手动给各个单细胞亚群命名
if(T){
  sce.all.int
  celltype=data.frame(ClusterID=0:6 ,
                      celltype= 0:6  ) 
  #定义细胞亚群        
  celltype[celltype$ClusterID %in% c( 4  ),2]='Tcells'      
  celltype[celltype$ClusterID %in% c( 6 ),2]='myeloids'  
  celltype[celltype$ClusterID %in% c( 3 ),2]='normal'  
  
  celltype[celltype$ClusterID %in% c( 0 ),2]='epi-1'
  celltype[celltype$ClusterID %in% c( 1 ),2]='epi-2' 
  celltype[celltype$ClusterID %in% c( 2 ),2]='epi-3'    
  celltype[celltype$ClusterID %in% c( 5 ),2]='epi-4'    
   
  head(celltype)
  celltype
  table(celltype$celltype)
  sce.all.int@meta.data$celltype = "NA"
  
  for(i in 1:nrow(celltype)){
    sce.all.int@meta.data[which(sce.all.int@meta.data$RNA_snn_res.0.1 == celltype$ClusterID[i]),'celltype'] <- celltype$celltype[i]}
  Idents(sce.all.int)=sce.all.int$celltype
  
  
} 
table(Idents(sce.all.int))

然后就可以看功能：

sce = sce.all.int
  sel.clust = "celltype"
  sce <- SetIdent(sce, value = sel.clust)
  table(sce@active.ident) 
  markers <- FindAllMarkers(object = sce, only.pos = TRUE, 
                            min.pct = 0.25, 
                            thresh.use = 0.25)
  #write.csv(markers,file='markers.csv')
  table(markers$cluster)
  library(dplyr)  
  library(clusterProfiler)  
  markers = markers[markers$p_val <0.01 & markers$avg_log2FC >1,]  
  table(markers$cluster)
  symbols_list = split(markers$gene,markers$cluster)
  source('../com_go_kegg_ReactomePA_human.R') 
  com_go_kegg_ReactomePA_human(symbols_list, pro='epi' )

上面的代码需要source一个外部的R代码文件哦，然后就可以愉快的做go和kegg注释并且一次性出很多图啦！很明显的可以看到确实是命名为t细胞和髓系是合理的：