大家好紊搪，這一次我們來分享一下有關(guān)做10X空間轉(zhuǎn)錄組的分析軟件SPATA，其實(shí)關(guān)于10X空間轉(zhuǎn)錄組做軌跡分析的文章之前分享過一個(gè)全景，文章在10X空間轉(zhuǎn)錄組的軌跡分析耀石，其中用到的軟件是stlearn，我們放一張stlearn分析出來的圖

圖片.png

今天我們來分享另外一個(gè)軟件爸黄，SPATA滞伟，文章在Inferring spatially transient gene expression pattern from spatial transcriptomic studies揭鳞，我們還是先分享文獻(xiàn)，最后分享示例代碼诗良。

Abstract

Spatial transcriptomic is a technology to provide deep transcriptomic profiling by preserving the spatial organization.（這個(gè)老生常談了）汹桦。Here, we present a framework for SPAtial Transcriptomic Analysis (SPATA, https://themilolab.github.io/SPATA), to provide a comprehensive characterization of spatially resolved gene expression, regional adaptation of transcriptional programs and transient dynamics along spatial trajectories.（運(yùn)用這個(gè)軟件來做空間軌跡分析鲁驶，講道理鉴裹，我很好奇！钥弯！）径荔。

Brief Communication（這個(gè)地方我們總結(jié)一下）

Deep transcriptional profiling of single cells by RNA-sequencing maps the cellular composition of tissue specimens regarding cellular origin, developmental trajectories and transcriptional programs（這是單細(xì)胞的優(yōu)勢），However, information determining the spatial arrangement of specific cell types or transcriptional programs are lacking and thus can only be predicted indirectly（空間信息的缺失）脆霎，which is a considerable drawback of this method.（也是眾所周知的缺點(diǎn)）总处。空間組織結(jié)構(gòu)傳統(tǒng)的方式是用成像技術(shù)提供高精度的信息，但是同時(shí)檢測的基因和蛋白的數(shù)量非常少睛蛛，當(dāng)然鹦马，已經(jīng)有了一些新的空間組織檢測技術(shù)，但受到空間分布細(xì)胞精度和深度的制約忆肾，Further, data integration, visualization and analysis of transcriptomic and spatial information remains challenging. Here, we present a software tool to provide a framework for integration of high-dimensional transcriptional data within a spatial context.這個(gè)軟件集友好的可視化界面荸频，分割或軌跡分析，提供了廣泛了不同分析方法客冈，In particular, we focus on transient changes of gene expression and aim to infer transcriptional programs that are dynamically regulated as a function of spatial organization.（這個(gè)軟件就是專門為10X空間轉(zhuǎn)錄組做準(zhǔn)備的）旭从。

In order to present an overview of possible analytic capabilities of the SPATA workflow,

圖片.png

• In a first step, we combine shared-nearest neighbor clustering and spatial pattern recognition by an external tool（SPARK）in order to determine genes with a defined spatially resolved expression pattern.（這個(gè)地方類似于Seurat計(jì)算空間差異基因）场仲。We found that the cortical layering is accurately reflected by our clustering approach（這個(gè)很自然）和悦。In order to gain insights into the spatial organization we provided a tool to compute the spatial distance within the defined layers or correspondent clusters. An increasing distance within individual clusters allows to differentiate between narrowly related or a widespread dispersion of spots within the cluster.(計(jì)算空間距離，這個(gè)在我分享的文章空間轉(zhuǎn)錄組細(xì)胞類型的距離分析之二---代碼實(shí)現(xiàn)渠缕，提到過,只不過注意鸽素，這里計(jì)算的是cluster內(nèi)部的距離)。
• Next, the spatial overlap of transcriptional programs or gene expression was analyzed using a Bayesian approach（有關(guān)貝葉斯的知識大家可以參考我的文章貝葉斯分類器(10X單細(xì)胞和10X空間轉(zhuǎn)錄組的基礎(chǔ)算法)）,resulting in an estimated correlation which quantifies the identical arrangement of expression in space.
• In a further step, we aimed to analyze dynamic changes, which were annotated using pseudotime estimation or RNA-velocity.(軌跡推斷)亦鳞。We directly implemented the pseudotime inference from the monocle3 package, but also allow the integration of any other tool such as “l(fā)atent-time” extracted from RNA-velocity (scVelo). （看來這個(gè)軟件的軌跡推斷需要借助單細(xì)胞的軌跡分析軟件）馍忽。Another option for dynamic gene expression analysis is the detection of defined transcriptional programs along a defined trajectory.（基因的動態(tài)表達(dá)）。我們在方法中重點(diǎn)關(guān)注一下這個(gè)蚜迅。
• Moreover, we provide the opportunity to screen for gene expression or transcriptional programs which transiently change along predefined trajectories by modelling gene expression changes in accordance to various biologically relevant behaviors.（軌跡分析的轉(zhuǎn)變基因）舵匾。
  
  圖片.png
  
  
  圖片.png
  
  這個(gè)地方我們稍微討論一下
• 第一步，SPARK識別空間pattern谁不，這個(gè)pattern需要人為指定坐梯，看來這個(gè)劃分了3個(gè)pattern
• 第二步，各個(gè)pattern在空間上的映射刹帕，這個(gè)地方主要對各個(gè)pattern做一下富集吵血，看看區(qū)域分布
• 第三步谎替，Copy-number analysis of the sample and cluster annotation（空間注釋）
• 第四步，Comparison of recognized pattern with known gene expression classification （識別模式與已知基因表達(dá)分類的比較 蹋辅，看看是否一致）
• 第五步钱贯，模式中顯著表達(dá)的途徑的表達(dá)（就是顯著富集的通路）
• 第六步，腫瘤區(qū)域的空間軌跡分析侦另。
• 第七步秩命，也是最為重要的一步，Change of z-scored geneset expression along the trajectory (i) and marker genes of microenvironmental alterations and inflammation（關(guān)鍵基因的變化）褒傅。

看來對空間上的動態(tài)變化基因還是要引起重視弃锐。不過這里空間SPOT做inferCNV有待商榷。

還有一個(gè)很有意思的地方殿托，Thus, we were able to show 116 that immune related genes from myeloid cells and reactive astrocytes were localized in a “glial-scars” resembling structure, sharply separating normal brain from tumor regions Figure 2h. We observed a transient increase of macrophage and microglia activation directed towards the tumor boarder,這個(gè)方法是不是值得一試呢霹菊？。

Methods支竹，我們關(guān)注一下其中的重點(diǎn)

Modeling of transient gene expression along spatial trajectories

圖片.png

Spatial distance measurement

圖片.png

Spatial overlap and correlation analysis

圖片.png

我們來分享一下示例代碼

這里的示例代碼很多旋廷，我們關(guān)注一下軌跡部分的代碼

Vsiualization
Spatial trajectories
Segmentation
DE - Analysis

空間軌跡推斷

With spatial trajectory analysis SPATA introduces a new approach to find, analyze and visualize differently expressed genes and gene-sets in a spatial context. While the classic differential gene expression analyzes this differences between experimental groups as a whole it neglects changes of expression levels that can only be seen while maintaining the spatial dimensions. Spatial trajectories allow to answer questions that include such a spatial component. E.g.:

1. In how far do expression levels change the more we move towards a region of interest?
2. Which genes follow the same pattern along these paths?

The spatial trajectory tools provided in SPATA enable new ways of visualization as well as new possibilities to screen for genes, gene-sets and patterns of interest.

load package

[library](https://rdrr.io/r/base/library.html)(SPATA)
[library](https://rdrr.io/r/base/library.html)(magrittr)
[library](https://rdrr.io/r/base/library.html)(ggplot2)
[library](https://rdrr.io/r/base/library.html)(patchwork)

load spata-object

spata_obj <- [loadSpataObject](https://themilolab.github.io/SPATA/reference/loadSpataObject.html)(input_path = "data/spata-obj-trajectory-analysis-tumor-example.RDS")

Draw trajectories

Spatial trajectories of a sample in a given spata-object can be drawn interactively using the function createTrajectories() as shown in the example below. createTrajectories() opens a mini-shiny application. This app allows one the one hand to investigate the sample with regards to spatial gene expression like plotSurfaceInteractive() does and on the other hand to draw trajectories through the areas of interest in four easy steps.

library(SPATA)
# open interactive application
spata_obj <- [createTrajectories](https://themilolab.github.io/SPATA/reference/createTrajectories.html)(object = spata_obj)

1. Drawing

Draw the trajectory interactively by clicking on the locations where you want the start and endpoints to be. If you click more than two times the trajectory is divided into sub-trajectory parts which will be displayed by subsequent plotting functions. Hence, it makes sense to stop everytime where you would expect a change as this will be highlighted by the transition from part to the subsequent one.

Figure 1.1. Interactive trajectory drawing

Figure 1.1. Interactive trajectory drawing

2. Highlighting

If start- and end-points are set adjust the trajectory width to and click on ‘Highlight’.

Figure 1.2. Highlight the trajectory

Figure 1.2. Highlight the trajectory

3. Naming

If you are satisfied with your trajectory enter a valid trajectory name and - if you want - a comment.

4. Saving

Click on ‘Save’ and on ‘Close application’ in order to return the updated spata-object.

Trajectories are stored inside the object itself. Therefore createTrajectories() returns an updated version of the object that was specified as the function’s input now containing the additional trajectories.

To obtain all trajectory names that have been drawn and saved in your spata-object use getTrajectoryNames(). Depending on the number of samples specified it returns a character vector or a list of such.

get trajectory names

getTrajectoryNames(object = spata_obj)

## [1] "tumor-layers"

Visualize trajectories

In order to visualize specific trajectories make use of plotTrajectory().

plotTrajectory(object = spata_obj,
               trajectory_name = "tumor-layers",
               color_to = "seurat_clusters",
               pt_size = 1.7,
               pt_clrp = "npg",
               display_image = FALSE) +
               theme(legend.position = "top") +

plotTrajectory(object = spata_obj,
               trajectory_name = "tumor-layers",
               color_to = "pseudotime",
               pt_size = 1.7,
               smooth = TRUE,
               smooth_span = 0.0125,
               display_image = FALSE) +
               theme(legend.position = "top")

Figure 2.1. Surface related visualization of the trajectory

Figure 2.1. Surface related visualization of the trajectory

Visualize trajectory trends

Trajectory trends via lineplots

To display the trend of certain variables along a trajectory via lineplots use one of the following three functions:

• `plotTrajectoryFeatures()
• `plotTrajectoryGeneSets()
• `plotTrajectoryGenes()

A call to one of those three functions will return a smoothed lineplot where the direction of the trajectory is mapped onto the x-axis and the values of the chosen aspects are mapped onto the y-axis. All functions follow the same concept and share the majority of their arguments. Still, they feature slight differences given the nature of their aspects to display.

The following three arguments have to be set in order to specify the trajectory that is to be analyzed:

• object
• trajectory_name
• of_sample (if more than one sample exists in the object, else it can be neglected)

The way values are smoothed and displayed can be manipulated via:

• smooth_method
• smooth_span
• smooth_se

Features along trajectories

Specify the features of interest via the argument features. Note that the function will compute the mean of all points on the axis orthogonal to the trajectory’s direction. As only numeric variables allow for such computation displaying the trend of categorical variables (e.g. cluster) is not possible via lineplots.

plotTrajectoryFeatures(object = spata_obj,
                       trajectory_name = "tumor-layers",
                       features = "pseudotime",
                       smooth_method = "loess",
                       smooth_span = 0.2,
                       smooth_se = TRUE)</pre>

圖片.png

Figure 3.1. Visualization of ontinuous features along a trajectory

In order to plot discrete features use `plotTrajectoryFeaturesDiscrete().

plotTrajectoryFeaturesDiscrete(object = spata_obj,
                               trajectory_name = "tumor-layers",
                               discrete_feature = "seurat_clusters",
                               clrp = "npg",
                               display_trajectory_parts = FALSE)

Figure 3.2. Visualization of discrete features along a trajectory

Figure 3.2. Visualization of discrete features along a trajectory

Genes and gene-sets along trajectories

Specify the genes or gene-sets of interest by using the respective functions and arguments.

# gene-set names
genes_of_interest <- [c](https://rdrr.io/r/base/c.html)("CALM1", "VIM", "GFAP")

# plot lineplot
plotTrajectoryGenes(object = spata_obj,
                    trajectory_name = "tumor-layers",
                    genes = genes_of_interest,
                    smooth_span = 0.2,
                    smooth_se = TRUE)

Figure 3.3. Visualization of gene expression along a trajectory

Figure 3.3. Visualization of gene expression along a trajectory

plotTrajectoryGeneSets(
  object = spata_obj,
  trajectory_name = "tumor-layers",
  gene_sets = "HM_HYPOXIA",
  display_trajectory_parts = FALSE) + # results in missing vertical lines 
  ggplot2::theme(legend.position = "top")

image

Analyze & model trajectories

To find variables that follow a certain trend along the trajectory of interest use assessTrajectoryTrends(). It fits the trajectory-trends of all specified variables to a variety of patterns and assesses their fit by calculating the area under the curve of the residuals.

all_genes <- getGenes(spata_obj)
all_gene_sets <- getGeneSets(spata_obj)

# obtain an assessed trajectory data.frame for all genes
atdf_gs <- assessTrajectoryTrends(object = spata_obj,
                                  trajectory_name = "tumor-layers",
                                  variables = all_genes)

# obtain an assessed trajectory data.frame for all gene-sets
atdf_gene_sets <- assessTrajectoryTrends(object = spata_obj,
                                         trajectory_name = "tumor-layers",
                                         variables = all_gene-sets)

# output example
atdf_gs

圖片.png

The resulting data.frame contains the evaluation of all fits for all provided variables. plotTrajectoryFit() visualizes the fitting results for a particular variable.

# compare the trend of a variable to different models
plotTrajectoryFit(object = spata_obj,
                  trajectory_name = "tumor-layers",
                  variable = "GFAP",
                  display_residuals = TRUE) +
  ggplot2::theme(legend.position = "top")

Figure 4.1. Visulization of the fitting results

Figure 4.1. Visulization of the fitting results

To get a summary of your trajectory assessment use examineTrajectoryAssessment() which displays the distribution of all areas under the curve measured. Use filterTrend() or a combination of dplyr::filter() and dplyr::pull() on your assessed trajectory data.frame in order to extract the variables following a specific trend to a certain quality-degree.

# example 1: extract all variables that follow the linear descending trend while moving towards the hypoxic area (See Figure 2.1)
# with an auc-evaluation equal to or lower than 2
descending_gene_sets <-
  filterTrends(atdf = atdf_gene_sets,
               limit = 2,
               trends = [c](https://rdrr.io/r/base/c.html)("Linear descending"),
               variables_only = FALSE) # return a data.frame

descending_gene_sets

圖片.png

# example 2: extract variables that featured a peak along the trajectory 
# with auc-evaluation equal to or lower than 4
peaking_genes <- filterTrends(atdf = atdf_gs,
                              limit = 4,
                              trends = c("Early peak", "Late peak", "One peak"),
                              variables_only = TRUE) # return a vector of variables

head(peaking_genes)

## [1] "CKB"    "METRN"  "TMSB4X" "LPL"    "CALM1"  "ITGA3"

Trajectory trends via heatmaps

Displaying the expression trends of too many variables in a lineplot gets unclear quite rapidly. In order to deal with several variables use plotTrajectoryHeatmap(). While the trajectory’s directory is still displayed from left to right the dynamic of the gene-sets’ expression values is now displayed by the gradual change of a row’s color instead of a line’s slope in a lineplot. The advantage of using heatmaps to display dynamic expression changes along trajectories becomes obvious when the number of variables to display rises.

descending_gene_sets <- dplyr::pull(descending_gene_sets, variables)

plotTrajectoryHeatmap(object = spata_obj,
                      trajectory_name = "tumor-layers",
                      variables = descending_gene_sets,
                      arrange_rows = "maxima",
                      show_row_names = FALSE, # gene set names are long...
                      split_columns = TRUE,
                      smooth_span = 0.5)

Figure 4.2) A heatmap displaying linear descending trends

Figure 4.2) A heatmap displaying linear descending trends

plotTrajectoryHeatmap(object = spata_obj,
                      trajectory_name = "tumor-layers",
                      variables = peaking_genes,
                      arrange_rows = "maxima",
                      show_row_names = FALSE, # gene set names are long...
                      split_columns = TRUE, # splits the heatmap 
                      smooth_span = 0.5)

Figure 4.3. A heatmap displaying the gradual change of peaking genes

Figure 4.3. A heatmap displaying the gradual change of peaking genes

這個(gè)方法光看示例好像沒什么優(yōu)勢，需要在實(shí)際中檢驗(yàn)一下了礼搁，感覺整體效果不如stlearn

圖片.png

生活很好饶碘，有你更好