空间轨迹向量场

作者,追风少年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。

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

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
禁止转载,如需转载请通过简信或评论联系作者。
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 214,444评论 6 496
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 91,421评论 3 389
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 160,036评论 0 349
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 57,363评论 1 288
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 66,460评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 50,502评论 1 292
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 39,511评论 3 412
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 38,280评论 0 270
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,736评论 1 307
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,014评论 2 328
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 39,190评论 1 342
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,848评论 5 338
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 40,531评论 3 322
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 31,159评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 32,411评论 1 268
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 47,067评论 2 365
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 44,078评论 2 352

推荐阅读更多精彩内容