单细胞免疫组库VDJ| 从零开始scRepertoire分析,解决真实场景中可能的问题(下)

本文涉及的产品
语种识别,语种识别 100万字符
图片翻译,图片翻译 100张
文档翻译,文档翻译 1千页
简介: 单细胞免疫组库VDJ| 从零开始scRepertoire分析,解决真实场景中可能的问题(下)

本文首发于“生信补给站”公众号   https://mp.weixin.qq.com/s/dqxGfDRS2jg8jV_y1BBonA



scRepertoire-VDJ + RNA


以上即为单细胞免疫组库的常规分析,我们更多的还是与scRNA-seq数据结合,就可以在umap上查看clone的分布情况 以及 基于cluster/celltype展示clonotype的分布情况。

3.1 Seurat标准流程

folders的样本文件夹内为cellranger的3个标准文件,的简单的走一遍标准流程单细胞工具箱|Seurat官网标准流程,此处仅为示例就不考虑批次问题了。

folders=list.files('./')
folders
#[1] "t2_Normal"  "t2_PBMC"    "t3_Center"  "t3_Normal"  "t3_PBMC"    "t4_Center"  "t4_Normal"  "t4_PBMC"    "UT1_Normal"
#批量读取10X单细胞转录组数据文件夹
sceList = lapply(folders,function(folder){ 
  CreateSeuratObject(counts = Read10X(folder), 
                     project = folder )
})
sce.all <- merge(sceList[[1]], 
                 y = c(sceList[[2]],sceList[[3]],sceList[[4]], 
                       sceList[[5]],sceList[[6]],sceList[[7]],sceList[[8]],sceList[[9]]),
                 add.cell.ids = folders, #添加样本名
                 project = "scRNA")
sce.all #查看
#线粒体比例
sce.all[["percent.mt"]] <- PercentageFeatureSet(sce.all, pattern = "^MT-")
#过滤
sce <- subset(sce.all, subset = nFeature_RNA > 500 & 
                nFeature_RNA < 5000 & 
                percent.mt < 20)
VlnPlot(sce, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
#标准化
sce <- NormalizeData(sce)
#高变基因
sce <- FindVariableFeatures(sce, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(sce)
#归一化
sce <- ScaleData(sce, features = all.genes)
#降维聚类
sce <- RunPCA(sce, npcs = 50)
sce <- FindNeighbors(sce, dims = 1:30)
sce <- FindClusters(sce, resolution = 0.5)
sce <- RunUMAP(sce, dims = 1:30)
sce <- RunTSNE(sce, dims = 1:30)
p1 <- DimPlot(sce, reduction = "umap", pt.size=0.5, label = F,repel = TRUE)
p2 <- DimPlot(sce, reduction = "umap",group.by = "orig.ident", pt.size=0.5, label = F,repel = TRUE)
p1 + p2

(1) 添加meta,处理barcode

读取文章提供的meta信息,如下图所示,红框中的不一致的部分需要tidyverse处理一下Tidyverse|数据列的分分合合,一分多,多合一


meta <- read.table("ccRCC_6pat_cell_annotations.txt",
                   sep = "\t",header = T)
meta2 <- meta %>% 
  separate(cell,into = c( "CB","sample","pos"),sep = "_") %>%
  mutate(cell = paste(sample,pos,CB,sep = "_"))
sce@meta.data <- sce@meta.data %>% 
  rownames_to_column("cell") %>%
  inner_join(meta2,by = "cell") %>% 
  column_to_rownames("cell")
head(sce@meta.data)
#                             orig.ident nCount_RNA nFeature_RNA percent.mt RNA_snn_res.0.5 seurat_clusters                 CB
#t2_Normal_AAACCTGCAACACCTA-1  t2_Normal       7863         2244   3.700878               9               9 AAACCTGCAACACCTA-1
#t2_Normal_AAACCTGCACCTCGTT-1  t2_Normal       2194         1188  10.027347               7               7 AAACCTGCACCTCGTT-1
#                             sample    pos   type region Sample   Sample2 cluster               cluster_name     UMAP1      UMAP2
#t2_Normal_AAACCTGCAACACCTA-1     t2 Normal Normal Normal     t2 t2_Normal      13                       cDC2 -9.836647 0.07936505
#t2_Normal_AAACCTGCACCTCGTT-1     t2 Normal Normal Normal     t2 t2_Normal      11 CD45- Vascular Endothelium -7.601253 9.00771236
#                                    Sample_name
#t2_Normal_AAACCTGCAACACCTA-1 Ipi/Nivo resistant
#t2_Normal_AAACCTGCACCTCGTT-1 Ipi/Nivo resistant

(2) sce的细胞数 和 meta信息 不一致

因为前面Seurat流程使用的自己的参数,而meta信息是文章给的。

解决:使用subset函数 提取 meta.data中的cell

sce@meta.data <- sce@meta.data %>% 
  rownames_to_column("Cell") %>% 
  mutate(Cell2 = Cell) %>% 
  column_to_rownames("Cell2")
sce <- subset(sce,Cell %in% sce@meta.data$Cell)


3.2 VDJ结果结合RNA

使用combineExpression函数将VDJ结果和 RNA结果结合。

注:combined的barcode 和  sce的rownames 一定要一致,有出入的使用tidyverse包进行数据处理。

seurat <- combineExpression(combined, sce, 
                            cloneCall="gene", 
                           # group.by = "Sample", 
                            proportion = FALSE, 
                            cloneTypes=c(Single=1, Small=5, Medium=20, Large=100, Hyperexpanded=500))
head(seurat)

因为很多barcode 并没有TCR clone ,因此cloneTypes 含有一些NA是正常的 。

3.3 clonoType umap分布

可以使用DimPlot函数 绘制 clonetype在UUMAP中的分布


colorblind_vector <- colorRampPalette(rev(c("#0D0887FF", "#47039FFF", 
                                            "#7301A8FF", "#9C179EFF", "#BD3786FF", "#D8576BFF",
                                            "#ED7953FF","#FA9E3BFF", "#FDC926FF", "#F0F921FF")))
p111 <- DimPlot(seurat, group.by = "cluster_name",label = T)
p222 <- DimPlot(seurat, group.by = "cloneType",label = T) + #NoLegend() +
  scale_color_manual(values=colorblind_vector(5)) + 
  theme(plot.title = element_blank())
p111 + p222

可以大概看出来CD8细胞中clone更多,符合文献的研究结果。

3.4 highlightClonotypes 高亮特定克隆型

使用highlightClonotypes函数 展示特定序列的克隆型分布

seurat <- highlightClonotypes(seurat, 
                              cloneCall= "aa", 
                              sequence = c("CAVGGRSNSGYALNF;CAASRNNNDMRF_CASSSLGNEQFF", 
                                           "CAASRNNNDMRF_CASSSLGNEQFF"))
DimPlot(seurat, group.by = "highlight") + 
  theme(plot.title = element_blank())

3.5 occupiedscRepertoire cloneType分布

查看CD8T和 CD4T细胞中cloneType的分布情况

seurat_T <- subset(seurat, cluster_name %in% c("CD8A+ Exhausted","CD8A+ Exhausted IEG", 
               "CD8A+ NK-like","CD8A+ Proliferating","CD8A+ Tissue-resident",
               "CD4+ Activated IEG","CD4+ Effector" ,"CD4+ Naive","CD4+ Proliferating" ,
               "CD4+ Treg") )
occupiedscRepertoire(seurat_T, x.axis = "cluster_name")

3.6 alluvialClonotypes

使用alluvialClonotypes函数绘制T细胞种各个celltype之间的clonetype信息的桑葚图


alluvialClonotypes(seurat_T, cloneCall = "gene", 
                   y.axes = c("cluster_name", "cloneType"), 
                   color = "cluster_name")

3.7 getCirclize

使用getCirclize函数绘制circos图

library(circlize)
library(scales)
circles <- getCirclize(seurat, 
                       cloneCall = "gene+nt",
                       groupBy = "cluster_name"
                       )

scRepertoire基于celltype


第二部分中是根据sample计算的clone的一系列指标,如克隆型分布,长度,空间稳态等 。那经过第三部分中和单细胞转录组结合后就可以根据celltype维度进行clone的上述所有统计,这里仅示例几个。

4.1 按照celltype拆分

首先使用expression2List函数将seurat按照cluster_name进行拆分,然后就可进行二中的所有分析 。

注意这里不是官网的split.by ,而要使用group



combined2 <- expression2List(seurat_T, 
                             group = "cluster_name")
length(combined2) 
p01 <- clonalDiversity(combined2, 
                cloneCall = "nt")
p02 <- clonalHomeostasis(combined2, 
                  cloneCall = "nt")
p03 <- clonalOverlap(combined2, 
              cloneCall="aa", 
              method="overlap")
p01
p02
#保存数据(后台获取),以备后用
save(combined,seurat_T,file = "combined_seurat_T.RData")

目录
打赏
0
0
0
0
21
分享
相关文章
R绘图实战|GSEA富集分析图
GSEA(Gene Set EnrichmentAnalysis),即基因集富集分析,它的基本思想是使用预定义的基因,将基因按照在两类样本中的差异表达程度排序,然后检验预先设定的基因集合是否在这个排序表的顶端或者底端富集。
2867 0
R绘图实战|GSEA富集分析图
Visium HD空间数据分析、可视化以及整合 (1)
【8月更文挑战第1天】Visium HD空间数据分析、可视化以及整合 (1)
Visium HD空间数据分析、可视化以及整合 (1)
生信分析代码之前还好好的,怎么就报错了 Error in Ops. data. frame(guide_loc, panel_loc) :'==' only defined for equally-sized data frames
执行 `DimPlot` 函数时遇到错误 `;Error in Ops. data. frame(g guides_loc, panel_loc) : &#39;==&#39; only defined for equally-sized data frames`。解决方案和办法
2008 0
生信分析代码之前还好好的,怎么就报错了 Error in Ops. data. frame(guide_loc, panel_loc) :'==' only defined for equally-sized data frames
scRNA分析|使用CellChat完成细胞通讯分析-简单且可视化出众,代码自取
scRNA分析|使用CellChat完成细胞通讯分析-简单且可视化出众,代码自取
2039 0
单细胞免疫组库VDJ| 从零开始scRepertoire分析,解决真实场景中可能的问题(上)
单细胞免疫组库VDJ| 从零开始scRepertoire分析,解决真实场景中可能的问题
1096 0
PCA 解释方差骤降图(ElbowPlot)拐点的量化识别
本文分享了一种在R语言里面实现量化识别 PCA方差贡献肘部图(ElbowPlot)中拐点 的简单方法,以供学习参考 定
1582 0
2-华大时空组学分析软件 Spateo 空转数据基础分析用法示例
本分分享了使用华大时空组学分析软件 Spateo进行聚类、DE等简单空间转录组分析的用法示例,以供参考
1419 1
转录组下游分析丨利用limma包进行差异表达分析,结果可视化绘制火山图和热图
转录组下游分析丨利用limma包进行差异表达分析,结果可视化绘制火山图和热图
AI助理

你好,我是AI助理

可以解答问题、推荐解决方案等

登录插画

登录以查看您的控制台资源

管理云资源
状态一览
快捷访问