10X空间转录组CNV和速率分析的空间图谱

hello，我们今天要继续扩展10X空间转录组的个性化分析，其中主要是CNV分析和速率分析的空间图谱，参考文献在Spatiotemporal heterogeneity of glioblastoma is dictated by microenvironmental interference，大家可以好好看一看，许多方法值得我们借鉴和运用。

首先来分享空间转录组的CNV分析，关键是扩展到空间的CNV图谱，其实关于空间转录组的CNV分析，大家可以参考我之前分享的文章10X空间转录组数据分析之CNV事件的spatial landscape、10X空间转录组数据分析之思路总结（针对肿瘤样本）、10X空间转录组数据推断基因拷贝数畸变（copy number aberrations）CNA等，而我们今天分享的CNV，是之前的延申。

inferCNV的分析会得到很多的文件，其中这里我们关注三个文件

（1）17_HMM_predHMMi3.hmm_mode-samples.pred_cnv_regions.dat，染色体区段的CNV事件

（2）17_HMM_predHMMi3.hmm_mode-samples.pred_cnv_genes.dat，针对基因的CNV事件。

（3）17_HMM_predHMMi6.rand_trees.hmm_mode-subclusters.cell_groupings，针对每个细胞的CNV标识

根据上述的三个文件，就可以拿到样本全面的CNV信息，从而构建CNV图谱。

看一看CNV分析的详细做法

通过将基因与其染色体位置对齐并对相对表达值应用移动平均值来估计拷贝数变异 (CNV)，每个染色体内有 100 个基因的滑动窗口。首先，使用 InferCNV包根据基因组各自的定位排列基因。作为非恶性细胞的参考集，我们使用了来自非恶性皮质样本的空间转录组数据集(这里采用空间转录组数据直接进行CNV分析)。为了提高速度和计算能力，下采样是可选的。为了避免任何特定基因对移动平均值的显着影响，我们通过使用 infercnv 包替换所有高于/低于 exp(i)=|2.6| 的值来限制相对表达值 [-2.6,2.6]。 如前所述，这仅在 CNV 估计的背景下进行。 导出的 .RDS 文件被重新导入并按估计的 CNV 改变的染色体平均值分组，并使用 SPATA对象的 fdata slot对齐到它们的空间位置。 使用 SPATA::joinWithFeatures() 函数进行聚类比较的提取。 此外，还实现了选择染色体中变化最大的基因的选项。

最后挑选显著的CNV事件，实现我们的空间CNV图谱

接下来我们第二个分析点，空间转录组的速率分析。

其实在空间转录组出现的时候，对于轨迹的分析就希望得到细胞不仅仅是时间的变化，更多的是想要空间上的迁移，大家可以参考我之前分享的内容10X空间转录组的轨迹分析，但是当时限于空间转录组的分析精度，很少再深入研究，不过随着分析的进展，单细胞空间联合方法、10X空间转录组速率分析（Velocyto）之SIRV等等，越来越趋向于空间转录组对于轨迹的判断，虽然很多时候细胞群之间不一定有分化关系，但是这样的转变确实体现了基因、通路的迁移。

空间转录组做RNA速率之前呢，大家最好先看看之前分享的内容，包括10X单细胞（10X空间转录组）基于马尔可夫过程和RNA速率的细胞命运预测（Cellrank代码分享）、10X单细胞（10X空间转录组）基于马尔可夫过程和RNA速率的细胞命运预测（Cellrank）、10X空间转录组数据分析之Pattern recognition and clustering，因为空间转录组的分析是建立在这里基础之上的。

其中空间上的坐标是一致的，也可以借用单细胞的方法，在UMAP图上展示轨迹，但是我最欣赏下面的降维展示方法：

关于PCA的意义，这里不再多说，大家可以参考我之前的文章10X单细胞10X空间转录组降维分析之PCA轴的秘密，每个轴都有特定的含义，这样展示的方式更加一目了然。

接下来具体的做法RNA velocity estimation

使用 CellRanger BAM 文件通过ready-to-use pipeline从 velocyto 包中分离拼接和未拼接读数的表达矩阵。生成的 .loom 文件被读入 scVelo Seurat。合并Seurat 对象并执行了批次效应的去除。数据整合后，具有外显子和内含子基因级 UMI 计数的 Seurat 对象被转换为 h5ad 格式。我们将 h5ad 文件读入到 AnnData 对象中。接下来，执行标准化并通过scVelo包选择了 2,000 个高变的基因。排除了外显子和内含子read中分配读数少于 20 的所有基因，并使用动力学模型估计了 RNA 速度和潜伏时间。数据将导出为 .csv 文件并实施到 SPATA 对象中以进行进一步的可视化。The explained pipeline is implemented into a SPATA wrapper for scVelo (SPATA::getRNA velocity, in the development branch).

可视化的具体做法，通过使用集成到 adata.obsm['X_umap'] slot中的 AnnData 对象中的前 3 个主成分 (PC1-3) 对所有肿瘤样本进行可视化。 velocity streams由参考“X_umap” slot的 pl.velocity_embedding_stream 函数计算。 在空间转录组数据中，目标是在添加velocity streams并保留空间结构。 我们将空间坐标从 SPATA对象迁移到 AnnData 对象到用于 pl.velocity_embedding_stream 函数的 adata.obsm['X_umap'] slot中。

Estimation of transient gene expression programs along RNA velocity streams，为了估计在空间（空间转录组学）和时间（RNA 速度估计）中动态调节的转录程序，使用计算出的velocity streams作为空间轨迹。 使用 SPATA::createTrajectories 函数，沿着最近描述的时空轨迹寻找遵循预定义动态的基因。

不过我还是很喜欢下面的展示方式

• 注：Diffusion plots (dimensional reduction, with RNA-velocity) indicate the difference of reactive and lineage differentiation along the major axis. Heatmaps of the lineage to reactive trajectory are illustrated. On the bottom of the heatmap, a barplot indicates the sample source. A lineplot at the top of the heatmap illustrates the gene set enrichment analysis of the reactive immune and NPC-like cell states

最后补充一个Spatial correlation analysis

为了绘制空间相关基因表达或基因集富集，使用了 z score的归一化表达值。 基因集的一个基因表达vector或富集vector用于沿空间轨迹对spot进行排序。 构建从最低排名到最高排名的点的轨迹（基于 z score输入vector）。 感兴趣的基因（与空间轨迹相关）通过来自 stats-package（R 软件）的 loess-fit 进行拟合，并与排名spot对齐并进行缩放。 通过 Pearson 的乘积矩相关系数进行相关分析。 对于热图说明，通过沿预定义的空间轨迹对黄土拟合表达的最大峰值进行排序来计算基因顺序。

好了，不知道大家学到了么

生活很好，有你更好