空转分析方法拾遗2

作者，追风少年i

周五了，一周的收官之战，马上国庆了，希望大家有个好的假期，一切都会好起来，周星驰有一句台词，虽然我们都是路人甲乙丙丁，但是一样是有生命，有灵魂的。

书接上文，空间多组学分析破译胶质母细胞瘤中的双向肿瘤-宿主相互依赖性（空间微环境），我们来进行方法拾遗并复现。

这里我们主要关注空间转录组的分析，关于空间代谢和空间蛋白，仍然需要不断的学习。

第一点，空间CNV分析。

• CNA estimation
  

  图片.png

通过将基因与其染色体位置对齐并对相对表达值应用移动平均值来估计拷贝数变异 (CNV)，每个染色体内有 100 个基因的滑动窗口（跟inferCNV的原理一致），首先，使用 InferCNV 包根据它们各自的基因组定位排列基因。 作为非恶性细胞的参考集，使用了来自非恶性皮层样本的空间转录组数据集（注意这里ref的选择）。 为了提高速度和计算能力，subsample是可选的。 为了避免任何特定基因对移动平均线的相当大的影响，将相对表达值限制在 [-2.6,2.6] 的范围内，将所有高于/低于 exp(i) = |2.6| 的值替换为 exp(i) = |2.6|。然后重新导入导出的 .RDS 输出文件，并按估计的 CNA 的染色体平均值进行分组，并使用 SPATA 对象的 fdata slot与它们的空间位置对齐。使用 SPATA::joinWithFeatures() 函数，然后提取了集群比较。空间数据的 CNA 分析是通过 SPATA2::runCnvAnalysis() 命令在 SPATA2 中实现的。在下一步中，估计的 CNA 值重新排列到定义的染色体bin中。这些 bin 由 SPATAwrappers 函数 Create.ref.bins() 创建，使用 SPATA 对象和给定 bin 的大小作为输入。在本研究中，使用了 1Mbp 的 bin 大小，产生 3847 个染色体 bin，每个 bin 平均覆盖 5.5 个基因。使用 10kbp 滑动窗口进行重新缩放和插值。对于标准化，使用了一个 loess 回归模型，用于确定来自 InferCNV 输出的拷贝数值。使用 SPATAwrappers::run-CNV.Normalization() 函数进行Interpolation和归一化。

SPATA分析CNV的代码如下

library(SPATA2)
spata_obj <- loadSpataObject("data/spata-obj-gbm275.RDS")###加载infercnv分析得到的rds

spata_obj <-
  runCnvAnalysis(
    object = spata_obj
    directory_cnv_folder = "data-gmb275/cnv-results" # example directory
    cnv_prefix = "Chr"
    )

# change input
spata_obj <- 
  runCnvAnalysis(
    object = spata_obj
    directory_cnv_folder = "data-gmb275/cnv-results" # example directory to 
    clear_noise_via_ref_mean_sd = list(sd_amplifier = 2)

# example on how to use the functions arguments if you want to specify the reference (ref_*)
runCnvAnalysis(
  object = spata_obj,
  directory_cnv_folder = "data-gmb275/cnv-results", # example directory to 
  ref_annotation = cnv_ref$annotation, 
  ref_mtr = cnv_ref$mtr, 
  ref_regions = cnv_ref$regions
)

##CNV-Results
cnv_results <- getCnvResults(spata_obj)

names(cnv_results)
##Visualization
plotCnvResults(object = spata_obj)

图片.png

plotSurface(object = spata_obj, color_by = "seurat_clusters")

plotCnvResults(object = spata_obj, across = "seurat_clusters")

图片.png

图片.png

###CNV-Heatmap
getCnvFeatureNames(object = spata_obj)

# are part of all feature names
getFeatureNames(object = spata_obj)
##            numeric            integer            numeric            numeric 
##   "nCount_Spatial" "nFeature_Spatial"       "percent.mt"       "percent.RB" 
##             factor          character             factor             factor 
##  "seurat_clusters"     "segmentation"      "Leiden_UMAP"          "cluster" 
##            logical            numeric            numeric            numeric 
##             "keep"             "Chr0"             "Chr1"            "Chr10" 
##            numeric            numeric            numeric            numeric 
##            "Chr11"            "Chr12"            "Chr13"            "Chr14" 
##            numeric            numeric            numeric            numeric 
##            "Chr15"            "Chr16"            "Chr17"            "Chr18" 
##            numeric            numeric            numeric            numeric 
##            "Chr19"             "Chr2"            "Chr20"            "Chr21" 
##            numeric            numeric            numeric            numeric 
##            "Chr22"            "Chr23"             "Chr3"             "Chr4" 
##            numeric            numeric            numeric            numeric 
##             "Chr5"             "Chr6"             "Chr7"             "Chr8" 
##            numeric 
##             "Chr9"
plotSurface(object = spata_obj, color_by = "Chr7", pt_clrsp = "Reds 3", c1 = 1)

plotSurface(object = spata_obj, color_by = "Chr10", pt_clrsp = "Blues 3", rev = FALSE, c1 = 1)

图片.png

CNA subclone analysis
为了根据 CNA 数据确定亚克隆，进行了 PCA 分析，然后对前 30 个组分进行 KMeans 聚类（CNV得到的矩阵进行分析）。为了确定最优的tree cut，我们估计了 Calinski-Harabasz 指数并使用了 k~CH 指数曲线的第一个峰值（CH 指数 vs Number of clusters）。 下一个目标是量化空间上不同的转录程序与亚克隆结构的空间相关性和重叠。 首先，我们通过拟合值 (f(c_k)) 估计给定 CNA 亚克隆的每个点 (S1, S2;.; Sk) 到 KMeans 质心 (c) 的距离 (Sdist)，如下所示：

图片.png

从进一步分析中去除了距离较近的点（Sdist < quantilðSdist ;0.9）。 然后，估计了所有给定点的主要空间不同转录程序，并量化了给定亚克隆中每个转录程序的表达百分比。

第二部分

• Prediction of tumor cell content
  为单细胞测序研究建立的分析和算法不一定适用于空间分辨转录组数据的分析。 这主要是由于每个点内的细胞组成不明确，分析这些数据的一个优势是近似每个点的真实组成。单细胞数据和空间转录组学的整合已经可以预测某些细胞类型的空间丰度，但在肿瘤中，由于样本内部和样本之间存在巨大的异质性，这种预测很困难。 与高度可变的基因表达相比，染色体改变是评估每个点的肿瘤细胞含量的相当稳健和合适的参数。 为了训练和验证预测模型，对同一供体进行了单细胞测序和相应的空间分辨转录组学。

推断了 scRNA-seq 数据集的拷贝数变化，并提取了肿瘤和非恶性细胞的特征，这些特征显示出很大的相似性。

Validation of CNA Analysis

接下来，验证了从 RNA-seq 数据调用的推断 CNA 反映 DNA 测序中检测到的 CNA 的程度。 因此，使用了最近发表的 slide-RNA-seq 和 slide-DNA-seq 数据集，并比较了两种空间测序方法检测到的 CNA。 两种验证方法都表明，从空间转录组学中获得的推断 CNA 反映了使用 10X Visium 平台获得的空间分辨转录组学中的实际染色体增益和损失。

第三点，Spatial autocorrelation Moran’s I，关于这个内容，用过monocle3的童鞋应该很熟悉。

空间自相关使用 Moran 的 I 统计值根据特征位置和属性值确定点与其邻居的空间依赖性。 这允许评估定义的基因表达或特征是分组的、分布的还是随机的。 如果 Z 值或 p 值表示统计显著性，则 Moran's I 指数值为正表示有聚集趋势，而 Moran's I 指数值为负表示有分散趋势。 自相关定义为：

图片.png

where N defines the number of spatial spots with the index x and y and the matrix of spatial weights as w_xy. The feature is indicated by exp and the mean value of neighboring features. The spatial weights were defined as a distance matrix of the cartesian space.W is defined as the sum of all w_xy.

最后来一张空间轨迹

空间轨迹

主要还是针对CNV的分析为主

最后我把下面的内容和视频中细胞定义的脚本（clustermole+singleR+cellID）列在下面，有需要的可以拿去研究。

单细胞数据分析之单样本分析篇
单细胞数据分析之多样本整合篇
单细胞数据分析之细胞注释篇
单细胞数据分析之轨迹分析篇
单细胞数据分析之细胞通讯篇
单细胞数据分析之CNV篇
单细胞数据分析之转录因子篇
空间数据的全部个性化分析PPT