SPOTlight

這個(gè)是官網(wǎng)示例數(shù)蒙兰,把自己的數(shù)據(jù)填進(jìn)去就行异吻。
這個(gè)好的一點(diǎn)是可以看出來每一個(gè)spot中包括的細(xì)胞類型伤提,我簡直太需要這個(gè)數(shù)據(jù)了题涨。


title: "Spatial Transcriptomics Deconvolution with SPOTlight"
date: "r BiocStyle::doc_date()"
author:

  • name: Marc Elosua-Bayes
    affiliation:
    • &CNAG-CRG National Center for Genomic Analysis - Center for Genomic Regulation
    • &UPF University Pompeu Fabra
      email: marc.elosua@cnag.crg.eu
  • name: Helena L. Crowell
    affiliation:
    • &IMLS Institute for Molecular Life Sciences, University of Zurich, Switzerland
    • &SIB SIB Swiss Institute of Bioinformatics, University of Zurich, Switzerland
      email: helena.crowell@uzh.ch
      abstract: >
      <p> Spatially resolved gene expression profiles are key to understand tissue organization and function. However, novel spatial transcriptomics (ST) profiling techniques lack single-cell resolution and require a combination with single-cell RNA sequencing (scRNA-seq) information to deconvolute the spatially indexed datasets. Leveraging the strengths of both data types, we developed SPOTlight, a computational tool that enables the integration of ST with scRNA-seq data to infer the location of cell types and states within a complex tissue. SPOTlight is centered around a seeded non-negative matrix factorization (NMF) regression, initialized using cell-type marker genes and non-negative least squares (NNLS) to subsequently deconvolute ST capture locations (spots).
      package: "r BiocStyle::pkg_ver('SPOTlight')"
      vignette: >
      %\VignetteIndexEntry{"SPOTlight"}
      %\VignettePackage{SPOTlight}
      %\VignetteEncoding{UTF-8}
      %\VignetteEngine{knitr::rmarkdown}
      output:
      BiocStyle::html_document
      editor_options:
      markdown:
      wrap: 80

<style type="text/css">
.smaller {
  font-size: 10px
}
</style>

For a more detailed explanation of SPOTlight consider looking at our
manuscript:

Elosua-Bayes M, Nieto P, Mereu E, Gut I, Heyn H.
SPOTlight: seeded NMF regression to deconvolute
spatial transcriptomics spots with single-cell transcriptomes.
Nucleic Acids Res. 2021;49(9):e50. doi: 10.1093

Load packages {.unnumbered}

library(ggplot2)
library(SPOTlight)
library(SingleCellExperiment)
library(SpatialExperiment)
library(scater)
library(scran)

Introduction

What is SPOTlight?

SPOTlight is a tool that enables the deconvolution of cell types and cell type
proportions present within each capture location comprising mixtures of cells.
Originally developed for 10X's Visium - spatial transcriptomics - technology, it
can be used for all technologies returning mixtures of cells.

SPOTlight is based on learning topic profile signatures, by means of an NMFreg
model, for each cell type and finding which combination of cell types fits best
the spot we want to deconvolute. Find below a graphical abstract visually
summarizing the key steps.

[圖片上傳失敗...(image-a6d87e-1670316450053)]

Starting point

The minimal unit of data required to run SPOTlight are:

  • ST (sparse) matrix with the expression, raw or normalized, where rows = genes
    and columns = capture locations.
  • Single cell (sparse) matrix with the expression, raw or normalized,
    where rows = genes and columns = cells.
  • Vector indicating the cell identity for each column in the single cell
    expression matrix.

Data inputs can also be objects of class r Biocpkg("SpatialExperiment") (SE),
r Biocpkg("SingleCellExperiment") (SCE) or r CRANpkg("Seurat") objects from
which the minimal data will be extracted.

Getting started

Data description

For this vignette, we will use a SE put out by 10X Genomics containing a
Visium kidney slide. The raw data can be accessed
here.

SCE data comes from the The Tabula Muris
Consortium
which contains
>350,000 cells from from male and female mice belonging to six age groups,
ranging from 1 to 30 months. From this dataset we will only load the kidney
subset to map it to the Visium slide.

Loading the data

Both datasets are available through Biocondcutor packages and can be loaded into R as follows.
`
Load the spatial data:

library(TENxVisiumData)
spe <- MouseKidneyCoronal()
# Use symbols instead of Ensembl IDs as feature names
rownames(spe) <- rowData(spe)$symbol

Load the single cell data. Since our data comes from the
Tabula Muris Sensis
dataset, we can directly load the SCE object as follows:

library(TabulaMurisSenisData)
sce <- TabulaMurisSenisDroplet(tissues = "Kidney")$Kidney

Quick data exploration:

table(sce$free_annotation, sce$age)

We see how there is a good representation of all the cell types across ages
except at 24m. In order to reduce the potential noise introduced by age and
batch effects we are going to select cells all coming from the same age.

# Keep cells from 18m mice
sce <- sce[, sce$age == "18m"]
# Keep cells with clear cell type annotations
sce <- sce[, !sce$free_annotation %in% c("nan", "CD45")]

Workflow

Preprocessing

If the dataset is very large we want to downsample it to train the model, both
in of number of cells and number of genes. To do this, we want to keep a
representative amount of cells per cluster and the most biologically relevant
genes.

In the paper we show how downsampling the number of cells per cell identity to
~100 doesn't affect the performance of the model. Including >100 cells per
cell identity provides marginal improvement while greatly increasing
computational time and resources. Furthermore, restricting the gene set to the
marker genes for each cell type along with up to 3.000 highly variable genes
further optimizes performance and computational resources. You can find a more
detailed explanation in the original
paper.

Feature selection

Our first step is to get the marker genes for each cell type. We follow the
Normalization procedure as described in
OSCA. We
first carry out library size normalization to correct for cell-specific biases:

sce <- logNormCounts(sce)

Variance modelling

We aim to identify highly variable genes that drive biological heterogeneity.
By feeding these genes to the model we improve the resolution of the biological structure and reduce the technical noise.

# Get vector indicating which genes are neither ribosomal or mitochondrial
genes <- !grepl(pattern = "^Rp[l|s]|Mt", x = rownames(sce))
dec <- modelGeneVar(sce, subset.row = genes)
plot(dec$mean, dec$total, xlab = "Mean log-expression", ylab = "Variance")
curve(metadata(dec)$trend(x), col = "blue", add = TRUE)
# Get the top 3000 genes.
hvg <- getTopHVGs(dec, n = 3000)

Next we obtain the marker genes for each cell identity. You can use whichever
method you want as long as it returns a weight indicating the importance of that
gene for that cell type. Examples include avgLogFC, AUC, pct.expressed,
p-value...

colLabels(sce) <- colData(sce)$free_annotation
# Compute marker genes
mgs <- scoreMarkers(sce, subset.row = genes)

Then we want to keep only those genes that are relevant for each cell identity:

mgs_fil <- lapply(names(mgs), function(i) {
    x <- mgs[[i]]
    # Filter and keep relevant marker genes, those with AUC > 0.8
    x <- x[x$mean.AUC > 0.8, ]
    # Sort the genes from highest to lowest weight
    x <- x[order(x$mean.AUC, decreasing = TRUE), ]
    # Add gene and cluster id to the dataframe
    x$gene <- rownames(x)
    x$cluster <- i
    data.frame(x)
})
mgs_df <- do.call(rbind, mgs_fil)

Cell Downsampling

Next, we randomly select at most 100 cells per cell identity. If a cell type is
comprised of <100 cells, all the cells will be used. If we have very
biologically different cell identities (B cells vs. T cells vs. Macrophages vs.
Epithelial) we can use fewer cells since their transcriptional profiles will be
very different. In cases when we have more transcriptionally similar cell
identities we need to increase our N to capture the biological heterogeneity
between them.

In our experience we have found that for this step it is better to select the
cells from each cell type from the same batch if you have a joint dataset from
multiple runs. This will ensure that the model removes as much signal from the
batch as possible and actually learns the biological signal.

For the purpose of this vignette and to speed up the analysis, we are going to
use 20 cells per cell identity:

# split cell indices by identity
idx <- split(seq(ncol(sce)), sce$free_annotation)
# downsample to at most 20 per identity & subset
# We are using 5 here to speed up the process but set to 75-100 for your real
# life analysis
n_cells <- 5
cs_keep <- lapply(idx, function(i) {
    n <- length(i)
    if (n < n_cells)
        n_cells <- n
    sample(i, n_cells)
})
sce <- sce[, unlist(cs_keep)]

Deconvolution

You are now set to run SPOTlight to deconvolute the spots!

Briefly, here is how it works:

  1. NMF is used to factorize a matrix into two lower dimensionality matrices
    without negative elements. We first have an initial matrix V (SCE count matrix),
    which is factored into W and H. Unit variance normalization by gene is performed
    in V and in order to standardize discretized gene expression levels,
    ‘counts-umi’. Factorization is then carried out using the non-smooth NMF method,
    implemented in the R package NMF r CRANpkg("NMF") (14). This method is
    intended to return sparser results during the factorization in W and H, thus
    promoting cell-type-specific topic profile and reducing overfitting during
    training. Before running factorization, we initialize each topic, column,
    of W with the unique marker genes for each cell type with weights. In turn, each
    topic of H in SPOTlight is initialized with the corresponding belongance of each
    cell for each topic, 1 or 0. This way, we seed the model with prior information,
    thus guiding it towards a biologically relevant result. This initialization also
    aims at reducing variability and improving the consistency between runs. \

  2. NNLS regression is used to map each capture location's transcriptome in V’
    (SE count matrix) to H’ using W as the basis. We obtain a topic profile
    distribution over each capture location which we can use to determine its
    composition. \

  3. we obtain Q, cell-type specific topic profiles, from H. We select all cells
    from the same cell type and compute the median of each topic for a consensus
    cell-type-specific topic signature. We then use NNLS to find the weights of each
    cell type that best fit H’ minimizing the residuals.

You can visualize the above explanation in the following workflow scheme:

[圖片上傳失敗...(image-7c6b20-1670316450053)]

res <- SPOTlight(
    x = sce,
    y = spe,
    groups = as.character(sce$free_annotation),
    mgs = mgs_df,
    hvg = hvg,
    weight_id = "mean.AUC",
    group_id = "cluster",
    gene_id = "gene")

Alternatively you can run SPOTlight in two steps so that you can have the trained model. Having the trained model allows you to reuse with other datasets you also want to deconvolute with the same reference. This allows you to skip the training step, the most time consuming and computationally expensive.

mod_ls <- trainNMF(
    x = sce,
    y = spe,
    groups = sce$type,
    mgs = mgs,
    weight_id = "weight",
    group_id = "type",
    gene_id = "gene")
 # Run deconvolution
res <- runDeconvolution(
    x = spe,
    mod = mod_ls[["mod"]],
    ref = mod_ls[["topic"]])

Extract data from SPOTlight:

# Extract deconvolution matrix
head(mat <- res$mat)[, seq_len(3)]
# Extract NMF model fit
mod <- res$NMF

Visualization

In the next section we show how to visualize the data and interpret
SPOTlight's results.

Topic profiles

We first take a look at the Topic profiles. By setting facet = FALSE we want to
evaluate how specific each topic signature is for each cell identity.
Ideally each cell identity will have a unique topic profile associated to it as
seen below.

plotTopicProfiles(
    x = mod,
    y = sce$free_annotation,
    facet = FALSE,
    min_prop = 0.01,
    ncol = 1) +
    theme(aspect.ratio = 1)

Next we also want to ensure that all the cells from the same cell identity share
a similar topic profile since this will mean that SPOTlight has learned a
consistent signature for all the cells from the same cell identity.

plotTopicProfiles(
    x = mod,
    y = sce$free_annotation,
    facet = TRUE,
    min_prop = 0.01,
    ncol = 6)

Lastly we can take a look at which genes the model learned for each topic.
Higher values indicate that the gene is more relevant for that topic.
In the below table we can see how the top genes for Topic1 are characteristic
for B cells (i.e. Cd79a, Cd79b, Ms4a1...).

library(NMF)
sign <- basis(mod)
colnames(sign) <- paste0("Topic", seq_len(ncol(sign)))
head(sign)
# This can be dynamically visualized with DT as shown below
# DT::datatable(sign, fillContainer = TRUE, filter = "top")

Spatial Correlation Matrix

See here
for additional graphical parameters.

plotCorrelationMatrix(mat)

Co-localization

Now that we know which cell types are found within each spot we can make a graph
representing spatial interactions where cell types will have stronger edges
between them the more often we find them within the same spot.

See here for additional
graphical parameters.

plotInteractions(mat, which = "heatmap", metric = "prop")
plotInteractions(mat, which = "heatmap", metric = "jaccard")
plotInteractions(mat, which = "network")

Scatterpie

We can also visualize the cell type proportions as sections of a pie chart for
each spot. You can modify the colors as you would a standard r CRANpkg("ggplot2").

ct <- colnames(mat)
mat[mat < 0.1] <- 0
# Define color palette
# (here we use 'paletteMartin' from the 'colorBlindness' package)
paletteMartin <- c(
    "#000000", "#004949", "#009292", "#ff6db6", "#ffb6db", 
    "#490092", "#006ddb", "#b66dff", "#6db6ff", "#b6dbff", 
    "#920000", "#924900", "#db6d00", "#24ff24", "#ffff6d")
pal <- colorRampPalette(paletteMartin)(length(ct))
names(pal) <- ct
plotSpatialScatterpie(
    x = spe,
    y = mat,
    cell_types = colnames(mat),
    img = FALSE,
    scatterpie_alpha = 1,
    pie_scale = 0.4) +
    scale_fill_manual(
        values = pal,
        breaks = names(pal))

With the image underneath - we are rotating it 90 degrees counterclockwise and mirroring across the horizontal axis to show how to align if the spots don't overlay the image.

plotSpatialScatterpie(
    x = spe,
    y = mat,
    cell_types = colnames(mat),
    img = FALSE,
    scatterpie_alpha = 1,
    pie_scale = 0.4, 
    # Rotate the image 90 degrees counterclockwise
    degrees = -90,
    # Pivot the image on its x axis
    axis = "h") +
    scale_fill_manual(
        values = pal,
        breaks = names(pal))

Residuals

Lastly we can also take a look at how well the model predicted the proportions
for each spot. We do this by looking at the residuals of the sum of squares for
each spot which indicates the amount of biological signal not explained by the model.

spe$res_ss <- res[[2]][colnames(spe)]
xy <- spatialCoords(spe)
spe$x <- xy[, 1]
spe$y <- xy[, 2]
ggcells(spe, aes(x, y, color = res_ss)) +
    geom_point() +
    scale_color_viridis_c() +
    coord_fixed() +
    theme_bw()

Session information

sessionInfo()
最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
  • 序言:七十年代末蜈彼,一起剝皮案震驚了整個(gè)濱河市筑辨,隨后出現(xiàn)的幾起案子,更是在濱河造成了極大的恐慌幸逆,老刑警劉巖棍辕,帶你破解...
    沈念sama閱讀 216,324評(píng)論 6 498
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件,死亡現(xiàn)場(chǎng)離奇詭異还绘,居然都是意外死亡楚昭,警方通過查閱死者的電腦和手機(jī),發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 92,356評(píng)論 3 392
  • 文/潘曉璐 我一進(jìn)店門拍顷,熙熙樓的掌柜王于貴愁眉苦臉地迎上來抚太,“玉大人,你說我怎么就攤上這事昔案∧蚱叮” “怎么了?”我有些...
    開封第一講書人閱讀 162,328評(píng)論 0 353
  • 文/不壞的土叔 我叫張陵踏揣,是天一觀的道長庆亡。 經(jīng)常有香客問我,道長捞稿,這世上最難降的妖魔是什么又谋? 我笑而不...
    開封第一講書人閱讀 58,147評(píng)論 1 292
  • 正文 為了忘掉前任,我火速辦了婚禮娱局,結(jié)果婚禮上彰亥,老公的妹妹穿的比我還像新娘。我一直安慰自己铃辖,他們只是感情好剩愧,可當(dāng)我...
    茶點(diǎn)故事閱讀 67,160評(píng)論 6 388
  • 文/花漫 我一把揭開白布。 她就那樣靜靜地躺著娇斩,像睡著了一般仁卷。 火紅的嫁衣襯著肌膚如雪。 梳的紋絲不亂的頭發(fā)上犬第,一...
    開封第一講書人閱讀 51,115評(píng)論 1 296
  • 那天锦积,我揣著相機(jī)與錄音,去河邊找鬼歉嗓。 笑死丰介,一個(gè)胖子當(dāng)著我的面吹牛,可吹牛的內(nèi)容都是我干的。 我是一名探鬼主播哮幢,決...
    沈念sama閱讀 40,025評(píng)論 3 417
  • 文/蒼蘭香墨 我猛地睜開眼带膀,長吁一口氣:“原來是場(chǎng)噩夢(mèng)啊……” “哼!你這毒婦竟也來了橙垢?” 一聲冷哼從身側(cè)響起垛叨,我...
    開封第一講書人閱讀 38,867評(píng)論 0 274
  • 序言:老撾萬榮一對(duì)情侶失蹤,失蹤者是張志新(化名)和其女友劉穎柜某,沒想到半個(gè)月后嗽元,有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體,經(jīng)...
    沈念sama閱讀 45,307評(píng)論 1 310
  • 正文 獨(dú)居荒郊野嶺守林人離奇死亡喂击,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 37,528評(píng)論 2 332
  • 正文 我和宋清朗相戀三年剂癌,在試婚紗的時(shí)候發(fā)現(xiàn)自己被綠了。 大學(xué)時(shí)的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片翰绊。...
    茶點(diǎn)故事閱讀 39,688評(píng)論 1 348
  • 序言:一個(gè)原本活蹦亂跳的男人離奇死亡佩谷,死狀恐怖,靈堂內(nèi)的尸體忽然破棺而出辞做,到底是詐尸還是另有隱情琳要,我是刑警寧澤,帶...
    沈念sama閱讀 35,409評(píng)論 5 343
  • 正文 年R本政府宣布秤茅,位于F島的核電站稚补,受9級(jí)特大地震影響,放射性物質(zhì)發(fā)生泄漏框喳。R本人自食惡果不足惜课幕,卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 41,001評(píng)論 3 325
  • 文/蒙蒙 一、第九天 我趴在偏房一處隱蔽的房頂上張望五垮。 院中可真熱鬧乍惊,春花似錦、人聲如沸放仗。這莊子的主人今日做“春日...
    開封第一講書人閱讀 31,657評(píng)論 0 22
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽诞挨。三九已至莉撇,卻和暖如春,著一層夾襖步出監(jiān)牢的瞬間惶傻,已是汗流浹背棍郎。 一陣腳步聲響...
    開封第一講書人閱讀 32,811評(píng)論 1 268
  • 我被黑心中介騙來泰國打工, 沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留银室,地道東北人涂佃。 一個(gè)月前我還...
    沈念sama閱讀 47,685評(píng)論 2 368
  • 正文 我出身青樓励翼,卻偏偏與公主長得像,于是被迫代替她去往敵國和親辜荠。 傳聞我的和親對(duì)象是個(gè)殘疾皇子汽抚,可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 44,573評(píng)論 2 353

推薦閱讀更多精彩內(nèi)容