ixxmu / mp_duty

抓取网络文章到github issues保存
https://archives.duty-machine.now.sh/
116 stars 30 forks source link

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

Closed ixxmu closed 6 days ago

ixxmu commented 6 days ago

https://mp.weixin.qq.com/s/avHS4ojBOkNDMwGFX0jdqQ

ixxmu commented 6 days ago

每个单细胞亚群取子集后继续降维聚类分群标准操作(以b细胞为例) by 生信技能树

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

值得注意是,我把文献里面的髓系免疫细胞区分成为巨噬细胞和中性粒细胞,还有mast细胞。然后就是把b细胞是可以区分成为浆细胞和普通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

可以看到文献给出来的b细胞是7708个,然后我给出来的是6309个b细胞以及1245个浆细胞,数量上是差不多的!

我给出来的第一层次降维聚类分群后的单细胞亚群命名后在各个样品的数量分布如下所示:

单细胞亚群命名后在各个样品的数量分布

这个时候b细胞和浆细胞因为在我们的第一层次降维聚类分群的umap图里面确实是泾渭分明,所以我就把它们拆开命名了。

以b细胞为例举例说明取子集后继续降维聚类分群标准操作

通常我们拿到了肿瘤相关的单细胞转录组的表达量矩阵后的第一层次降维聚类分群通常是:

  • immune (CD45+,PTPRC),
  • epithelial/cancer (EpCAM+,EPCAM),
  • stromal (CD10+,MME,fibro or CD31+,PECAM1,endo)

参考我前面介绍过 CNS图表复现08—肿瘤单细胞数据第一次分群通用规则,这3大单细胞亚群构成了肿瘤免疫微环境的复杂。绝大部分文章都是抓住免疫细胞亚群进行细分,包括淋巴系(T,B,NK细胞)和髓系(单核,树突,巨噬,粒细胞)的两大类作为第二次细分亚群。但是也有不少文章是抓住stromal 里面的 fibro 和endo进行细分,并且编造生物学故事的。

前面我们已经介绍了心肝脾肺肾等多个器官的上皮细胞的细分亚群, 以及免疫细胞里面的髓系和B细胞细分亚群:

既然是前面的第一层次降维聚类分群后有   2-harmony/sce.all_int.rds 文件,以及对应的 phe.Rdata 注释信息。工作文件夹里面新建了  sub-cluster/sub-epi-inferCNV 这样的文件夹结构,然后重新开始新的r项目,所以需要如下所示的代码提取b细胞:

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() 

可以看到,这个时候取了子集后的单细胞转录组数据对象,仍然是进行了harmony整合哦。因为harmony操作针对的是pca分析结果,而我们取子集后相当于是一个从零开始的表达量矩阵了没有做过pca也没有任何降维聚类分群操作。也就是说,harmony操作并不会改变表达量矩阵本身。

最后仍然是手动命名,前面我们提到了:上皮细胞里面混入了淋巴系和髓系免疫细胞呢,其实做每个单细胞亚群的子集同样的操作都会有其它干扰亚群的混入。大家可以详细的阅读 这个胃癌单细胞数据集GSE163558对应的文献“Revealing the transcriptional heterogeneity of organ-specific metastasis in human gastric cancer using single-cell RNA Sequencing”,就可以看到b细胞里面其实是有t细胞的混入 :

b细胞里面其实是有t细胞的混入

大家处理这个数据集,也可以很容易的复现出来b细胞里面其实是有t细胞的混入,而且甚至是有髓系哦 :

甚至是有髓系哦

这个时候还需要B细胞细分亚群里面的基因列表来区分naive和memory这两种b细胞哦,在免疫学中,B细胞是适应性免疫系统的关键组成部分,负责产生抗体来识别和中和外来病原体。以下是关于B细胞的几个不同阶段和类型的描述,以及它们的区别:

  1. Naive B细胞(未经激活的B细胞)
  • 这些是尚未遇到其特定抗原的B细胞。它们在骨髓中成熟,并在没有遇到病原体的情况下循环于血液中。
  • Memory B细胞(记忆B细胞)
    • 记忆B细胞是在先前的免疫反应中遇到特定抗原后形成的。它们具有长期存活的能力,并能够快速响应再次出现的相同抗原,提供更快速和有效的免疫反应。
  • Germinal Center (GC) B细胞
    • GC B细胞是在淋巴器官的生发中心经历快速增殖和突变的B细胞。它们参与体细胞高频突变(somatic hypermutation, SHM)和抗体类别转换(class-switch recombination, CSR),以产生更高亲和力的抗体。
  • Plasma Cells(浆细胞)
    • 浆细胞是B细胞激活和分化的最终阶段,主要功能是大量产生和分泌抗体。它们可以在骨髓中长期存活,并作为长期的抗体来源。

    区别

    • 成熟阶段:Naive B细胞是成熟的早期阶段,而GC B细胞和记忆B细胞是成熟过程中的更晚阶段。
    • 功能:Naive B细胞尚未发挥功能,GC B细胞参与抗体的亲和力成熟,记忆B细胞提供长期免疫记忆,浆细胞负责抗体的产生。
    • 持久性:记忆B细胞和浆细胞在体内有较长的寿命,而Naive B细胞和GC B细胞的寿命相对较短。
    • 分布:Naive B细胞主要在血液和脾脏中循环,GC B细胞主要在淋巴结的生发中心,记忆B细胞和浆细胞可以在骨髓、脾脏和其他组织中找到。
    • 抗原经验:Naive B细胞没有遇到抗原的经验,而记忆B细胞和浆细胞是由先前接触过抗原的B细胞发展而来。

    了解这些B细胞类型的区别对于研究免疫反应、疫苗开发和某些免疫相关疾病的治疗非常重要。

    根据文献积累,我们大概有下面的这些基因来区分不同的b细胞,如下所示:

    区分naive和memory这两种b细胞

    我给出来了对应关系,然后就做出来了比较合理的单细胞亚群注释 ;

      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 

    这些基因就可以去做go和kegg数据库的注释;

    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('../')

    同理,接下来就可以看不同亚群的比例,针对这个细分后的结果进行拟时序分析,转录因子分析等等。这些分析都是在任意单细胞亚群是通用的代码。

    写在文末

    如果你也想做单细胞转录组数据分析,最好是有自己的计算机资源哦,比如我们的2024的共享服务器交个朋友福利价仍然是800,而且还需要有基本的生物信息学基础,也可以看看我们的生物信息学马拉松授课(买一得五) ,你的生物信息学入门课。而且下周六日我们在广州线上授课哦:千呼万唤,让我们广州线下约起