參考教程鏈接：
- https://github.com/dmcable/spacexr
- 官方教程： https://github.com/dmcable/spacexr/tree/master/vignettes
- 文獻(xiàn)解讀：10X單細(xì)胞空間聯(lián)合分析之十（RCTD） http://www.reibang.com/p/92a6df0dcb1b
RCTD可以將單細(xì)胞類型或細(xì)胞類型混合分配到空間轉(zhuǎn)錄組spots上。RCTD 有三種模式：
- doublet mode：它為每個(gè)Spot分配1-2種細(xì)胞類型，推薦用于高空間分辨率的技術(shù)陪拘，如Slide-seq和MERFISH
- full mode：它為每個(gè)Spot分配任意數(shù)量的細(xì)胞類型，推薦用于空間分辨率較低的技術(shù)苍糠，如100微米分辨率的Visium
- multi mode：doublet mode的擴(kuò)展，可以在每個(gè)點(diǎn)上發(fā)現(xiàn)兩種以上的細(xì)胞類型啤誊，作為full mode的替代選項(xiàng)
每種模態(tài)都有一個(gè)對(duì)應(yīng)的教程
今天先來學(xué)習(xí)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
穩(wěn)健細(xì)胞類型分解(Robust Cell Type Decomposition岳瞭，簡(jiǎn)稱RCTD)是一種從空間轉(zhuǎn)錄組數(shù)據(jù)中學(xué)習(xí)細(xì)胞類型的統(tǒng)計(jì)方法。
在本次示例中蚊锹，我們將為小腦Slide-seq數(shù)據(jù)集反卷積注釋細(xì)胞類型瞳筏。教程使用帶注釋的snRNA-seq小腦數(shù)據(jù)集定義細(xì)胞類型。
Note：參考也可以是單細(xì)胞數(shù)據(jù)集或細(xì)胞類型特定的bulkRNA-seq數(shù)據(jù)集牡昆。
數(shù)據(jù)預(yù)處理
RCTD需要兩個(gè)數(shù)據(jù)：帶有注釋的單細(xì)胞轉(zhuǎn)錄組參考數(shù)據(jù)姚炕，需要反卷積注釋的空間轉(zhuǎn)錄組數(shù)據(jù)。
reference數(shù)據(jù)
首先是單細(xì)胞參考數(shù)據(jù)集丢烘。reference使用RCTD包中Reference函數(shù)柱宦，這個(gè)函數(shù)需要三個(gè)參數(shù)：
- counts：為矩陣或dgCmatrix對(duì)象。行名是基因播瞳，列名代表barcode/細(xì)胞名掸刊。counts應(yīng)該是未轉(zhuǎn)換的原始counts數(shù)據(jù)
- cell_types：為帶有細(xì)胞命名的(通過細(xì)胞barcode)細(xì)胞類型因子。因子的“水平”是可能的細(xì)胞類型標(biāo)識(shí)赢乓。
- nUMI：可選忧侧，具有barcode的每個(gè)像素中總counts or UMI列表。如果沒有提供骏全，nUMI將被假定為出現(xiàn)在每個(gè)像素上的總數(shù)。
reference可能來自各種數(shù)據(jù)類型尼斧，但需要將其加載到R環(huán)境中姜贡。
在本教程中，我們的reference作為兩個(gè)csv文件存儲(chǔ)在' reference /Vignette '文件夾中:
- meta_data.csv：CSV文件(有3列棺棵，標(biāo)題為“barcode”楼咳、“cluster”和“nUMI”)熄捍，包含每個(gè)細(xì)胞的nUMI和單元格類型分配。
- dge.csv：標(biāo)準(zhǔn)10x格式的基因表達(dá)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已經(jīng)構(gòu)建好。
空間轉(zhuǎn)錄組數(shù)據(jù)
接下來苹熏，加載空間數(shù)據(jù)為SpatialRNA對(duì)象碟贾。SpatialRNA函數(shù)需要三個(gè)參數(shù)：
- coords：數(shù)據(jù)框或者矩陣對(duì)象，為空間像素坐標(biāo)轨域，行名為每個(gè)像素的barcode袱耽，列為x和y
- counts：矩陣或者dgCmatrix對(duì)象，行名是基因干发，列名是每個(gè)像素的barcode朱巨，counts值沒有進(jìn)行轉(zhuǎn)換。
- nUMI：可選
構(gòu)建SpatialRNA有不同的方法枉长，本教程為兩個(gè)文件：
- BeadLocationsForR.csv：CSV文件冀续，三列，列名為"barcodes”, “xcoord”, and “ycoord”必峰，包含每個(gè)像素的空間坐標(biāo)洪唐；
- MappedDGEForR.csv：每個(gè)像素的表達(dá)矩陣文件
##=============================== 讀取空間數(shù)據(jù)
# 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))
結(jié)果圖：可以看到空間數(shù)據(jù)的每個(gè)像素中UMI的數(shù)值分布。
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')
結(jié)果圖：使用x自点，y坐標(biāo)展示每個(gè)像素中的UMI分布
創(chuàng)建RCTD對(duì)象
使用create.RCTD函數(shù)創(chuàng)建桐罕。以下為幾個(gè)有用的參數(shù)：
- max_cores：對(duì)于并行處理，使用的線程數(shù)
- gene_cutoff, fc_cutoff, gene_cutoff_reg, fc_cutoff_reg：差異表達(dá)基因篩選閾值
- UMI_min, UMI_max：每個(gè)像素點(diǎn)的閾值篩選
運(yùn)行RCTD
可以使用run.RCTD函數(shù)桂敛，這個(gè)函數(shù)等價(jià)于：fitBulk,choose_sigma_c, andfitPixels這三個(gè)函數(shù)功炮。doublet_mode可以選擇上面提到的三種反卷積模式。
myRCTD <- create.RCTD(puck, reference, max_cores = 10)
myRCTD <- run.RCTD(myRCTD, doublet_mode = 'doublet')
str(myRCTD)
myRCTD結(jié)構(gòu)：
RCTD結(jié)果
RCTD結(jié)果保存在myRCTD對(duì)象中的@results中
- @results$weights：是一個(gè)數(shù)據(jù)框术唬，為full mode模式下每一個(gè)spots（即像素）中細(xì)胞類型的權(quán)重
- @results$results_df：doublet_mode的結(jié)果
- spot_class：因子類型薪伏，RCTD的分類，singlet表示每個(gè)spot有一種細(xì)胞類型粗仓，doublet_certain表示每個(gè)spot有兩種細(xì)胞類型嫁怀，doublet_uncertain表示spot上有2種細(xì)胞類型，但只確定有1種借浊，reject表示沒有預(yù)測(cè)到結(jié)果塘淑。
- first_type：因子類型，預(yù)測(cè)到的第一種細(xì)胞類型
- second_type：因子類型蚂斤，預(yù)測(cè)到的第一種細(xì)胞類型
- @results$weights_doublet：doublet_mode下每種細(xì)胞類型的權(quán)重
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結(jié)果示例如下：
接下來對(duì)結(jié)果可視化繪圖：
# make the plots
# 繪制full_mode模式下每種細(xì)胞類型的可信權(quán)重 (saved as
# 'results/cell_type_weights_unthreshold.pdf')
plot_weights(cell_type_names, spatialRNA, resultsdir, norm_weights)
# 繪制full_mode模式下每種細(xì)胞類型的權(quán)重
# 這里每種細(xì)胞類型一幅圖存捺，點(diǎn)表示空間上的一個(gè)像素或者spot，顏色為權(quán)重 (saved as
# 'results/cell_type_weights.pdf')
plot_weights_unthreshold(cell_type_names, spatialRNA, resultsdir, norm_weights)
# 繪制full_mode模式下每種細(xì)胞類型預(yù)測(cè)到的spots數(shù) (saved as
# 'results/cell_type_occur.pdf')
plot_cond_occur(cell_type_names, resultsdir, norm_weights, spatialRNA)
# 繪制doublet_mode模式下每種細(xì)胞類型的權(quán)重 (saved as
# 'results/cell_type_weights_doublets.pdf')
plot_weights_doublet(cell_type_names, spatialRNA, resultsdir, results$weights_doublet, results$results_df)
其中，繪制full_mode模式下每種細(xì)胞類型預(yù)測(cè)到的spots數(shù)結(jié)果圖如下：
從這里可以看得出捌治，空間數(shù)據(jù)預(yù)測(cè)出來的主要為細(xì)胞類型10岗钩，其次是18。
# 所有細(xì)胞類型的map結(jié)果 (saved as
# 'results/all_cell_types.pdf')
plot_all_cell_types(results$results_df, spatialRNA@coords, cell_type_names, resultsdir)
所有細(xì)胞類型的map結(jié)果圖：
繪制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)
# 對(duì)每種細(xì)胞類型繪制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出現(xiàn)的堆積柱狀圖 (saved as
# 'results/doublet_stacked_bar.pdf')
plot_doub_occur_stack(doub_occur, resultsdir, cell_type_names)
doublet_stacked_bar結(jié)果圖：
感覺這個(gè)軟件的結(jié)果圖不是很好看肖油，可調(diào)整度也不高兼吓，我覺得跟教程的數(shù)據(jù)有關(guān)系，后面重新找個(gè)數(shù)據(jù)看看森枪。
教程的空間數(shù)據(jù)是slide-seq的數(shù)據(jù)：
找了個(gè)10x Visium的文獻(xiàn)來看看：
ref1：Cell2location maps fine-grained cell types in spatial transcriptomics ( https://doi.org/10.1038/s41587-021-01139-4 )
