Zenn
🐙

【Visium】隣接するSpotの同種接続/異種接続をスコアリング

2024/08/08に公開
R
graph
bioinformatics
singlecell
visium
tech

次の論文では、Visium spotのアノテーションの連続性についてスコアリングを行っている。同じ細胞種spotが連続しているかどうかをHomotypic score、異なる細胞種同士での接続でどの細胞種とどの細胞種が共起/排他関係なのかをHeterotypic scoreで表現している。
https://www.sciencedirect.com/science/article/pii/S1550413121003636#undfig1

Homotypic score

Visium spotがどれだけ同じ成分（クラスターなど）で繋がっているかをスコアリング。

  1. 同じ細胞アノテーション同士の接続（同種接続）を作ってその次数の平均を計算
    下図Cのdegreeが1つのspotが何個接続を持つか（次数）を示している。
  2. ランダムにアノテーションをシャッフルした場合での同種接続の次数平均を計算
    これを50回繰り返して50回の平均を代表値として使用する。
  3. 1と2との差をスコアとする。

    原著 Figure.2

Heterotypic score

隣接するVisium spotの異なる細胞間での接続をスコアリング。

  1. 異種細胞間での接続（異種接続）数を集計して隣接行列を作成
    隣接行列は細胞種数*細胞種数の行列で、どの細胞とどの細胞が何回隣接しているかをまとめたものである。
  2. ランダムにアノテーションをシャッフルした場合での異種接続数を集計
    これを50回繰り返してその平均と標準偏差を計算する。
  3. 1の隣接行列を、2で求めて平均とSDでZ score化する。

    原著 Figure.3
    図中ではM2マクロファージと脂肪細胞の前駆細胞の共起性をピックアップしている。



GitHubには計算コードを含めかなり詳細に情報を公開してくれている。論文中のVisiumデータもdepositされている。有難い。

GitHub: https://github.com/lfranzen/scwat-st
deposit data: https://data.mendeley.com/datasets/3bs5f8mvbs/1

Homotypic score、Heterotypic score計算機能が定義: https://github.com/lfranzen/scwat-st/blob/master/scripts/nbs_analysis_functions.R
それら機能を使った処理スクリプトの例: https://github.com/lfranzen/scwat-st/blob/master/scripts/other_nbs_analysis_10x_data.R


準備

STutilityパッケージ

解析にはSeuratに加えて、Visiumデータ特有の画像や座標の処理にSTutilityパッケージの機能を多用する。現在は更新が終わって後継のsemlaパッケージが公開されているが、STutilityもまだ使用可能であった。

https://ludvigla.github.io/STUtility_web_site/Installation.html

devtools::install_github("jbergenstrahle/STUtility")


deposit data

Mendeley dataページのDownload Allからダウンロードし、zipを解凍しておく。10検体分のVisiumデータがある。

https://data.mendeley.com/datasets/3bs5f8mvbs/1

本記事のデモでは、process済みのSeuratオブジェクト「se-object.visium_baseline.rds」を使用するが、Visium画像や座標の読み込みには個々のSpaceRanger出力が必要となるため全てのデータを落としておく。（GitHubでは、git cloneして「data/visium」フォルダにdeposit dataを配置するように書いてある。）
SpaceRanger出力から開始したい場合は、次のRmdファイルで原著の処理をなぞることができる。

https://github.com/lfranzen/scwat-st/blob/master/scripts/visium_baseline.main_analysis.Rmd


データ読み込み

library(Seurat)
library(dplyr)
library(ggplot2)
library(STutility)

処理済みのSeuratオブジェクトを読み込む。※ パスは自分の環境に合わせる。

sc <- readRDS("data/visium/se-object.visium_baseline.rds/se-object.visium_baseline.rds")

既にサンプル間のバッチ補正やspotのアノテーション情報が登録されている。

DimPlot(sc, group.by = "cluster_anno", reduction = "umapICA") + 
    DimPlot(sc, group.by = "sample_id", reduction = "umapICA") & 
    theme(aspect.ratio = 1)


画像の登録

ここではGitHubコードに倣ってSTutility::LoadImages()機能で画像を登録する。読み込んだSeuratオブジェクトに既に画像へのファイルパスが登録されているのだが、著者らの環境でのファイルパスになっているので上記でDownloadした各症例のSpaceRangerフォルダが配置されたパスで上書きしておく。

画像のファイルパス情報を上書き
sc@tools[["Staffli"]]@imgs <- 
  gsub(pattern = "/Users/lovisa.franzen/Documents/ST_Adipose/scwat-st",
     replacement = "C:/Users/...", # <-- 自分の環境に合わせる
     sc@tools[["Staffli"]]@imgs)

パスが間違えでなければ次の1行で画像情報などがSeuratオブジェクトに登録される。

sc <- STutility::LoadImages(sc)

sc@tools$Staffli@meta.dataに各spotの座標情報などが追加される。

sample列が連番になっているのでファイル名中のサンプル名を追記しておく。

meta <- sc@tools$Staffli@meta.data
# サンプル名をキャプチャー Sと数字2文字の箇所
samplenames <- sub(pattern = ".*(S\\d{2}).*",
                   replacement = "\\1",
                   sc@tools$Staffli@imgs)

meta$sample_name <- samplenames[as.numeric(meta$sample)]
sc@tools$Staffli@meta.data <- meta



以下で隣接関係をスコア化するコードを紹介していく。

1昔前のVisiumは大きめのspotサイズが等間隔に並んでおり、接続の最大数は6つである。また接続距離も一定であるため距離の閾値設定は容易。組織の辺縁でなければ多くのspotは6つの隣接spotを持つ。
論文ではDBSCANパッケージのkNN()機能で6近傍を見つけて距離閾値を満たす接続を返している。STutility::GetSpatNet()機能でできる。

一方で、昨今のVisium HDやXeniumでは1細胞レベルのセグメンテーションデータを扱うことができる。組織中の細胞は整列していないため、接続距離や接続最大数は一定ではない。どこまで遠くの位置まで隣接と捉えるかの距離閾値は検討しながら自分で決める必要がある。

本記事では昨今の潮流に対応するためドロネー三角形分割を用いた接続探索でスクリプトを更新する。この方法であればVisium spotでも細胞セグメンテーションデータでもどちらにも対応できる。

【 ドロネー三角形分割 】

deldirパッケージを使用する。

library(deldir) # install.packages("deldir")

X座標、Y座用のdata.frameを用意するだけで進む。デモとして1症例分の処理例を示す。

xy <- sc@tools$Staffli@meta.data[,c("pixel_x","pixel_y", "sample_name")]
xy <- split(xy, f=xy$sample_name) # 症例で分ける

sample_tmp <- xy[[1]] # デモ用に1症例取り出す

# 1症例のXY座標情報を渡す。
res <- deldir(x = sample_tmp$pixel_x,
              y = sample_tmp$pixel_y)

出力のリストの内、delsgsの箇所がspot間の接続を示している。ind1、ind2が接続元、接続先の行番号である。この情報でspot barcodeを紐づけられる。


接続の確認

ドロネー三角形分割では全ての座標がどこかしらに繋がっており、非常に離れた箇所でも繋がり得る。deldir()の出力をplot()コマンドに渡すと接続が簡単に確認できる。

組織の端と端が繋がってしまっていることがわかる。

plot(res)


ドロネー図（実線部分。点線はボロノイ図）


距離閾値の設定

遠い接続をフィルタリングするための閾値を設定する。

そのためにまず2 spot間のユークリッド距離を計算する。

edges <- res$delsgs

edges$dist <- sqrt((edges$x1-edges$x2)^2 + (edges$y1-edges$y2)^2)


ユークリッド距離を追加



距離閾値を検証する例を紹介する。
座標情報はピクセル座標であるため、何ピクセル離れれば何μm離れているかは解像度情報が無ければ算出できない。適当に値を入れながら閾値を決定すると良い。特に細胞レベルのセグメンテーションの場合は適切な距離閾値を目で確認するのは大事なステップである。

ドロネー三角形分割結果のdelsgsを上書きしてplot()に渡す。

threshold <- 30

tmp <- res
tmp$delsgs <- edges[edges$dist < threshold, ]
plot(tmp)

実線の繋がりが変に遠くに繋がらないようにフィルタリングする。（点線はボロノイ図の線なので気にしなくてよい。）


接続距離の閾値でフィルタリングしたドロネー図

threshold <- 30
edges <- edges[edges$dist < threshold, ]


フィルタリング前後の接続数


無向グラフのエッジ関係に修正

ドロネー三角形分割結果は2つのspot間の接続につき1つのind1とind2がある。有向グラフであれば接続元（ind1）と接続先（ind2）のままでも良いが、ここでは無向グラフであり、接続元と接続先に順序は無い。
そこでind1とind2を逆にした情報も追加する。その結果、接続数は2倍になる。

edges_inv <- edges
colnames(edges_inv) <- c("x2", "y2", "x1", "y1", "ind2", "ind1", "dist")
edges <- rbind(edges,edges_inv[,colnames(edges)])


データ成型

最後に、index番号をspotのバーコードに修正しておく。

edges$from <- rownames(sample_tmp)[edges$ind1]
edges$to <- rownames(sample_tmp)[edges$ind2]


最後に関数化して、複数症例で実施しておく。

GetNetwork <- function(
  x,y,threshold,barcodes
){
  res <- deldir(x,y)
  edges <- res$delsgs
  edges$dist <- sqrt((edges$x1-edges$x2)^2 + (edges$y1-edges$y2)^2)
  edges <- edges[edges$dist < threshold, ]

  edges_inv <- edges
  colnames(edges_inv) <- c("x2", "y2", "x1", "y1", "ind2", "ind1", "dist")
  edges <- rbind(edges,edges_inv[,colnames(edges)])
  
  edges$from <- barcodes[edges$ind1]
  edges$to <- barcodes[edges$ind2]
  
  return(edges)
}

xy <- sc@tools$Staffli@meta.data[,c("pixel_x","pixel_y", "sample_name")]
edge_list <- split(xy, f=xy$sample_name)
edge_list <- lapply(edge_list, function(xy){
  GetNetwork(x = xy$pixel_x,
             y = xy$pixel_y,
             threshold = 30,
             barcodes = rownames(xy))
})


10症例分のedge情報リスト

【 Homotypic Score 】

1つのVisium spotのアノテーションが、隣接するspotにも見られるとき、その接続を同種接続と呼ぶこととする。一般的に1つのノード（ここではvisium spot）が持つ接続（エッジ）数を次数（degree）と呼ぶ。なお本記事ではspotを細胞、clusterをアノテーションや細胞種と記載することもある。

平均次数とは言葉通り1つのノードが持つ次数の平均である。Homotyipic scoreの計算では全次数ではなく、同種接続の次数の平均を計算する。サンプルごとに細胞種の存在量比は異なるため、その補正が必要である。そこで1つのサンプル内で細胞種の存在量比は一定でその座標だけランダムにシャッフルした場合に得られる平均同種次数を計算して、補正に用いる。（ある細胞種が沢山いれば、平均同種次数も高く出る。その場合、ランダムシャッフルしても平均同種次数は高く得られる。一方、意味を持って集合を作っていればランダムシャッフルでの平均同種次数との差が大きくなる。）


同種細胞間の平均次数を計算

原著は初めから複数サンプル用になっているので、ここでは1サンプル内で行う例をまずは紹介する。各spotの細胞アノテーション情報を使用して、同じアノテーション同士の接続をまずは算出していく。

edge <- edge_list[[1]]

接続元、接続先のspot barcodeからアノテーション名を登録する。

column.cluster.id <- "cluster_anno"

edge$cluster_from <- sc@meta.data[edge$from, column.cluster.id]
edge$cluster_to <- sc@meta.data[edge$to, column.cluster.id]

各クラスターのspotを抜き出し、1つのspotが持つ同種接続数の平均を計算する。

cluster_level <- sort(unique(sc@meta.data[, column.cluster.id]))
k_avg <- numeric(length = length(cluster_level))
names(k_avg) <- cluster_level

for( cluster in cluster_level){
  # 接続元ノードがclusterのエッジを取得
  edge_subset <- subset(edge, cluster_from == cluster)
  # さらに接続先も同じclusterのものを取得
  edge_subset_homo <-subset(edge_subset, cluster_to == cluster)
  
  
  # 1つの接続元ノードに同種接続が何回あるか集計
  ki <- edge_subset_homo %>% 
    dplyr::group_by(from) %>% 
    dplyr::count(name = "ki") # 集計結果の列名はki
  
  # ユニークな接続元ノード数
  N <- length(unique(edge_subset$from))  
  
  # 平均次数を取得 同種接続数で割るのではなく、1つのクラスターの合計spot数
  # --> あるクラスターの1 spotに平均何個の同種接続があるか。  
  k_avg_cluster <- (sum(ki$ki))/N
  
  k_avg[cluster] <- k_avg_cluster
  }

クラスターごとに同種接続を集計して、そのクラスターのspot数で割って平均を取得している。これが1つの症例内での同種接続の平均次数である。

par(mar = c(12,3,1,1))
barplot(k_avg, las=2)



これを関数化して全症例で回しておく。

CalculateAvgDegree <- function (
    seuratObject,
    edge,
    column.cluster.id
){
  cluster_level <- sort(unique(seuratObject@meta.data[, column.cluster.id]))
  
  edge$cluster_from <- seuratObject@meta.data[edge$from, column.cluster.id]
  edge$cluster_to <- seuratObject@meta.data[edge$to, column.cluster.id]
  
  cluster_level <- sort(unique(seuratObject@meta.data[, column.cluster.id]))
  k_avg <- numeric(length = length(cluster_level))
  names(k_avg) <- cluster_level
  
  for( cluster in cluster_level){
    # 接続元ノードがclusterのエッジを取得
    edge_subset <- subset(edge, cluster_from == cluster)
    # さらに接続先も同じclusterのものを取得
    edge_subset_homo <-subset(edge_subset, cluster_to == cluster)
    
    
    # 1つの接続元ノードに同種接続が何回あるか集計
    ki <- edge_subset_homo %>% 
      dplyr::group_by(from) %>% 
      dplyr::count(name = "ki") # 集計結果の列名はki
    
    # ユニークな接続元ノード数
    N <- length(unique(edge_subset$from))  
    
    # 平均次数を取得 1つのクラスターの合計spot数で割る
    # --> あるクラスターの1 spotに平均何個の同種接続があるか。  
    k_avg_cluster <- (sum(ki$ki))/N
    
    k_avg[cluster] <- k_avg_cluster
  }
  
  return(k_avg)
}
全症例に繰り返し処理
k_ave_list <- lapply(edge_list, function(x){
  CalculateAvgDegree(seuratObject = sc, edge = x, column.cluster.id = "cluster_anno")
})

参考までに全症例の同種接続の平均次数を1つにまとめてヒートマップで視覚化してみる。

k_ave_df <- do.call(rbind, k_ave_list)

pheatmap::pheatmap(k_ave_df)

症例ごとにspot数は違うので症例間を直接比較するのは尚早だが、この時点でもPreadipocyte3やVascular2の同種接続が症例によって偏っているのが分かる。


ランダムサンプルデータで平均同種次数を計算

どのspotとどのspotが接続しているかの情報は変わらない。spotのアノテーション情報をランダムシャッフルして、平均同種次数を計算する。

1症例で処理例を示す。

edge <- edge_list[[1]]

1症例が持つアノテーション列を抽出し、ランダムにシャッフルする。この時、アノテーションの存在量比は変わらないので、サンプル固有のアノテーションの偏りは保持される。

column.cluster.id <- "cluster_anno"

# アノテーションベクトルを取得
annotations <- sc@meta.data[unique(edge$from), column.cluster.id]
# シャッフル
annotations <- sample(annotations)
# アノテーションベクトルの名前属性にシャッフル前のバーコードを登録
names(annotations) <- unique(edge$from)

# 接続元、接続先のspotのアノテーション情報にシャッフルしたアノテーションを登録
edge$cluster_from <- annotations[edge$from]
edge$cluster_to <- annotations[edge$to]

後は先述した平均同種次数の計算と同じである。

cluster_level <- sort(unique(sc@meta.data[, column.cluster.id]))
k_avg <- numeric(length = length(cluster_level))
names(k_avg) <- cluster_level

for( cluster in cluster_level){
  # 接続元ノードがclusterのエッジを取得
  edge_subset <- subset(edge, cluster_from == cluster)
  # さらに接続先も同じclusterのものを取得
  edge_subset_homo <-subset(edge_subset, cluster_to == cluster)
  
  
  # 1つの接続元ノードに同種接続が何回あるか集計
  ki <- edge_subset_homo %>% 
    dplyr::group_by(from) %>% 
    dplyr::count(name = "ki") # 集計結果の列名はki
  
  # ユニークな接続元ノード数
  N <- length(unique(edge_subset$from))  
  
  # 平均次数を取得 同種接続数で割るのではなく、1つのクラスターの合計spot数
  # --> あるクラスターの1 spotに平均何個の同種接続があるか。  
  k_avg_cluster <- (sum(ki$ki))/N
  
  k_avg[cluster] <- k_avg_cluster
  }


原著ではこれを50回繰り返して、その平均値を補正用データとして使用している。繰り返しも含めて関数化しておく。返り値には平均と標準偏差も含めたリストにしている。

ランダムシャッフルしたアノテーションでの平均同種次数を取得する機能を定義
RandomClustersAvgDegree <- function (
    seuratObject,
    edge,
    column.cluster.id,
    nperm=50
){
  cluster_level <- sort(unique(seuratObject@meta.data[, column.cluster.id]))
  
  k_avg_list <- list()
  
  for( i in 1:nperm){
    annotations <- seuratObject@meta.data[unique(edge$from), column.cluster.id]
    annotations <- sample(annotations)
    names(annotations) <- unique(edge$from)
    
    edge$cluster_from <- annotations[edge$from]
    edge$cluster_to <- annotations[edge$to]
    
    k_avg <- numeric(length = length(cluster_level))
    names(k_avg) <- cluster_level
    
    for( cluster in cluster_level){
      # 接続元ノードがclusterのエッジを取得
      edge_subset <- subset(edge, cluster_from == cluster)
      # さらに接続先も同じclusterのものを取得
      edge_subset_homo <-subset(edge_subset, cluster_to == cluster)
      
      
      # 1つの接続元ノードに同種接続が何回あるか集計
      ki <- edge_subset_homo %>% 
        dplyr::group_by(from) %>% 
        dplyr::count(name = "ki") # 集計結果の列名はki
      
      # ユニークな接続元ノード数
      N <- length(unique(edge_subset$from))  
      
      # 平均次数を取得 1つのクラスターの合計spot数で割る
      # --> あるクラスターの1 spotに平均何個の同種接続があるか。  
      k_avg_cluster <- (sum(ki$ki))/N
      
      k_avg[cluster] <- k_avg_cluster
    }
    k_ave_list[[i]] <- k_avg
  }
  
  k_ave_df <- do.call(rbind,k_ave_list)
  k_ave_mean <- colMeans(k_ave_df)
  k_ave_sd <- matrixStats::colSds(k_ave_df)
  k_ave_sd[ k_ave_sd == 0] <- mean(k_ave_sd) # もしSDが0の箇所があったら補正

  return(list(k_ave_df=k_ave_df, mean=k_ave_mean, sd=k_ave_sd))
}
1サンプルでの実行例
random_res <- RandomClustersAvgDegree(
  seuratObject = sc,
  edge = edge,
  column.cluster.id = "cluster_anno",
  nperm = 50
)

Adipocyte1やUnspecific1、Unspecific2はランダムシャッフルで同種接続が高く得られやすいことがわかる。



これを全症例で実行する。

randam_homo_list <- lapply(edge_list, function(x){
  RandomClustersAvgDegree(seuratObject = sc,
                          edge = x,
                          column.cluster.id = "cluster_anno",
                          nperm = 50)
})


10症例分のランダムシャッフルでの平均同種次数のリスト


Homotypic scoreの計算

原著の図では平均同種次数とランダムシャッフルした平均同種次数の差でHomotyipic scoreを求めている。GitHubのコードではZ scoreで求めるオプションも用意されていた。

まずは1症例の結果を示す。

k_ave <- k_ave_list[[1]]
random_res <- randam_homo_list[[1]]

差でスコア化する場合は、ランダムシャッフルの平均値を引くだけである。

homotyipic_score <- k_ave - random_res$mean

par(mar = c(12,3,1,1))
barplot(homotyipic_score, las=2)

このサンプルではB細胞やAdipocyte3の固まった局在がランダムシャッフルでは起こりえないような頻度で起きているようだ。



Zスコア版も見ておく。

homotyipic_score_z <- (k_ave - random_res$mean)/random_res$sd

par(mar = c(12,3,1,1))
barplot(homotyipic_score_z, las=2)

差のスコアと同様の傾向だが、B細胞の同種接続がより強調される結果となった。



最後に全症例でスコアを計算しておく。

homotypic_score_list <- list()
for( i in 1:length(k_ave_list)){
  k_ave <- k_ave_list[[i]]
  random_res <- randam_homo_list[[i]]
  
  homotypic_score <- k_ave - random_res$mean
  homotypic_score_list[[i]] <- homotypic_score
}
names(homotypic_score_list) <- names(k_ave_list)

homotypic_Zscore_list <- list()
for( i in 1:length(k_ave_list)){
  k_ave <- k_ave_list[[i]]
  random_res <- randam_homo_list[[i]]
  
  homotypic_score_z <- (k_ave - random_res$mean)/random_res$sd
  homotypic_Zscore_list[[i]] <- homotypic_score_z
}
names(homotypic_Zscore_list) <- names(k_ave_list)



homotypic_score_df <- do.call(rbind, homotypic_score_list)
pheatmap::pheatmap(homotypic_score_df)


homotypic_Zscore_df <- do.call(rbind, homotypic_Zscore_list)
pheatmap::pheatmap(homotypic_Zscore_df)

Preadipocyte3やVascular2の同種接続が症例によって偏っているのが改めて確認できた。


【 Heterotypic score 】

1つのspotが繋がるspotの細胞種を集計することを繰り返して細胞種数*細胞種数の行列にまとめる。各セルは2種類の細胞種間の合計接続数になる。異種接続に着目しているだけで処理上は同種接続も集計できる。
原著ではSTutility::RegionNeighbours()で各spotの隣接spotのアノテーションを集計しているが、ここではカスタムして、上述した同種接続での計算例を全組み合わせに拡張している。原著通りの集計値にはならないが、補正側も同じ計算で行うので、最終的なスコアは凡そ同様になる。


全ての接続組み合わせを集計

1症例での例を示す。あるclusterのspotを取得し、その接続先をtable()で集計することをclusterの数だけ繰り返している。

edge <- edge_list[[1]]

column.cluster.id <- "cluster_anno"
cluster_level <- sort(unique(sc@meta.data[, column.cluster.id]))

edge$cluster_from <- factor(sc@meta.data[edge$from, column.cluster.id],
                            levels = cluster_level)
edge$cluster_to <- factor(sc@meta.data[edge$to, column.cluster.id],
                          levels = cluster_level)

neighbor_clusters_list <- list()
for( cluster in cluster_level){
  # 接続元ノードがclusterのエッジを取得
  edge_subset <- subset(edge, cluster_from == cluster)
  # 接続先のclusterを集計
  neighbor_clusters_list[[cluster]] <- table(edge_subset$cluster_to)
}

AdjMat <- do.call(rbind,neighbor_clusters_list)


各接続を集計した隣接行列



関数化して全症例で回しておく。

CreateAdjMatrix()機能の定義
CreateAdjMatrix <- function(
  seuratObject,
  edge,
  column.cluster.id
){
  cluster_level <- sort(unique(seuratObject@meta.data[, column.cluster.id]))

  edge$cluster_from <- factor(sc@meta.data[edge$from, column.cluster.id],
                              levels = cluster_level)
  edge$cluster_to <- factor(sc@meta.data[edge$to, column.cluster.id],
                            levels = cluster_level)
  
  neighbor_clusters_list <- list()
  for( cluster in cluster_level){
    # 接続元ノードがclusterのエッジを取得
    edge_subset <- subset(edge, cluster_from == cluster)
    # 接続先のclusterを集計
    neighbor_clusters_list[[cluster]] <- table(edge_subset$cluster_to)
  }
  
  AdjMat <- do.call(rbind,neighbor_clusters_list)
  
  return(AdjMat)
  }
全症例で繰り返し処理
AdjMat_list <- lapply(edge_list, function(x){
  CreateAdjMatrix(seuratObject = sc,
                  edge = x,
                  column.cluster.id = "cluster_anno")
})


ランダムシャッフルして全ての接続組み合わせを集計

アノテーション情報をシャッフルして、接続組み合わせごとにその接続数を集計する。上述してきたことの組み合わせになるので解説は割愛させていただく。同種接続は細胞種数の長さのベクトルだが、全組み合わせは細胞種数*細胞種数の行列で集計する点は異なる。

次のように関数を定義し、全症例で回しておく。

RandomClustersAdjMat <- function(
    seuratObject,
    edge,
    column.cluster.id,
    nperm=50
){
  cluster_level <- sort(unique(seuratObject@meta.data[, column.cluster.id]))
  
  random_res_list <- list()
  for( i in 1:nperm){
    annotations <- factor(
      seuratObject@meta.data[unique(edge$from), column.cluster.id],
      levels = cluster_level
    )
    annotations <- sample(annotations) # シャッフル
    names(annotations) <- unique(edge$from)
    
    edge$cluster_from <- annotations[edge$from]
    edge$cluster_to <- annotations[edge$to]
    
    neighbor_clusters_list <- list()
    for( cluster in cluster_level){
      # 接続元ノードがclusterのエッジを取得
      edge_subset <- subset(edge, cluster_from == cluster)
      # 接続先のclusterを集計
      neighbor_clusters_list[[cluster]] <- table(edge_subset$cluster_to)
    }
    
    random_res_list[[i]] <- do.call(rbind,neighbor_clusters_list)
  }
  
  ## 各接続の平均とSDを取得
  mean_mat <- matrix(nrow = length(cluster_level),ncol = length(cluster_level))
  rownames(mean_mat) <- cluster_level
  colnames(mean_mat) <- cluster_level
  sd_mat <- mean_mat
  
  for( i in 1:nrow(mean_mat)){
    for( j in 1:ncol(mean_mat)){
      # 各セルの値をベクトルとして取得
      values <- sapply(random_res_list, function(x) x[i,j])
      
      mean_mat[i,j] <- mean(values)
      sd_mat[i,j] <- sd(values)
    }
  }
  
  sd_mat[ sd_mat == 0] <- mean(as.vector(sd_mat)) # もしSDが0の箇所があったら補正

  return(list(random_res_list=random_res_list, mean=mean_mat, sd=sd_mat))
}

randam_adj_list <- lapply(edge_list, function(x){
  RandomClustersAdjMat(seuratObject = sc,
                       edge = x,
                       column.cluster.id = "cluster_anno",
                       nperm = 50)
})


10症例分のランダムシャッフルでの隣接行列リスト


Heterotypic scoreの計算

ランダムな隣接行列を使って、隣接行列をZ score化する。

# 結果の入れ先
heterotypic_Zscore_list <- list()

for( i in 1:length(AdjMat_list)){
  adjmat <-AdjMat_list[[i]]
  random_adjmat <- randam_adj_list[[i]]
  
  heterotypic_score_z <- (adjmat - random_adjmat$mean)/random_adjmat$sd
  heterotypic_Zscore_list[[i]] <- heterotypic_score_z
}

names(heterotypic_Zscore_list) <- names(AdjMat_list)


Z score化した隣接行列。10症例分



同種接続と異種接続では、どうしても同種接続の方が多く存在する。2つを同時に評価しようとすると値のスケールが異なるので見づらい。

そこで同種接続の箇所はここではNAにしておく。（0でも良い。）

heterotypic_Zscore_list <- lapply(heterotypic_Zscore_list,function(x){
  diag(x) <- NA # 対角成分（同種接続）はNA
  return(x)
})



最後にHeterotypic scoreを可視化する。論文と凡そ同じスケール感の値が得られている。

pheatmap::pheatmap(heterotypic_Zscore_list[[1]])


pheatmap::pheatmap(heterotypic_Zscore_list[[10]])

この症例ではマクロファージとPreadipocyte3の接続が特徴的。


全症例まとめてheterotypic scoreを計算

論文のデータは全症例で代表的なheterotypic scoreを計算しているように見える。参考までに全症例を1つの集合として計算する例も紹介しておく。



spot間の接続は症例内で計算済みである。このedge_listを1つのedge集合にまとめてCreateAdjMatrix()機能に渡す。

adjmat_all <- CreateAdjMatrix(seuratObject = sc,
                              edge = do.call(rbind,edge_list),
                              column.cluster.id = "cluster_anno")



ランダムデータでは、先にedge集合を作ってからランダムシャッフルすると、異なる症例のアノテーション存在量比がバイアスを生む可能性がある。症例内でアノテーションの存在量比が変わらないようにシャッフルしてから結合する。

別途関数として定義したものを用意した。

症例内でアノテーションシャッフルして1つのedge集合を返す機能
RandomClustersAdjMat_allsample <- function(
    seuratObject,
    edge_list,
    column.cluster.id,
    nperm=50
){
  cluster_level <- sort(unique(seuratObject@meta.data[, column.cluster.id]))

    random_res_list <- list()
  for( i in 1:nperm){
    # 症例内でアノテーションをシャッフルしたedge情報を作成
    shuffled_edge_list <- lapply(edge_list, function(edge){
      annotations <- factor(
        seuratObject@meta.data[unique(edge$from), column.cluster.id],
        levels = cluster_level
      )
      annotations <- sample(annotations) # シャッフル
      names(annotations) <- unique(edge$from)
      
      edge$cluster_from <- annotations[edge$from]
      edge$cluster_to <- annotations[edge$to]

      return(edge)
    })
    # 全症例のedge情報を連結
    edge_merge <- do.call(rbind, shuffled_edge_list)
    
    neighbor_clusters_list <- list()
    for( cluster in cluster_level){
      # 接続元ノードがclusterのエッジを取得
      edge_subset <- subset(edge_merge, cluster_from == cluster)
      # 接続先のclusterを集計
      neighbor_clusters_list[[cluster]] <- table(edge_subset$cluster_to)
    }
    
    random_res_list[[i]] <- do.call(rbind,neighbor_clusters_list)
  }
  
  ## 各接続の平均とSDを取得
  mean_mat <- matrix(nrow = length(cluster_level),ncol = length(cluster_level))
  rownames(mean_mat) <- cluster_level
  colnames(mean_mat) <- cluster_level
  sd_mat <- mean_mat
  
  for( i in 1:nrow(mean_mat)){
    for( j in 1:ncol(mean_mat)){
      # 各セルの値をベクトルとして取得
      values <- sapply(random_res_list, function(x) x[i,j])
      
      mean_mat[i,j] <- mean(values)
      sd_mat[i,j] <- sd(values)
    }
  }
  
  sd_mat[ sd_mat == 0] <- mean(as.vector(sd_mat)) # もしSDが0の箇所があったら補正

  return(list(random_res_list=random_res_list, mean=mean_mat, sd=sd_mat))
}
ランダムシャッフルした場合の隣接行列の平均とSDを取得
random_adjmat_all <- RandomClustersAdjMat_allsample(
  seuratObject = sc,
  edge = edge_list,
  column.cluster.id = "cluster_anno",
  nperm = 50)

最後にZスコアを計算する。

heterotypic_score_z <- (adjmat_all - random_adjmat_all$mean)/random_adjmat_all$sd
diag(heterotypic_score_z) <- NA

pheatmap::pheatmap(heterotypic_score_z)

血管2種類の接続が最も強いのも論文通りの結果である。

Discussion