单细胞空间转录分析之Scanpy

Posted 2023-03-21

参考技术A

将空间位置信息和转录组分析相结合，对于癌症、免疫、肿瘤免疫相互作用，组织微环境，神经和发育等领域，有着令人期待的应用前景。

单细胞的一切分析，加前缀Spatial 都是一个新的分析点，因此Scanpy 扩展后也可用于空间转录组数据分析。 https://scanpy-tutorials.readthedocs.io/en/latest/spatial/integration-scanorama.html

python 包Scanpy很多函数是借鉴了R包Seurat，所以这两种方法分析结果差异不大，大家可以对照Seurat分析，上面网址也提供了Seurat包处理单细胞空间转录分析过程。

和分析单细胞转录组数据一样，单细胞空间转录组主要包括了：质控（QC），标准化（Normalization），降维聚类（Dimensional reduction and clustering），Cluster marker genes， Spatially variable genes

为了和Seurat结果比较，我们使用了相同的一套数据集： https://support.10xgenomics.com/spatial-gene-expression/datasets/1.1.0/V1_Mouse_Brain_Sagittal_Anterior ， 新鲜的冷冻小鼠脑组织， 前牙矢状切面，可以参考前面讲述的ABA大脑图谱： https://www.jianshu.com/p/5d087fffeb35
导入相关包

读取数据

预处理

我们根据总counts和表达的genes对spots进行一些基本的过滤:

这儿和Seurat得到的QC小提琴图一样，只是形式不同。

标准化，HVG

PCA，聚类，可视化

每一簇marker genes

空间特异性genes

保存数据

我们可以发现与Seurat相比，分类结果还是有差异，不过大的区域识别两种方法都没有什么问题。

10X单细胞（10X空间转录组）多样本批次效应去除分析之RCA2

参考技术A

Data processing can also be carried out with Seurat. Here is an example how you can combine a RCA analysis with data preprocessed in Seurat.

Using the same 10x data as before, we generate a Seurat object and perform an initial analysis:

To run RCA, no further processing steps would be needed. However, we want to also compare the RCA result to the Seurat based clustering, therefore we first go on with a Seurat based analysis:

Based on the Elbowplot (not shown here), we use 20 PCs for further analysis.

We generate a UMAP of the data stored in the Seurat object using the umap R package:

We use the RCA function createRCAObject to generate a RCA object from the raw and optionally also the normalized data stored in our Seurat object.

Next, we can compute the projection, cluster the data, and estimate the most likely cell type for each cell as above:

Using the RCA cell type labels, RCA and Seurat clusters, we generate two new UMAPs whose coordinates are based on the PCs derived from HVGs and that are colored according to RCA clusters and cell type labels.

The RCA clusters show a high concordance to the Seurat clusters shown in the previous UMAP.

For greater convenience the results of RCA can be saved within the Seurat object for further analysis.

Also, a UMAP reduction based on the projection space can be added to the Seurat object:

RNA velocity describes the rate of gene expression change for an individual gene at a given time point based on the ratio of its spliced and unspliced messenger RNA (mRNA). Here, we describe how one can use the scvelo package, in Python, to visualize RNA velocity on the RCA generated result.

To transfer spliced RNA counts to scvelo, first transpose the raw RCA data matrix to get a cells x genes matrix, and export it to a CSV file.

In addition, export the RCA projection and UMAP embeddings to respective CSV files too.

Create an iPython notebook in the same folder and import the required packages as below.

Then, create a Scanpy object using the raw counts from the CSV file.

Populate the PCA slot in the Scanpy object as the projection data from RCA.

Populate the UMAP slot in the Scanpy object as the umap coordinates from RCA.

Load the unspliced loom object generated by velocyto .

Then, merge the spliced and unspliced objects together as described below:

As recommended by the scvelo tutorial, perform the following steps to compute RNA velocity:

It is possible that not all barcodes had sufficient quality of both spliced and unspliced reads, and thus some cells may have been discarded during the merging process. To ensure your cell type labels are still maintained, export the merged data observations from the merged scvelo object to a CSV file.

In R, load this CSV file in and extract the RCA labels and filter only those which were considered in the merged data by scvelo.

Note: If your cell names have underscores in them, scanpy will automatically split the cell name into barcode and sample_batch.

In this case, replace the last line of the above block of code with the following:

Now export these cluster labels to a CSV file.

Back in the scvelo iPynb, load this RCA cluster annotation table and set it as the observation slot of your merged data.