10X scRNA免疫治療學(xué)習(xí)筆記-4-細(xì)胞亞群的生物學(xué)命名

劉小澤寫于19.10.15
筆記目的：根據(jù)生信技能樹的單細(xì)胞轉(zhuǎn)錄組課程探索10X Genomics技術(shù)相關(guān)的分析
課程鏈接在：http://jm.grazy.cn/index/mulitcourse/detail.html?cid=55
第二單元第8講：細(xì)胞亞群的生物學(xué)命名

上一次我們好不容易得到了這個(gè)1.9G的RData，也就是作者自己做出來的Seurat對象，那么我們要怎么利用它去進(jìn)一步探索呢？

加載上次運(yùn)行的結(jié)果

首先還是三步走

rm(list = ls()) 
options(warn=-1)
suppressMessages(library(Seurat))

然后加載進(jìn)來之前保存的5.1G PBMC對象

start_time <- Sys.time()
load('./patient1.PBMC.output.Rdata')
end_time <- Sys.time()
end_time - start_time
# Time difference of 12.6741 secs

重溫上次的聚類結(jié)果

# 使用Seurat 2.3.4版本
colP<-c('green4', 
        'pink', 
        '#FF7F00', 
        'orchid', 
        '#99c9fb', 
        'dodgerblue2', 
        'grey30', 
        'yellow', 
        'grey60', 
        'grey', 
        'red', 
        '#FB9A99', 
        'black'
)
TSNEPlot(PBMC, 
         colors.use =  colP,
         do.label = T)

開始新的探索

看看作者整合的數(shù)據(jù)對批次的處理

也就是把四個(gè)時(shí)間點(diǎn)映射到上面的tsne坐標(biāo)中，并且理論上應(yīng)該是：每群細(xì)胞都覆蓋到四個(gè)時(shí)間點(diǎn)

TSNEPlot(PBMC,group.by = "TimePoints")

再用table對比一下

> table(PBMC@meta.data$TimePoints,PBMC@ident)
               
                  0   1   2   3   4   5   6   7   8   9  10  11  12
  PBMC_ARD614   665 726 572 559 420 302 457 313 283 123  17  11  68
  PBMC_EarlyD27  43 173 245  85 120 110 543  59  91  29   7   3  84
  PBMC_Pre      369 527 197  93 146 393   4  76  17  48  25 187   0
  PBMC_RespD376 800 433 555 677 636 516 119 324 204 200 170  11  39

可視化一些marker基因

這些marker基因也是來源于文章的Supp Fig.7，基于他們對免疫知識的了解

allGenes = row.names(PBMC@raw.data)
markerGenes <- c(
  "CD3D",
  "CD3E",
  "TRAC",
  "IL7R",
  "GZMA",
  "FCGR3A",
  "CD14",
  "MS4A1",
  "FCER1A" 
)
# 判斷這些marker是不是存在表達(dá)矩陣中
> markerGenes %in% allGenes
[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
# 選擇性保存pdf格式
# pdf('patient1_pBMC_marker_FeaturePlot.pdf', width=10, height=15)
FeaturePlot(object = PBMC, 
            features.plot =markerGenes, 
            cols.use = c("grey", "blue"), 
            reduction.use = "tsne")

利用上面marker基因在不同細(xì)胞群的特殊表達(dá)，其實(shí)就能得到每個(gè)群的細(xì)胞命名，不過這些marker基因的獲得以及和細(xì)胞名稱的對應(yīng)，是需要一段時(shí)間的研究才能得到的。我們這里只是展示如何操作

賦予每個(gè)cluster細(xì)胞類型

根據(jù)我們得到的cluster和原文的命名

就可以對應(yīng)得到一個(gè)細(xì)胞名稱列表：

cat >celltype-patient1-PBMC.txt
0 "B cells"
1 "CD4+ T cells"
2 "Naive memory T cells"
3 "Classical monicytes"
4 "CD8+ effector T cells"
5 "NK cells"
6 "rm1"
7 "Non-classical monocytes"
8 "Dendritic cells"
9 "rm2"
10  "CD8+ cytotoxic T cells"
11  "Myeloid cells"
12  "rm3"

假如我們是先有了這個(gè)列表，核心就是要將分群的clsuter數(shù)字與細(xì)胞名稱和顏色對應(yīng)起來

利用Seurat V3

版本3的優(yōu)勢是可以直接提取出seurat對象的cluster數(shù)字，并且提供了RenameIdents函數(shù)，直接進(jìn)行轉(zhuǎn)換

a=read.table('celltype-patient1-PBMC.txt')
new.cluster.ids <- as.character(a[,2])
names(new.cluster.ids) <- levels(PBMC_V3)
PBMC_V3 <- RenameIdents(PBMC_V3, new.cluster.ids)

DimPlot(PBMC_V3, reduction = "tsne", label = TRUE, pt.size = 0.5, cols = colP) + NoLegend()
利用Seurat V2

版本2就需要自己去對應(yīng)：首先要了解它的分群信息存儲在PBMC@ident中，然后要將全部12874個(gè)細(xì)胞的分群編號與celltype-patient1-PBMC.txt表中的第一列編號對應(yīng)（只有這樣，才能和表中的第二列對應(yīng)上）

用到對應(yīng)關(guān)系時(shí)，首先思考match能不能做到；如果要做，需要準(zhǔn)備什么；對應(yīng)關(guān)系搞清楚

# 先與表中第一列對應(yīng)
match(as.numeric(as.character(PBMC@ident)),a[,1])

# 第一點(diǎn)需要注意的是：為什么先用as.character后用as.numeric，而不是直接用as.numeric？
# 原因就是PBMC@ident存儲的是因子型變量，直接取只會得到它們的位置信息，而不是真實(shí)的分群信息
> head(as.numeric(PBMC@ident))
[1] 2 1 2 2 2 6
> head(as.character(PBMC@ident))
[1] "1" "0" "1" "1" "1" "5"
> head(as.numeric(as.character(PBMC@ident)))
[1] 1 0 1 1 1 5

# 第二點(diǎn)需要注意的是：match函數(shù)的規(guī)則是，A要在B中找到對應(yīng)位置，那么就是 match(A,B)

接下來就可以得到對應(yīng)的第二列，也就是細(xì)胞名稱

labels=a[match(as.numeric(as.character(PBMC@ident)),a[,1]),2]

# 檢查一下，主要看數(shù)量的對應(yīng)
> table(labels)
labels
                B cells            CD4+ T cells  CD8+ cytotoxic T cells 
                   1877                    1859                     219 
  CD8+ effector T cells     Classical monicytes         Dendritic cells 
                   1322                    1414                     595 
          Myeloid cells                NK cells    Naive memory T cells 
                    212                    1321                    1569 
Non-classical monocytes                     rm1                     rm2 
                    772                    1123                     400 
                    rm3 
                    191 
> table(PBMC@ident)

   0    1    2    3    4    5    6    7    8    9   10   11   12 
1877 1859 1569 1414 1322 1321 1123  772  595  400  219  212  191

添加到metadata中，方便后面使用

PBMC@meta.data$labels=labels

TSNEPlot(PBMC, group.by = 'labels',
         colors.use =  colP,
         do.label = T)

但很明顯，顏色標(biāo)記和原文不同。因?yàn)橹爸皇菍?yīng)了分群的編號和細(xì)胞名稱，此外還需要修改顏色的順序

# 修改顏色順序，思想就是：將現(xiàn)在的細(xì)胞名與表中的細(xì)胞名對應(yīng)一下，然后這個(gè)順序就是顏色出現(xiàn)的順序
colP=colP[match(levels(as.factor(labels)),a[,2])]
TSNEPlot(PBMC, group.by = 'labels',
         colors.use =  colP,
         do.label = T)

再按時(shí)間拆分分群結(jié)果

做出文章的這張圖

首先得到我們的四個(gè)時(shí)間點(diǎn)
> TimePoints = PBMC@meta.data$TimePoints
> table(TimePoints)
TimePoints
  PBMC_ARD614 PBMC_EarlyD27      PBMC_Pre PBMC_RespD376 
         4516          1592          2082          4684
然后用一個(gè)函數(shù)SubsetData

舉一個(gè)例子，繪制Pre時(shí)期的分群圖

# SubsetData函數(shù)其實(shí)也是利用TRUE/FALSE取子集
PBMC_Pre  = SubsetData(PBMC,TimePoints =='PBMC_Pre')
TSNEPlot(PBMC_Pre, 
         colors.use = c('green4', 'pink', '#FF7F00', 'orchid', '#99c9fb', 'dodgerblue2', 'grey30', 'yellow', 'grey60', 'grey', 'red', '#FB9A99', 'black'),
         do.label = F)
# ggsave('PBMC_Pre_tSNE.pdf')
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請結(jié)合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡書系信息發(fā)布平臺，僅提供信息存儲服務(wù)。

友情鏈接更多精彩內(nèi)容