SIMBA方法解讀

fig1

預處理

scRNA-seq

過濾掉在少于三個細胞中表達的基因。原始計數按文庫大小標準化，然后進行對數轉換。可選地，可以執行HVG選擇以刪除非信息性基因并加速訓練過程。在將特征輸入限制為通過HVG選擇確定的特征輸入時，未觀察到所得細胞嵌入的顯著差異，但不會生成非可變基因的 SIMBA 嵌入，因為它們未在圖中編碼。

scATAC-seq

過濾掉少于三個細胞中的peak。或者，實施一個可擴展的基于截斷 SVD 的程序來選擇峰，作為初步步驟，以額外過濾非信息峰并加速訓練過程。首先，選擇前 k 個主成分 (PC)，其中 k 的選擇基于方差圖。然后，對于前 k 個 PC 中的每一個，使用由“kneed”實現的拐點檢測算法根據載荷自動選擇peak。最后，將為每個 PC 選擇的峰組合起來并表示為“變量峰”。與使用 scRNA-seq 數據的觀察結果類似，變量峰選擇的可選步驟對生成的細胞嵌入的影響可以忽略不計。盡管對生成的嵌入的影響微乎其微，但此特征選擇步驟在減少訓練過程時間方面具有顯著的實際優勢。

使用 JASPAR2020 中的“Biostrings”和“motifmatchr”包執行 k-mer 和motif掃描。SIMBA 的實現中包含一個方便的 R 命令行腳本“scan_for_kmers_motifs.R”，它將peak列表（格式為 bed 文件）轉換為稀疏的peaks-by-k mers和peaks-by-motifs矩陣，該矩陣存儲為 hdf5 格式的文件。

圖構建（5種場景）

scRNA-seq分析

在構建細胞和基因graph時，如果基因在給定細胞中表達，則在細胞和基因之間添加一條邊。為了區分每條邊的強度，提出了一種bins方法，將基因表達值分為不同的級別，同時保留原始分布。不同級別的基因表達由不同類型的關系編碼。具體來說，首先使用基于 k-means 的程序近似歸一化基因表達矩陣中非零值的分布。首先，將連續的非零值分箱到 n 個區間（默認情況下，n = 5）。使用一維 k 均值聚類定義箱寬，其中每個箱中的值分配給相同的聚類中心。然后將連續矩陣轉換為離散矩陣，其中 1、…、n 用于表示 n 個基因表達級別。零值保留在此矩陣中。然后，通過將兩種類型的實體（細胞和基因）編碼為節點，將具有 n 個不同權重的關系（即 n 個基因表達級別）編碼為邊來構建圖。這 n 個關系權重的范圍從 1.0 到 5.0，步長為 5 / n，表示基因表達水平（最低：1.0，最高：5.0），因此與高表達水平相對應的邊對嵌入的影響比中等或低表達水平的邊更大。正如預期的那樣，觀察到，隨著bins數量的增加，離散化分布接近原始分布。然而，表達分辨率的增加對生成的嵌入影響不大。此離散化是在 SIMBA 包中使用函數“si.tl.discretize()”實現的。

除了關系類型權重外，SIMBA 還支持在構建圖時將基因表達值直接編碼為邊權重。此過程會生成與分箱過程類似的嵌入。這進一步表明離散化bins在捕獲生物信息方面是有效的。這種對邊權重的支持是在 SIMBA 包中使用函數“si.tl.gen_graph(add_edge_weights=True)”實現的。

scATAC-seq分析

peak-by-cell矩陣被二值化：“1”表示峰內至少有一個read，否則分配“0”。該圖是通過將兩種類型的實體（細胞和峰）編碼為節點，將它們之間的關系（表示給定峰在細胞中的存在）編碼為邊來構建的。單個關系類型的權重為 1.0。當 DNA 序列特征可用時，它們被使用 k-mer 和motif實體作為節點編碼到圖中。這是通過首先將peak-by-k mer或peak-by-motif矩陣二值化，然后使用peak、k-mer 和motif作為節點，并使用peak內這些實體的存在作為這些額外節點和峰節點之間的邊來構建原始peak-by-cell圖的擴展。k-mer 和峰之間的關系被分配了 0.02 的權重，而 TF 基序之間的關系被分配了 0.2 的權重。值得注意的是，根據具體的分析任務，k-mers 和motif可以彼此獨立地用作圖的節點輸入。

多模態分析

將上述使用 scRNA-seq 和 scATAC-seq 數據構建圖的策略結合起來，構建了多組學圖。

批次整合

按照“scRNA-seq分析”中所述構建每個批次的圖。通過基于截斷隨機 SVD 的程序推斷不同批次細胞之間的邊緣，以鏈接不同批次的不相交圖。更具體地說，在 scRNA 序列數據的情況下，考慮兩個基因表達矩陣 $X1_{n_{1}\times m}$ 和 $X2_{n_{2}\times m}$ ，其中 $n_{1}$ 和 $n_{2}$ 分別是兩個批次的細胞數量， $m$ 是gene數量。

然后計算： $X=X1\times X2^{T}$ 隨后對 $X$ 執行截斷隨機 SVD： $X=U\times \Sigma\times V^{T}$ 其中， $U$ 是 $n_{1}\times d$ 的矩陣， $\Sigma$ 是 $d\times d$ 的矩陣， $V$ 是 $n_{2}\times d$ 的矩陣，默認 $d = 20$ 。

$U$ 和 $V$ 都進一步進行了 L2 歸一化。對于 $U$ 中的每個細胞，我們在 $V$ 中搜索 k 個最近鄰居，反之亦然（默認情況下，k = 20）。最終，只有 $U$ 和 $V$ 之間的相互最近鄰居被保留為細胞之間的邊（注意是推斷的邊）。推斷不同批次細胞之間的邊的過程在 SIMBA 包中的函數“si.tl.infer_edges()”中實現。

對于多個批次，SIMBA 可以靈活地推斷任意一對batch-pair之間的邊。然而，在實踐中，邊是在最大的數據集或包含最完整預期細胞類型集的數據集與其他數據集之間推斷的。

多模態整合

scRNA-seq 和 scATAC-seq 圖分別按照“scRNA-seq 分析”和“scATAC-seq 分析”中的步驟構建。為了推斷 scRNA-seq 和 scATAC-seq 細胞之間的邊，首先計算 scATAC-seq 數據的基因活性分數（gene activity score）。更具體地說，對于每個基因，考慮 TSS （轉錄起始位點）上游和下游 100 kb 內的peak。與基因體區域重疊或在基因體上游 5 kb 內的peak的權重為 1.0。否則，使用指數衰減函數根據peak value與 TSS 的距離對其進行加權： $exp(\frac{-distance}{5000})$ 。隨后，將每個基因的gene score計算為所考慮峰值的加權和。然后將這些基因得分縮放到相應的基因大小。這些步驟由 SIMBA 中的函數“si.tl.gene_scores()”實現。為了方便用戶，SIMBA 包整理了幾個常用參考基因??組的基因注釋，包括 hg19、hg38、mm9 和 mm10。一旦獲得基因得分，就執行“批次整合”中描述的相同程序，使用 SIMBA 中的函數“si.tl.infer_edges()”推斷 scRNA-seq 和 scATAC-seq 分析的細胞之間的邊。

生成圖的過程在 SIMBA 包中的函數“si.tl.gen_graph()”中實現。

圖學習

在構建生物實體之間的多關系圖之后，作者采用了知識圖譜和推薦系統中的圖嵌入技術來為這些實體構建無監督表示。

提供一個input無向圖 $G = (V, E)$ ，其中 $V$ 是一組實體（節點） $E$ 是一組邊，在源實體 $u$ 和目標實體 $v$ 之間存在通用邊 $e = (u, v)$ 。進一步假設每個實體都有不同的已知類型（例如，細胞或peak）。

圖嵌入方法通過隨機梯度下降優化edge預測目標，為每個 $v \in V$ 學習一個 $D$ 維嵌入向量，其中實驗中使用 $D = 50$ 。實體 $v$ 的embedding記為 $\theta_{v}$ 。

對于edge $e = (u, v)$ ，記 $s_{e}=\theta_{u}\cdot\theta_{v}$ 為 $e$ 的得分，損失為： $L_{e}=-log\frac{exp(s_{e})}{\sum_{e'\in N}exp(s_{e'})}w_{e}$ 其中， $N$ 是通過破壞 $e$ 生成的一組“負樣本”候選邊， $w_e$ 是邊權重，默認情況下是關系權重，但在每種關系類型中可能因邊而異。例如，細胞和基因之間的邊可以編碼為具有不同邊權重的單一關系，這些邊權重編碼標準化的基因表達水平（見“scRNA-seq分析”）。

通過將目標邊 $e = (u, v)$ 中的源實體或目標實體替換為隨機采樣的實體來構建負樣本。因此，例如，對于cell-peak的邊，僅對cell和peak實體之間的負候選樣本進行采樣。這種設置至關重要，因為大多數隨機選擇的邊是無效的（例如，峰-峰）。

使用了 PyTorch-BigGraph 框架，該框架可以高效計算多種實體類型的多關系圖嵌入，并且可以擴展到包含數百萬或數十億個實體的圖。對于 130 萬個細胞，PyTorch-BigGraph 訓練本身僅需大約 1.5 小時，使用 12 個 CPU 核心，無需 GPU。

SIMBA空間中查詢實體

信息豐富的 SIMBA 嵌入空間可用作實體（包括細胞和特征）的數據庫。為了在“SIMBA 數據庫”中查詢給定細胞或特征的鄰近實體，我們首先根據其 SIMBA 嵌入構建所有實體的 k-d 樹。然后，使用歐幾里得距離在樹中搜索最近的鄰居。為此，SIMBA 查詢可以在指定半徑內執行 k 最近鄰居 (KNN) 或最近鄰居搜索。SIMBA 還提供了將搜索限制為某些類型實體的選項，當某種類型的實體數量遠遠超過其他實體時，這很有用。例如，給定細胞的 k 個最近特征可能都是峰值，而基因是感興趣的特征。在這種情況下，SIMBA 允許用戶添加“過濾器”以確保在指定類型的實體內執行最近鄰居搜索。此過程在函數“st.tl.query()”中實現，其可視化在 SIMBA 包中的函數“st.pl.query()”中實現。

識別TF-target genes

為了推斷給定主調節因子的靶基因，我們假定，在共享的 SIMBA 嵌入空間中，（1）靶基因靠近 TF 基序和 TF 基因，表明靶基因的表達與 TF 的表達和 TF 基序的可及性高度相關，并且以細胞類型特異性的方式呈現；（2）靶基因位點附近的可及區域（峰）必須靠近 TF 基序和靶 TF 基因，表明靶基因位點附近的順式調控元件的可及性與 TF 的表達和 TF 基序的可及性高度相關，并且以細胞類型特異性的方式呈現。

給定一個主調節因子，通過比較 SIMBA 共嵌入空間中 TF 基因、TF 基序和候選靶基因基因組位點附近的峰的位置來識別其靶基因。

更具體地說，我們首先分別搜索該主調節因子的基序（TF 基序）和基因（TF 基因）周圍的 k 個最近鄰基因（默認 k = 200）。這些鄰居基因的并集就是初始的候選靶基因集。然后根據以下標準對這些基因進行篩選：假定靶基因 TSS 上游和下游 100 kb 內的開放區域（峰）必須包含 TF 基序。

接下來，對于每個候選靶基因，我們計算了 SIMBA 嵌入空間中的四種距離：(1) 候選靶基因與 TF 基因的嵌入之間的距離；(2) 候選靶基因與 TF 基序的嵌入之間的距離；(3) 候選靶基因與 TF 基序的基因組位點附近的峰之間的距離；以及 (4) 候選靶基因與其基因組位點附近的峰之間的距離。所有距離（默認為歐幾里得距離）都轉換為所有基因或所有峰之間的等級，以使距離在不同的主調節器之間具有可比性。

最終的靶基因列表由計算出的排名決定，使用兩個標準：（1）TF 基因或 TF 基序最近的峰值中至少有一個在預定范圍內；（2）候選靶基因的平均排名在預定范圍內。此過程在 SIMBA 中的函數“st.tl. find_target_genes ()”中實現。