免疫性血小板低下症患者的单细胞层面变化

生信技能树

发布于 2023-02-27 21:28:04

51300

代码可运行

运行总次数：0

代码可运行

今天是大年初一，年前的7天已经处理了60多个单细胞数据集并且发给了有需求的粉丝们。其实本来呢春节前后就14天，一天哪怕是干10个数据集，也很难完成全部的委托了，好在前面两个月培养的一百多个在线实习生，毕竟教了大家R语言，转录组，以及单细胞转录组。大家一起努力，才有可能把我们公众号活动坚持下来，详见：春节期间单细胞转录组数据分析全免费。

现在要分享的是数据集GSE196676，关于Immune thrombocytopenia (ITP), 中文翻译可以是免疫性血小板低下症，对应发表的文章也很新，Deciphering transcriptome alterations in bone marrow hematopoiesis at single-cell resolution in immune thrombocytopenia. Signal Transduct Target Ther 2022 Oct 7;7(1):347. PMID: 36202780

首先读取作者提供的表达量矩阵文件

这个数据集在：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE196676

可以看到作者提供的是一个单一的表达量矩阵文件：GSE196676_BM_HSPCs_gene_counts.csv.gz


rm(list=ls())
options(stringsAsFactors = F) 
source('scRNA_scripts/lib.R')

###### step1:导入数据 ######   
# 付费环节 800 元人民币
 
ct=fread( 'GSE196676_BM_HSPCs_gene_counts.csv.gz',data.table = F)
ct[1:4,1:4]
ct[nrow(ct),1:4]
rownames(ct)=ct[,1]
ct=ct[,-1]
ct[1:4,1:4]

sce.all =CreateSeuratObject(counts =  ct ) 
as.data.frame(sce.all@assays$RNA@counts[1:10, 1:2])
head(sce.all@meta.data, 10)
table(sce.all$orig.ident)
rm(ct)

library(stringr)
phe=str_split( colnames(sce.all),'[-_]',simplify = T)
head(phe) 
table(phe[,2]) 
sce.all$orig.ident=paste0('p',phe[,2]) 
table(phe[,1]) 
sce.all$group= phe[,1]

如果使用上面的代码，看到的是4个病人的两个分组，可能会误解为是同一个病人两次不同疾病状态取样：

> table(sce.all$group,sce.all$orig.ident)
     
        p1   p2   p3   p4
  HC  6971 6384 8289 6161
  ITP 5451 7844 8866 6346

对应的是 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE196676 页面信息：

GSM5898060 BM, CD34+HSPCs, ITP, 1, 10X
GSM5898061 BM, CD34+HSPCs, ITP, 2, 10X
GSM5898062 BM, CD34+HSPCs, ITP, 3, 10X
GSM5898063 BM, CD34+HSPCs, ITP, 4, 10X
GSM5898064 BM, CD34+HSPCs, HC, 1, 10X
GSM5898065 BM, CD34+HSPCs, HC, 2, 10X
GSM5898066 BM, CD34+HSPCs, HC, 3, 10X
GSM5898067 BM, CD34+HSPCs, HC, 4, 10X

另外，页面写的是 bone marrow (BM) CD34+ hematopoietic stem and progenitor cells (HSPCs) from four newly diagnosed treatment-naïve ITP patients and four healthy donors

所以应该是8个病人，但是我读入后做了一个操作是：sce.all$orig.ident=paste0('p',phe[,2])

其实是不太合理的，真正的sce.all$orig.ident 应该是phe的第一列和第二列的合并。不过问题不大。

然后质量控制并且降维聚类分群即可

质量控制其实每个数据集不一样的，取决于单细胞转录组来源的技术，特殊情况下需要个性化的修改 scRNA_scripts/qc.R 里面的代码，尤其注意的是本文研究的就是Immune thrombocytopenia (ITP), 所以这个时候千万不能过滤血小板相关的基因列表啦：

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

但是多样品整合，以及后续的降维聚类分群这个代码在任意单细胞转录组数据集都是大概率不需要修改的：

sp='human'

###### step3: harmony整合多个单细胞样品 ######
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.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()

然后就是合理的生物学命名了

前面强调了研究团队做单细胞转录组的时候，筛选了免疫细胞，所以降维聚类分群后主要是免疫细胞亚群进行细分，包括淋巴系（T,B,NK细胞）和髓系（单核，树突，巨噬，粒细胞）的两大类作为第二次细分亚群，但是我发现里面仍然是有少量的上皮细胞，内皮细胞，成纤维细胞。

这样的话，就跟我给大家准备的基因列表主要是针对肿瘤单细胞的第一层次降维聚类分群，是：

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

参考我前面介绍过 CNS图表复现08—肿瘤单细胞数据第一次分群通用规则，这3大单细胞亚群构成了肿瘤免疫微环境的复杂。

如果使用我们的基因会发现基本上不太可能给出合理的生物学名字：

不太可能给出合理的生物学名字

起码都是免疫细胞没有问题，成熟的t细胞和b细胞很少，髓系就更缺乏了，勉强可以看到一点点肥大细胞。这就需要看文章了：

造血干细胞(HSC)
多能祖细胞(MPP)
粒细胞-巨噬细胞祖细胞(GMP)
嗜中性粒细胞祖细胞(NeuP)
单核细胞-树突细胞祖细胞(MDP)
嗜酸性粒细胞-嗜碱性粒细胞-肥大细胞祖细胞(EBMP)
常见淋巴祖细胞(CLP)
三个前 B 细胞群(preB1, preB2和 preB3)
自然杀伤/T 细胞祖细胞(NK/Tp)
巨核细胞-红细胞祖细胞(MEP)
两个巨核细胞祖细胞亚群(MkP1和 MkP2)
红细胞祖细胞(EryP)

而且一般来说，文章也会给出来具体的各个单细胞亚群的特异性基因，如下所示：

各个单细胞亚群的特异性基因

所以只需要把这些基因在我们的前面的降维聚类分群里面可视化出来，就可以肉眼给这先单细胞亚群生物学名字了。


genes_to_check =list(
  HSC=c("AVP","KLF2","CRHBP","HOPX","BST2","SPINK2"),
  MPP=c("C1QTNF4","SMIM24","KIAA0125","SELL"),
  GMP=c("MPO","ELANE","AZU1","PRTN3","CALR","CFD","MGST1","PRSS57"),
  NeuP=c("LYZ","LGALS1","CTSG","ANXA2"),
  MDP=c("IRF7","IRF8","JCHAIN","PLD4","LILRA4","GZMB","CCDC50","SPIB","SCT","UGCG"),
  EBMP=c("CLC","TPSAB2","MS4A3","TPSAB1","MS4A2","LMO4","HDC"),
  CLP=c("IL7R","ADA","TRBC2","SATB1"),
  preB1=c("HIST1H4C","HMGB2","NUSAP1","MK167","STMN1","CCND3","HMGN2"),
  preB2_3=c("CD9","CD79B","VPREB3","C10orf10","VPREB1","DNTT","ARPP21","AKAP12","NEIL","CD79A","MZB1","MDM4"),
  #preB3=c("GATA1","GATA2","TESPA1","CTNNBL1"),
  NKT_Tp=c("IL32","KLRB1","TRAC","CD3E","TRBC1","CD3D","CCL5","CD7","NKG7"),
  MEP=c("APOC1","HBD","CNR1P1","KLF1","BLVRB","LINC00152","TMEM14C","ITGA2B","MYC"),
  MKp1=c("GATA2","FCER1A","PBX1","CYTL1","SLC40A1"),
  MKp2=c("PPBP","PF4","NRGN","PLEK","RGS18","SDPR","PDLIM1"),
  Eryp=c("HBB","HBA2","HBA1","AHSP","CA1","HBG2","FAM178B","PRDX2")
) 

th=theme(axis.text.x = element_text(angle = 45, 
                                    vjust = 0.5, hjust=0.5)) 
genes_to_check=lapply(genes_to_check,str_to_upper) #列表的大小写转换
p <- DotPlot(sce.all.int, features = genes_to_check,
             group.by = 'RNA_snn_res.0.1',
             assay='RNA'  )   +th

p
ggsave('check_blood_markers_by_RNA_snn_res.0.1.pdf',
       height = 6,width = 20)

p <- DotPlot(sce.all.int, features = genes_to_check,
             group.by = 'RNA_snn_res.0.8',
             assay='RNA'  )   +th

p
ggsave('check_blood_markers_by_RNA_snn_res.0.8.pdf',
       height = 12,width = 20)

如下所示，这些亚群确实是可以通过作者提供的特异性基因列表进行合理的区分和命名：

合理的区分和命名

就是工作量有点大，需要人为的把这些亚群合并。

单细胞春节福利继续

如果你有自己的单细胞转录组数据需要处理，可以直接参加我们的活动，详见：春节期间单细胞转录组数据分析全免费，如果你感兴趣我们处理过的数据集，可以直接填表留下联系方式即可，我们会在元宵节之前统一发放全部的代码和数据给大家的：

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2023-01-22，如有侵权请联系 cloudcommunity@tencent.com 删除

数据分析

本文分享自生信技能树微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

数据分析

登录后参与评论

0 条评论

热度