空间转录组2022||空间数据反卷积RCTD分析：细胞类型鉴定

2022-11-19 本文已影响0人信你个鬼

参考教程链接：

https://github.com/dmcable/spacexr
官方教程： https://github.com/dmcable/spacexr/tree/master/vignettes
文献解读：10X单细胞空间联合分析之十（RCTD） https://www.jianshu.com/p/92a6df0dcb1b

RCTD可以将单细胞类型或细胞类型混合分配到空间转录组spots上。RCTD 有三种模式：

doublet mode：它为每个Spot分配1-2种细胞类型，推荐用于高空间分辨率的技术，如Slide-seq和MERFISH
full mode：它为每个Spot分配任意数量的细胞类型，推荐用于空间分辨率较低的技术，如100微米分辨率的Visium
multi mode：doublet mode的扩展，可以在每个点上发现两种以上的细胞类型，作为full mode的替代选项

每种模态都有一个对应的教程

今天先来学习Doublet mode: spatial transcriptomics vignette. Also, most other vignettes use doublet mode.

包安装

# install.packages("devtools")
devtools::install_github("dmcable/spacexr", build_vignettes = FALSE)
library(spacexr)
library(Matrix)

Introduction

稳健细胞类型分解(Robust Cell Type Decomposition，简称RCTD)是一种从空间转录组数据中学习细胞类型的统计方法。
在本次示例中，我们将为小脑Slide-seq数据集反卷积注释细胞类型。教程使用带注释的snRNA-seq小脑数据集定义细胞类型。
Note：参考也可以是单细胞数据集或细胞类型特定的bulkRNA-seq数据集。

数据预处理

RCTD需要两个数据：带有注释的单细胞转录组参考数据，需要反卷积注释的空间转录组数据。

reference数据

首先是单细胞参考数据集。reference使用RCTD包中Reference函数，这个函数需要三个参数：

counts：为矩阵或dgCmatrix对象。行名是基因，列名代表barcode/细胞名。counts应该是未转换的原始counts数据
cell_types：为带有细胞命名的(通过细胞barcode)细胞类型因子。因子的“水平”是可能的细胞类型标识。
nUMI：可选，具有barcode的每个像素中总counts or UMI列表。如果没有提供，nUMI将被假定为出现在每个像素上的总数。

reference可能来自各种数据类型，但需要将其加载到R环境中。

在本教程中，我们的reference作为两个csv文件存储在' reference /Vignette '文件夹中:

meta_data.csv：CSV文件(有3列，标题为“barcode”、“cluster”和“nUMI”)，包含每个细胞的nUMI和单元格类型分配。
dge.csv：标准10x格式的基因表达CSV文件。

### Load in/preprocess your data, this might vary based on your file type
# directory for the reference
refdir <- system.file("extdata",'Reference/Vignette',package = 'spacexr') 
# load in counts matrix
counts <- read.csv(file.path(refdir,"dge.csv")) 
# Move first column to rownames
rownames(counts) <- counts[,1]
counts[,1] <- NULL
counts[1:4,1:4]


# load in meta_data (barcodes, clusters, and nUMI)
meta_data <- read.csv(file.path(refdir,"meta_data.csv")) 
head(meta_data)

# create cell_types named list
cell_types <- meta_data$cluster
names(cell_types) <- meta_data$barcode 
# convert to factor data type
cell_types <- as.factor(cell_types) 
head(cell_types)

# create nUMI named list
nUMI <- meta_data$nUMI
names(nUMI) <- meta_data$barcode 
head(nUMI)

### Create the Reference object
reference <- Reference(counts, cell_types, nUMI)
str(reference)

## Examine reference object (optional)
#observe Digital Gene Expression matrix
print(dim(reference@counts)) 
#> [1] 384 475

#number of occurences for each cell type
table(reference@cell_types) 

## Save RDS object (optional)
saveRDS(reference, file.path(refdir,'SCRef.rds'))

至此，reference已经构建好。

空间转录组数据

接下来，加载空间数据为SpatialRNA对象。SpatialRNA函数需要三个参数：

coords：数据框或者矩阵对象，为空间像素坐标，行名为每个像素的barcode，列为x和y
counts：矩阵或者dgCmatrix对象，行名是基因，列名是每个像素的barcode，counts值没有进行转换。
nUMI：可选

构建SpatialRNA有不同的方法，本教程为两个文件：

BeadLocationsForR.csv：CSV文件，三列，列名为"barcodes”, “xcoord”, and “ycoord”，包含每个像素的空间坐标；
MappedDGEForR.csv：每个像素的表达矩阵文件

##=============================== 读取空间数据
# directory for sample Slide-seq dataset
datadir <- system.file("extdata",'SpatialRNA/Vignette',package = 'spacexr') 
counts <- read.csv(file.path(datadir,"MappedDGEForR.csv"))
rownames(counts) <- counts[,1]; counts[,1] <- NULL
counts[1:4,1:4]

coords <- read.csv(file.path(datadir,"BeadLocationsForR.csv"))
rownames(coords) <- coords$barcodes
coords$barcodes <- NULL
head(coords)

# In this case, total counts per pixel is nUMI
nUMI <- colSums(counts) 
head(nUMI)

### Create SpatialRNA object
puck <- SpatialRNA(coords, counts, nUMI)
str(puck)


### Create SpatialRNA object
puck <- SpatialRNA(coords, counts, nUMI)
str(puck)

## Examine SpatialRNA object (optional)
print(dim(puck@counts))

# histogram of log_2 nUMI
hist(log(puck@nUMI,2))

结果图：可以看到空间数据的每个像素中UMI的数值分布。

image-20221111222523662.png

print(head(puck@coords)) # start of coordinate data.frame
barcodes <- colnames(puck@counts) # pixels to be used (a list of barcode names). 

# This list can be restricted if you want to crop the puck e.g. 
# puck <- restrict_puck(puck, barcodes) provides a basic plot of the nUMI of each pixel
# on the plot:
p <- plot_puck_continuous(puck, barcodes, puck@nUMI, ylimit = c(0,round(quantile(puck@nUMI,0.9))), title ='plot of nUMI')

结果图：使用x，y坐标展示每个像素中的UMI分布

image-20221111222947585.png

创建RCTD对象

使用create.RCTD函数创建。以下为几个有用的参数：

max_cores：对于并行处理，使用的线程数
gene_cutoff, fc_cutoff, gene_cutoff_reg, fc_cutoff_reg：差异表达基因筛选阈值
UMI_min, UMI_max：每个像素点的阈值筛选

运行RCTD

可以使用run.RCTD函数，这个函数等价于：fitBulk,choose_sigma_c, andfitPixels这三个函数。doublet_mode可以选择上面提到的三种反卷积模式。

myRCTD <- create.RCTD(puck, reference, max_cores = 10)
myRCTD <- run.RCTD(myRCTD, doublet_mode = 'doublet')
str(myRCTD)

myRCTD结构：

image-20221112182253019.png

RCTD结果

RCTD结果保存在myRCTD对象中的@results中

@results$weights：是一个数据框，为full mode模式下每一个spots（即像素）中细胞类型的权重
@results$results_df：doublet_mode的结果
- spot_class：因子类型，RCTD的分类，singlet表示每个spot有一种细胞类型，doublet_certain表示每个spot有两种细胞类型，doublet_uncertain表示spot上有2种细胞类型，但只确定有1种，reject表示没有预测到结果。
- first_type：因子类型，预测到的第一种细胞类型
- second_type：因子类型，预测到的第一种细胞类型
@results$weights_doublet：doublet_mode下每种细胞类型的权重

results <- myRCTD@results

# normalize the cell type proportions to sum to 1.
norm_weights = normalize_weights(results$weights) 

#list of cell type names
cell_type_names <- myRCTD@cell_type_info$info[[2]] 
spatialRNA <- myRCTD@spatialRNA

## you may change this to a more accessible directory on your computer.
resultsdir <- 'RCTD_Plots' 
dir.create(resultsdir)

results$results_df结果示例如下：

image-20221112184204705.png

接下来对结果可视化绘图：

# make the plots 
# 绘制full_mode模式下每种细胞类型的可信权重 (saved as 
# 'results/cell_type_weights_unthreshold.pdf')
plot_weights(cell_type_names, spatialRNA, resultsdir, norm_weights) 

# 绘制full_mode模式下每种细胞类型的权重
# 这里每种细胞类型一幅图，点表示空间上的一个像素或者spot，颜色为权重 (saved as
# 'results/cell_type_weights.pdf')
plot_weights_unthreshold(cell_type_names, spatialRNA, resultsdir, norm_weights) 

# 绘制full_mode模式下每种细胞类型预测到的spots数 (saved as 
# 'results/cell_type_occur.pdf')
plot_cond_occur(cell_type_names, resultsdir, norm_weights, spatialRNA)

# 绘制doublet_mode模式下每种细胞类型的权重 (saved as 
# 'results/cell_type_weights_doublets.pdf')
plot_weights_doublet(cell_type_names, spatialRNA, resultsdir, results$weights_doublet, results$results_df)

其中，绘制full_mode模式下每种细胞类型预测到的spots数结果图如下：

从这里可以看得出，空间数据预测出来的主要为细胞类型10，其次是18。

image-20221112192742832.png

# 所有细胞类型的map结果 (saved as 
# 'results/all_cell_types.pdf')
plot_all_cell_types(results$results_df, spatialRNA@coords, cell_type_names, resultsdir)

所有细胞类型的map结果图：

image-20221112194400757.png

绘制doublets：

# doublets
# obtain a dataframe of only doublets
doublets <- results$results_df[results$results_df$spot_class == "doublet_certain",] 

# 绘制所有的doublets (saved as 
# 'results/all_doublets.pdf')
plot_doublets(spatialRNA, doublets, resultsdir, cell_type_names) 

# 对每种细胞类型绘制doublets (saved as 
# 'results/all_doublets_type.pdf')
plot_doublets_type(spatialRNA, doublets, resultsdir, cell_type_names) 

# a table of frequency of doublet pairs 
doub_occur <- table(doublets$second_type, doublets$first_type) 

# 绘制doublet出现的堆积柱状图 (saved as   
# 'results/doublet_stacked_bar.pdf')
plot_doub_occur_stack(doub_occur, resultsdir, cell_type_names)

doublet_stacked_bar结果图：

image-20221112194946431.png

感觉这个软件的结果图不是很好看，可调整度也不高，我觉得跟教程的数据有关系，后面重新找个数据看看。

教程的空间数据是slide-seq的数据：

img

找了个10x Visium的文献来看看：

ref1：Cell2location maps fine-grained cell types in spatial transcriptomics ( https://doi.org/10.1038/s41587-021-01139-4 )

image-20221112222920831.png