【原】胃癌單細胞數(shù)據(jù)集GSE163558復現(xiàn)(十一)：T細胞亞群細分

健明 2024-08-09 發(fā)布于廣東

展開全文

前言

Hello小伙伴們大家好，我是生信技能樹的小學徒”我才不吃蛋黃“。今天是胃癌單細胞數(shù)據(jù)集GSE163558復現(xiàn)系列第十一期。第十期我們使用通過CopyKAT分析計算CNV評估上皮細胞良惡性。本期，我們將對T細胞亞群進行細分。

1.背景介紹

德國哲學家萊布尼茨說過:“世上沒有兩片完全相同的樹葉。”物種是有其多樣性的。對于細胞來說同樣如此。從理論上講，我們可以對細胞亞群進行無限分群。在初次分群時，可以不斷調(diào)整分辨率，同樣在細胞注釋之后，我們也可以進行亞群再細分。其實分群也是考驗我們對分群“度”的把控以及生物學背景知識的儲備。細胞注釋前我們要確保分群數(shù)量足夠多，而這個數(shù)量通常是要超過該組織理論的細胞類型數(shù)量。在進行亞群細分時，更加需要強大的生物學背景知識儲備以及對分辨率和再分群次數(shù)的“度”的精準把握。細胞亞群的分群的依據(jù)包括：1.細胞異質(zhì)性（每個細胞都有獨特的表達模式和功能，都有自己特有的基因）；2.細胞共性（同一類型的細胞都有近似的表達模式）；3.生物學基礎知識（基于已有的知識，對細胞進行分類鑒定）。

T淋巴細胞（T lymphocyte）簡稱T細胞，是由來源于骨髓的淋巴干細胞，在胸腺中分化、發(fā)育成熟后，通過淋巴和血液循環(huán)而分布到全身的免疫器官和組織中發(fā)揮免疫功能。T細胞分類方法有多種：按細胞表面分化抗原 (CD)的不同,可分為CD4+和CD8+兩大亞群;按T細胞表面受體 (TCR)的不同,可分為αβT細胞γδT細胞;按功能可分為輔助性T細胞 (Th 細胞)、抑制性T細胞 (Ts細胞)、細胞毒T細胞 (CTL或Tc細胞)和遲發(fā)型超敏反應T細胞 (TDTH細胞);按對抗原應答的不同,分為初始T細胞 (naive T cell)、活化的T細胞 (activated T cell)和記憶性T細胞 (memory T cell);以及區(qū)別于傳統(tǒng)T細胞的NKT細胞等等。在第一次分群注釋中，我們只注釋到了“初始的”T細胞。在這里，我們對T細胞亞群進行了粗“細分”。

2.數(shù)據(jù)分析

2.1 導入數(shù)據(jù)

首先提取第三期亞群注釋后的T細胞數(shù)據(jù)：

rm(list=ls())
options(stringsAsFactors = F)
library(Seurat)
library(ggplot2)
library(clustree)
library(cowplot)
library(dplyr)
dir.create("9-T")
setwd("9-T/") 
getwd()
set.seed(12345)
sce.all=readRDS( "../3-Celltype/sce_celltype.rds")
table(sce.all$celltype)
sce1 = sce.all[, sce.all$celltype %in% c( 'T' )]

2.2 降維分群聚類

細分亞群的流程很簡單，和大致分群一樣，就是對亞群的Seurat對象再次走降維分群聚類。

Seurat V5標準流程：

LayerData(sce1, assay = "RNA", layer = "counts")
sce1 <- JoinLayers(sce1)
as.data.frame(sce1@assays$RNA$counts[1:10, 1:2])
head(sce1@meta.data, 10)
table(sce1$orig.ident) 
sce = sce1
sce <- NormalizeData(sce, normalization.method =  "LogNormalize",  
                     scale.factor = 1e4)
GetAssay(sce,assay = "RNA")
sce <- FindVariableFeatures(sce, 
                            selection.method = "vst", nfeatures = 2000)  
sce <- ScaleData(sce) 
sce <- RunPCA(object = sce, pc.genes = VariableFeatures(sce))

DimHeatmap(sce, dims = 1:12, cells = 100, balanced = TRUE)

ElbowPlot(sce)

選取0.1分辨率，對T細胞進行聚類分群（為什么選取0.1，歡迎大家在評論區(qū)或者微信群留言討論）:

sce <- FindNeighbors(sce, dims = 1:15)
sce <- FindClusters(sce, resolution = 0.1)
table(sce@meta.data$RNA_snn_res.0.1)  
set.seed(321)
sce <- RunTSNE(object = sce, dims = 1:15, do.fast = TRUE)

sce <- RunUMAP(object = sce, dims = 1:5, do.fast = TRUE)

mycolors <-c('#E64A35','#4DBBD4' ,'#01A187'  ,'#3C5588'  ,'#F29F80'  ,
             '#8491B6','#91D0C1','#7F5F48','#AF9E85','#4F4FFF','#CE3D33',
             '#739B57','#EFE685','#446983','#BB6239','#5DB1DC','#7F2268','#6BD66B','#800202','#D8D8CD','pink'
)

p = DimPlot(sce,reduction = "tsne",label=T,cols = mycolors) 
p

2.3 SingleR注釋

SingleR是一種用于單細胞 RNA 測序（scRNAseq）數(shù)據(jù)的自動注釋方法。給定具有已知標簽的樣本（單細胞或塊）的參考數(shù)據(jù)集，它根據(jù)與參考的相似性標記測試數(shù)據(jù)集中的新細胞。

singleR使用步驟：

需要一個數(shù)據(jù)庫文件，構(gòu)建SingleR進行單細胞亞群命名的參考數(shù)據(jù)庫

使用SingleR包里面的SingleR函數(shù)即可把數(shù)據(jù)庫里面的細胞亞群注釋信息映射到需要命名的單細胞轉(zhuǎn)錄組數(shù)據(jù)集里面。

成功的運行了SingleR包里面的SingleR函數(shù)之后，就可以拿到每個單細胞的具體的身份信息

詳情介紹請見公眾號前期推文：SingleR及數(shù)據(jù)庫資源包celldex簡介。

SingleR包的安裝方式有以下兩種：

使用devtools包進行安裝：

devtools::install_github('dviraran/SingleR')

或者下載本地安裝：

# 安裝依賴包
install.packages(c("outliers", "pbmcapply", "doFuture"))
BiocManager::install(c("GSVA","singscore"))
# 安裝SingleR包
install.packages("~/SingleR.tar.gz", repos = NULL, type = "source")

celldex包是singleR自動注釋需要的數(shù)據(jù)庫，目前已經(jīng)上傳到Bioconductor，安裝方法如下：

BiocManager::install("celldex")

加載R包：

library(SingleR)
library(celldex)
library(dplyr)
library(stringr)
library(pheatmap)
library(ReactomeGSA)
library(ggplot2)
library(singleseqgset)
library(devtools)

開始SingleR注釋：

getwd()
load('/home/data/t020505/GSE163558-GC代碼版/9-T/hpca.RData')
load('/home/data/t020505/GSE163558-GC代碼版/9-T/bpe.RData')
unique(hpca.se$label.main)
unique(hpca.se$label.fine)
unique(bpe.se$label.main)
unique(bpe.se$label.fine)
#bpe.se <- BlueprintEncodeData()
#save(bpe.se,file = 'bpe.RData')
 str(sce)
anno <- SingleR(sce@assays$RNA$data,
                ref = list(BP=bpe.se,HPCA=hpca.se),
                labels = list(bpe.se$label.fine,hpca.se$label.main),
                clusters = sce@meta.data$seurat_clusters
)

plotScoreHeatmap(anno,clusters = anno@rownames,show_colnames = T)
table(anno$labels)

使用DimPlot可視化亞群細分之后的情況:

celltype = data.frame(ClusterID=rownames(anno), 
                      celltype=anno$labels, 
                      stringsAsFactors = F) 

sce@meta.data$singleR = celltype[match(sce@meta.data$seurat_clusters,celltype$ClusterID),'celltype']
table(sce$singleR)
p1 = DimPlot(sce, reduction = "tsne",
        group.by = "singleR",label = T,cols = mycolors,label.box=T) 
p1

我們可以發(fā)現(xiàn)，Treg由3群和4群組成，明顯還可以細分。此外，由于這里是從T細胞中注釋到了一個B細胞（6群），可以看看標記基因表達情況：

Idents(sce) = sce$singleR
VlnPlot(sce,
        features = c("CD3D","CD3E","CD4","CD8A",'MS4A1','IGHG1','MZB1', 'CD79A'),
        pt.size = 0,
        ncol = 4,
        cols=mycolors)