【原】百萬(wàn)細(xì)胞的Seurat對(duì)象依據(jù)第一層次降維聚類分群拆分成為子集

健明 2024-12-29 發(fā)布于廣東

展開(kāi)全文

前面在百萬(wàn)細(xì)胞舍我其誰(shuí)（一晚上解決戰(zhàn)斗），我們得到了張澤民老師的2024的文章：《Spatiotemporal single-cell analysis decodes cellular dynamics underlying different responses to immunotherapy in colorectal cancer》的百萬(wàn)單細(xì)胞數(shù)據(jù)集的第一層次降維聚類分群，并且保存完了r對(duì)象文件：

7.3G 12月 26 06:11 sce.all_int.rds

但其實(shí)完全沒(méi)有必要，因?yàn)樽髡咭呀?jīng)是給出來(lái)了每個(gè)單細(xì)胞的身份，我們?yōu)槭裁醋鼋稻S聚類分群并且umap可視化，就是為了輔助我們給細(xì)胞一個(gè)身份。既然結(jié)果都給出來(lái)了，那么過(guò)程就不重要了。接下來(lái)就應(yīng)該是每個(gè)單細(xì)胞亞群進(jìn)行細(xì)致的探索。仍然是寫(xiě)一個(gè)r腳本，內(nèi)容如下所示：

library(Seurat)
sce.all.int = readRDS('2-harmony/sce.all_int.rds')
sp='human'
tmp=sce.all.int@meta.data
colnames(sce.all.int@meta.data) 
table(sce.all.int$RNA_snn_res.0.8) 
Sys.time() 
id = unique(sce.all.int$MajorCellType)
lapply(id, function(x){
  sce=sce.all.int[,sce.all.int$MajorCellType %in% x]
  sce.all=CreateSeuratObject(
    counts = sce@assays$RNA$counts,
    meta.data = sce@meta.data
  )
  save(sce.all,file = paste0(x,'.sce.all.Rdata'))
  Sys.time()
})

假如這個(gè)腳本文件名字是（tmp.R ），接下來(lái)就可以在shell里面提交它到后臺(tái)運(yùn)行：

nohup Rscript tmp.R &

這個(gè)Seurat對(duì)象本來(lái)就是接近100萬(wàn)個(gè)細(xì)胞，但是我進(jìn)行了簡(jiǎn)單的質(zhì)量控制，過(guò)濾后剩下90萬(wàn)左右，前面的降維聚類分群過(guò)程就省略了，參考在百萬(wàn)細(xì)胞舍我其誰(shuí)（一晚上解決戰(zhàn)斗），現(xiàn)在我們直接看作者的結(jié)果：

[1]  36027 975275
[1]  36027 895488

      B     Epi     ILC     Mye Stromal       T 
 161935   96179   48537   68250   40433  480154

可以看到，一多半是t細(xì)胞然后四分之一左右是b細(xì)胞，其它的都占比很少很少。

 484M 12月 26 11:58 B.sce.all.Rdata
 296M 12月 26 11:47 Epi.sce.all.Rdata
 133M 12月 26 12:01 ILC.sce.all.Rdata
 236M 12月 26 12:00 Mye.sce.all.Rdata
 157M 12月 26 11:48 Stromal.sce.all.Rdata
 1.3G 12月 26 11:55 T.sce.all.Rdata

可以看到，因?yàn)閠細(xì)胞數(shù)量最多，所以導(dǎo)出它這個(gè)子集耗時(shí)最長(zhǎng)：

[1] "2024-12-26 11:48:34 CST"
[1] "2024-12-26 11:55:43 CST"
[1] "2024-12-26 11:58:40 CST"
[1] "2024-12-26 12:00:11 CST"
[1] "2024-12-26 12:01:11 CST"

接下來(lái) 同樣的流程降維聚類分群

在百萬(wàn)細(xì)胞舍我其誰(shuí)（一晚上解決戰(zhàn)斗），我們舉例了，所有的r代碼寫(xiě)完了，就提交到后臺(tái)后等半天就ok了，全程都是代碼自動(dòng)化處理而已，不需要自己干啥子就可以拿到海量的結(jié)果！一切的難點(diǎn)，就在于不同的單細(xì)胞公共數(shù)據(jù)集的讀取，只需要成為了r編程語(yǔ)言里面的Seurat對(duì)象，就可以全流程自動(dòng)化處理~

那么，我們這個(gè)時(shí)候把前面的百萬(wàn)單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)集拆分成為了第一層次是6個(gè)大的亞群：six major cell types: T cells, B cells, innate lymphoid cells (ILCs), myeloid cells, stromal cells, and epithelial：

每個(gè)具體的亞群，都可以是當(dāng)做是一個(gè)獨(dú)立的單細(xì)胞轉(zhuǎn)錄組項(xiàng)目，很簡(jiǎn)單的load進(jìn)去就是一個(gè)r編程語(yǔ)言里面的Seurat對(duì)象，就可以全流程自動(dòng)化處理~

load('../Stromal.sce.all.Rdata')

sce.all.int$celltype=sce.all.int$SubCellType
table(sce.all.int$celltype)   
Sys.time()

if("SubCellType" %in% colnames(sce.all.int@meta.data ) ){
  
  sel.clust = "SubCellType"
  sce.all.int <- SetIdent(sce.all.int, value = sel.clust)
  table(sce.all.int@active.ident) 
  
  dir.create('check-by-SubCellType')
  setwd('check-by-SubCellType')
  source('../scRNA_scripts/check-all-markers.R')
  setwd('../') 
  getwd()
  phe=sce.all.int@meta.data
  save(phe,file = 'phe.Rdata')
  pdf('SubCellType-vs-orig.ident.pdf',width = 10)
  gplots::balloonplot(table(sce.all.int$SubCellType,sce.all.int$orig.ident))
  dev.off()
}

Sys.time()

上面的r代碼是需要自己保證準(zhǔn)確無(wú)誤的，然后寫(xiě)入到文件：step1-load-by-Seurat-v5.R，接下來(lái)就可以nohup的在后臺(tái)運(yùn)行它啦~

在我們的服務(wù)器的Linux的診斷運(yùn)行：nohup Rscript step1-load-by-Seurat-v5.R &

可以看到，stromal確實(shí)是非常復(fù)雜，但是張澤民課題組處理的非常好！平滑肌細(xì)胞（Smooth Muscle Cells, SMCs）和周細(xì)胞（Pericytes），這個(gè)混合的單細(xì)胞亞群可以被稱為“混合間充質(zhì)細(xì)胞群”（Mixed Mesenchymal Cell Cluster）或“血管相關(guān)間充質(zhì)細(xì)胞群”（Vascular-Associated Mesenchymal Cell Cluster），具體命名取決于細(xì)胞的來(lái)源和所處的組織微環(huán)境。