阿狗工具

作者，追风少年i~

国庆前的最后一弹，分享一个简单的内容，空间轨迹向量场。

其中关于空间轨迹，我也写了很多，文章放在下面，供大家参考

• 时空轨迹分析导论
• 空间转录组之空间基因和细胞轨迹
• 单细胞个性化分析之轨迹分析篇

图片.png

首先我们来解读以下这个图片，这个地方类似于基因、细胞类型或者通路的区域转换（细胞迁移）。为了探索代谢改变区域中迁移基因表达特征的富集，确定了特定基因表达特征的低富集和高富集之间的定向梯度的空间方向。 简化后，每个点的方向向量是基于其局部邻域中所研究的基因表达特征的分级富集。这些向量场计算使我们能够近似空间基因表达轨迹，从而能够识别空间上相反的转录途径。基于这些矢量场计算，报告缺氧响应和迁移特征显示反向空间轨迹（上图C、D）。 总之，研究结果为代谢变化和氧化应激是基因组多样性的潜在互惠驱动因素提供了证据，从而导致 GBM 中的克隆进化。

其中我们要实现的部分在

空间轨迹向量场

话不多说，我们直接来

library(ggplot2)
library(Seurat)
library(SPATA2)
library(dplyr)

source('runVectorFields.R')
source('plotVectorFields.R')

####脚本放在最后

x = readRDS(Seurat.spatial.rds)
y = transformSeuratToSpata(seurat_object = x, sample_name = 'FT')

####这个时候我们把需要分析的轨迹基因、通路或者细胞类型进行分析，这里以CD3D为例

data = runVectorFields(y,'CD3D')

head(data)

data

有了data，就可以进行向量场绘图

p = plotVectorFields(data,'CD3D')
pdf('eg.pdf')
print(p)
dev.off()

就会得到向量场图。

图片.png

其中的颜色，点的大小都可以更改，选择自己喜欢的搭配，当然了，我这里是拿一个基因作为展示，更为有生物学意义的是细胞类型和信号通路，照猫画虎就可以了（就把对应一个的基因值替换成你想要的细胞类型分数或者通路得分）。

附上source的代码

• runVectorFields.R

runVectorFields <- function(object,
                            features,
                            cut_off=NULL,
                            normalize=T,
                            smooth=T,
                            smooth_span=NULL,
                            dist.spot=10,
                            run.mcor=T,
                            workers=8,
                            ram=50){


  # Get the data:
  df <- SPATA2::hlpr_join_with_aes(object,
                                   df = SPATA2::getCoordsDf(object),
                                   color_by = features,
                                   normalize = normalize,
                                   smooth = smooth,
                                   smooth_span=smooth_span)

  if(!is.null(cut_off)){df[df[,features]<cut_off, features]=0}


  # Prepare data ------------------------------------------------------------
  NN.file <- SPATAwrappers::getSurroundedSpots(object)


  if(run.mcor==T){
    base::options(future.fork.enable = TRUE)
    future::plan("multisession", workers = workers)
    future::supportsMulticore()
    base::options(future.globals.maxSize = ram *100* 1024^2)
    message("... Run multicore ... ")

  }

  if(dist.spot<min(NN.file %>% filter(distance!=0) %>% pull(distance))){
    dist.spot <-min(NN.file %>% filter(distance!=0) %>% pull(distance))
    message(paste0("The distance was adopted for the minimal distance: ",dist.spot, "px"))}


  VF <- furrr::future_map(.x=1:nrow(df), .f=function(i){

    #Spot Def.
    bc <- df[i, c("barcodes")]
    cc <- df[i, c("x", "y")]
    #Neighbour Spots
    NN <-
      NN.file %>%
      dplyr::filter(xo < cc$x+dist.spot & xo > cc$x-dist.spot) %>%
      dplyr::filter(yo < cc$y+dist.spot & yo > cc$y-dist.spot) %>%
      dplyr::pull(bc_destination)

    #Filter input DF
    NN.df <- df %>% dplyr::filter(barcodes %in% NN) %>% as.data.frame()
    parameter <- features

    #Create Vector
    V <- -c(as.numeric(cc) - c(NN.df$x[which.max(NN.df[,parameter])], NN.df$y[which.max(NN.df[,parameter])]))

    if(length(V)==0){out <- data.frame(barcodes=bc, t.x=0, t.y=0)}else{out <- data.frame(barcodes=bc, t.x=V[1], t.y=V[2])}

    return(out)



  }, .progress = T) %>%
    do.call(rbind, .) %>%
    as.data.frame()


  out <- cbind(df, VF)
  out[is.na(out)] <- 0

  return(out)




}

• plotVectorFields.R

plotVectorFields <- function(VF, parameter, pt.size=6,pt.alpha=0.8,color.extern=NULL,skip=1){

  VF <-
    VF %>%
    dplyr::select(x,y,{{parameter}}, t.x, t.y) %>%
    dplyr::rename("parameter":=!!sym(parameter))


  color.points <- VF$parameter
  if(!is.null(color.extern)){color.points <- color.extern}

  if(color.points %>% class()=="factor"){
  p <-
    ggplot2::ggplot(data=VF, aes(x,y))+
    ggplot2::geom_point(data=VF , mapping=aes(x,y, color=color.points), size=pt.size, alpha=pt.alpha)+
    metR::geom_vector(aes(dx = t.x, dy = t.y),skip=skip) +
    metR::scale_mag()+
    ggplot2::theme_void()+
    Seurat::NoLegend()
  }else{
    p <-
      ggplot2::ggplot(data=VF, aes(x,y))+
      ggplot2::geom_point(data=VF , mapping=aes(x,y, color=color.points), size=pt.size, alpha=pt.alpha)+
      ggplot2::scale_color_viridis_c(guide = "none")+
      metR::geom_vector(aes(dx = t.x, dy = t.y),skip=skip) +
      ggplot2::theme_void()+
      Seurat::NoLegend()+
      metR::scale_mag()
        }

  return(p)
}

关于轨迹分析，推荐大家采用SPATA2。

生活很好，有你更好,祝打赏的老板们国庆快乐