作者：yuanzan
鏈接：dyno使用教程--1個(gè)R包實(shí)現(xiàn)59種單細(xì)胞軌跡推斷分析
)
來(lái)源：微信公眾號(hào)
著作權(quán)歸作者所有次乓，任何形式的轉(zhuǎn)載都請(qǐng)聯(lián)系作者吓歇。

在上一期的《單細(xì)胞軌跡分析知多少--擬時(shí)間分析比較》中我們介紹了45種單細(xì)胞軌跡推斷分析軟件方法在以下4個(gè)方面的比較：

1. 準(zhǔn)確性
2. 可擴(kuò)展性
3. 穩(wěn)定性
4. 可用性

得出了幾項(xiàng)重要結(jié)論：

• 軌跡推斷（TI，trajectory inference）軟件方法輸入和輸出接口的標(biāo)準(zhǔn)化對(duì)于TI方法的廣泛應(yīng)用非常關(guān)鍵
• 不同TI方法的組合存在一定的互補(bǔ)性票腰，能夠高概率的得出與實(shí)際生物學(xué)模型相符的軌跡推斷結(jié)果
• 在1個(gè)項(xiàng)目的數(shù)據(jù)集上嘗試多種TI分析算法城看，一是可以提供多種算法證據(jù)支持，二是如果多種算法的結(jié)果有沖突有利于我們深度思考可能的具體生物學(xué)原因

A comparison of single-cell trajectory inference methods(Saelens et al., 2019)的作者做了一個(gè)R包--dyno為終端用戶提供完整的TI分析流程杏慰，dyno特點(diǎn)如下：

1. 統(tǒng)一59種TI方法的輸入輸出接口
2. 提供交互式指南工具测柠，可幫助用戶選擇最合適的TI方法
3. 簡(jiǎn)化了軌跡的解釋和可視化，包括根據(jù)基因表達(dá)或者cluster著色
4. 還可以進(jìn)行下游分析缘滥，例如潛在marker gene 的鑒定

dyno安裝

dyno是一個(gè)R包轰胁，需要R/Rstudio運(yùn)行環(huán)境，目前的存放地址是在github- https://github.com/dynverse/dyno完域，所以我們用devtools進(jìn)行安裝软吐，打開(kāi)Rstudio瘩将，輸入以下命令：

#install.packages("devtools")
devtools::install_github("dynverse/dyno")

在這個(gè)過(guò)程中會(huì)安裝一系列dynverse R包(dynwrap, dynplot, dynmethods, …)吟税。

如果是在linux系統(tǒng)使用，需要安裝udunits和ImageMagick

• Debian / Ubuntu / Linux Mint:
  sudo apt-get install libudunits2-dev imagemagick
• Fedora / CentOS / RHEL:
  sudo dnf install udunits2-devel ImageMagick-c++-devel

運(yùn)行TI方法需要用到Docker或者Singularity (version ≥ 3.0)姿现，如果電腦系統(tǒng)是Windows或者M(jìn)acOS肠仪，推薦使用docker，如果是在集群备典，推薦使用 Singularity异旧。

如果是win10系統(tǒng)，可以運(yùn)行Docker CE提佣。Singularity的安裝請(qǐng)參考https://www.sylabs.io/docs/吮蛹，如果已經(jīng)安裝了anaconda荤崇，可以執(zhí)行conda create -n singularity -c conda-forge singularity=3.0.1安裝

安裝過(guò)程遇到的問(wèn)題

API rate limit exceeded

因?yàn)閐yno存放在了Github，在dyno的安裝過(guò)程會(huì)使用GitHub API潮针，默認(rèn)的API限制是60次請(qǐng)求术荤，所以會(huì)遇到API rate limit exceeded的問(wèn)題，解決方式如下：

• 在Rstudio命令行執(zhí)行 usethis::browse_github_pat()建議一個(gè)GitHub token
• 執(zhí)行usethis::edit_r_environ() 添加以下環(huán)境變量GITHUB_PAT = 'your_github_token
• 重啟R/Rstudio (重新載入 GITHUB_PAT)
• 再次重新安裝:devtools::install_github("dynverse/dyno")

其他安裝問(wèn)題

如果沒(méi)有devtools是不能安裝dyno的每篷，所以有些時(shí)候devtools的安裝也會(huì)出現(xiàn)一些問(wèn)題瓣戚，例如網(wǎng)絡(luò)不穩(wěn)定等，這時(shí)可以采用國(guó)內(nèi)的源進(jìn)行安裝：install.packages('devtools', repos='https://mirrors.tuna.tsinghua.edu.cn/CRAN/')

在linux集群焦读，如果沒(méi)有root權(quán)限子库，安裝使用dyno會(huì)非常困難，我準(zhǔn)備了解決方案會(huì)在文章最后一節(jié)介紹矗晃。

windows10/MACOS用戶快速體驗(yàn)方式

dyno的整個(gè)安裝過(guò)程會(huì)非常漫長(zhǎng)仑嗅，如果是windows10或者M(jìn)ACOS用戶可以在安裝docker之后用我構(gòu)建的R鏡像：seqyuan/seqyuan-r:v0.0.1（R version 3.6.1）進(jìn)行體驗(yàn)，這個(gè)R鏡像安裝了dyno张症、tidyverse及Seurat等包无畔，可以免去安裝過(guò)程快速體驗(yàn)dyno的強(qiáng)大，具體使用會(huì)在文章最后一節(jié)介紹吠冤。

dyno的使用

dyno的官網(wǎng)https://dynverse.org/有詳細(xì)的使用步驟介紹浑彰，在這里我們就不再重復(fù)官網(wǎng)測(cè)試數(shù)據(jù)的結(jié)果，重點(diǎn)介紹Seurat分析的10X單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)的結(jié)果怎么樣和dyno對(duì)接拯辙。

主要步驟

主要步驟如下

準(zhǔn)備數(shù)據(jù)

為數(shù)據(jù)集選擇最佳方法

運(yùn)行TI方法

軌跡可視化

從生物學(xué)角度解釋軌跡

• Rooting
• Milestone labelling

基因表達(dá)在軌跡上的展示

• A global overview of the most predictive genes
• Lineage/branch markers
• Genes important at bifurcation points

準(zhǔn)備數(shù)據(jù)

首先打開(kāi)Rstudio（如果是linux郭变，打開(kāi)R命令行），載入以下R包

library(dyno)
library(tidyverse)
library(Matrix)
library(Seurat)

dyno對(duì)數(shù)據(jù)的包裝主要通過(guò)dynwrap包來(lái)實(shí)現(xiàn)

Gene expression data

dynwrap要求raw counts和normalised(log2)表達(dá)數(shù)據(jù)涯保，低表達(dá)的細(xì)胞诉濒、雙細(xì)胞、壞細(xì)胞應(yīng)該被提前過(guò)濾夕春。使用Seurat處理后的數(shù)據(jù)一般都包含了這些步驟未荒。

wrap_expression要求raw counts和normalised表達(dá)為 sparse matrix (dgCMatrix）(列為genes/features，行為細(xì)胞)

#讀入seurat處理后的rds文件
sdata <- readRDS(file = "/Users/yuanzan/Desktop/tmp/seurat_project.rds")

#添加raw counts和normalised expression
#seurat的矩陣需要進(jìn)行行列轉(zhuǎn)換及志，以使行為細(xì)胞片排，列為基因
dataset <- wrap_expression(
  counts = t(sdata@assays$RNA@counts),
  expression = t(sdata@assays$RNA@data)
)

#添加先驗(yàn)信息，這里添加的是開(kāi)始轉(zhuǎn)換的“細(xì)胞id”速侈，后期可視化可以根據(jù)具體的軌跡推斷結(jié)果進(jìn)行調(diào)整
dataset <- add_prior_information(
  dataset,
  start_id = "TTTGTTGCAACTCATG.wt"
)

#添加數(shù)據(jù)的cluster信息率寡，這里我們直接用“seurat_clusters”即可
dataset <- add_grouping(
   dataset,
   sdata$seurat_clusters
)

為數(shù)據(jù)集選擇最佳方法

如果是MACOS，可以執(zhí)行以下命令以獲取推薦的最佳TI方法

guidelines <- guidelines_shiny(dataset)
methods_selected <- guidelines$methods_selected

執(zhí)行guidelines <- guidelines_shiny(dataset)命令之后會(huì)彈出一個(gè)shiny頁(yè)面倚搬，如下圖
[圖片上傳失敗...(image-17d1d0-1579139138858)]
因?yàn)槲覀兘o了dataset冶共，所以guidelines頁(yè)面會(huì)自動(dòng)填寫(xiě)細(xì)胞數(shù)和基因數(shù)參數(shù)，我們?cè)侔杨A(yù)期的運(yùn)行時(shí)間和運(yùn)行內(nèi)存選擇設(shè)置一下，可以擴(kuò)大選擇TI方法的選擇范圍捅僵。

通過(guò)選擇對(duì)于推斷軌跡的預(yù)期以及左側(cè)關(guān)于內(nèi)存家卖、運(yùn)行時(shí)間等各項(xiàng)參數(shù)，guidelines為我推薦了幾個(gè)TI算法的組合庙楚，點(diǎn)擊右上角的Close & use關(guān)閉這個(gè)shiny頁(yè)面篡九。

methods_selected <- guidelines$methods_selected這個(gè)命令是把guidelines推薦給我們的方法賦值給methods_selected，這是一個(gè)字符串向量醋奠。

1. 通過(guò)命令dim(dataset$counts)獲取cell number和feature/gene number
2. 打開(kāi)網(wǎng)址http://guidelines.dynverse.org
3. 在guidelines網(wǎng)頁(yè)填寫(xiě)細(xì)胞數(shù)、gene數(shù)以及其他設(shè)置
4. 打?qū)吹腡I方法就是推薦給我們的塞祈，如下圖

運(yùn)行TI方法

model_slingshot <- infer_trajectory(dataset, methods_selected[1])
model_paga_tree <- infer_trajectory(dataset, methods_selected[2])
model_paga <- infer_trajectory(dataset, methods_selected[3])
model_mst <- infer_trajectory(dataset, methods_selected[4])

作者把相應(yīng)TI方法的源碼封裝在了docker鏡像里金刁，調(diào)用相應(yīng)方法的后臺(tái)操作其實(shí)是：先拉取TI方法對(duì)應(yīng)的鏡像，然后再R進(jìn)程內(nèi)部啟動(dòng)docker容器執(zhí)行算法议薪，所以如果你使用的方法事先沒(méi)有經(jīng)過(guò)docker pull尤蛮，那么算法的執(zhí)行時(shí)間其實(shí)包含了拉取鏡像所需要的時(shí)間，建議把常用到的TI方法先pull到本地斯议，這樣就能使腳本代碼的執(zhí)行更快产捞。

目前作者封裝的59種TI方法列表見(jiàn)網(wǎng)址：https://github.com/dynverse/dynmethods#list-of-included-methods，用到代碼中的調(diào)用一般是單詞的所有字母都小寫(xiě)哼御，單詞之間加下劃線坯临。

guidelines推薦給我的methods_selected包含4個(gè)TI方法，本文開(kāi)頭介紹了《單細(xì)胞軌跡分析知多少--擬時(shí)間分析比較》得出了幾項(xiàng)重要結(jié)論恋昼，理論指導(dǎo)實(shí)踐看靠，所以我打算把4中TI方法都試一下，看看結(jié)果怎么樣液肌。

簡(jiǎn)單可視化

model <- model_slingshot
plot_dimred(
  model, 
  expression_source = dataset$expression, 
  grouping = dataset$grouping
)

slingshot

paga

mst

paga_tree

從生物學(xué)角度解釋軌跡

Rooting

大多數(shù)方法沒(méi)有直接的方法來(lái)推斷軌跡的方向性挟炬。在這種情況下，應(yīng)該使用一些外部信息來(lái)“確定”軌跡嗦哆，例如使用一組marker基因

添加rooting gene列表谤祖，這里需要根據(jù)具體的生物學(xué)問(wèn)題來(lái)定，為了更快速的演示吝秕，我不再深究泊脐，這里僅給一個(gè)示例作為演示用

# 初步來(lái)看我認(rèn)為slingshot方法得到的結(jié)果和我的預(yù)期比較相近，后續(xù)將以此方法得到的model進(jìn)行后續(xù)演示
model <- model_slingshot
model <- model %>% 
  add_root_using_expression(c("KLRB1"), dataset$expression)

Milestone labelling

可以使用marker基因標(biāo)記里程碑烁峭，這些標(biāo)簽隨后可用于后續(xù)分析和可視化

# 通過(guò)model$milestone_ids我們知道了有4個(gè)里程碑
# 這個(gè)示例我們標(biāo)記3個(gè)里程碑，此處僅為示例，請(qǐng)忽略標(biāo)簽意義
model <- label_milestones_markers(
  model,
  markers = list(
    Tn = c("CCR7"),
    NK = c("KLRD1"),
    CTL = c("GZMK")
  ),
  dataset$expression
)

軌跡可視化

dyno提供了降維约郁、軌跡模型和細(xì)胞聚類等可視化展示方法

細(xì)胞亞群在軌跡上的展示

model <- model %>% add_dimred(dyndimred::dimred_mds, expression_source = dataset$expression)
plot_dimred(
  model, 
  expression_source = dataset$expression, 
  grouping = dataset$grouping,
  label_milestones = TRUE
)

通過(guò)基礎(chǔ)的4個(gè)細(xì)胞亞群的軌跡分布缩挑，我們明顯的看到3亞群分為了2部分，所以我打算后續(xù)調(diào)整第3群為2個(gè)亞群再做軌跡分析鬓梅，以便更好的體現(xiàn)它們的差異供置。

我們?cè)O(shè)置了label_milestones = TRUE，所以里程碑被標(biāo)注到了軌跡模型上绽快，model一共有4個(gè)里程碑我們只標(biāo)注了3個(gè)所以上圖只顯示3個(gè)芥丧。

基因表達(dá)在軌跡上的展示

這里我們做CCR7和CCL5兩個(gè)基因表達(dá)示例

plot_dimred(
  model, 
  expression_source = dataset$expression, 
  color_cells = "feature",
  feature_oi = "CCR7",
  color_density = "grouping",
  grouping = dataset$grouping,
  label_milestones = TRUE
)

在展示基因表達(dá)的同時(shí)，不同細(xì)胞亞群以細(xì)胞的背景顏色表示坊罢，以下官方示例展示的更清楚一些

預(yù)測(cè)和可視化感興趣的基因

dyno整合了幾種方法來(lái)從軌跡中提取候選標(biāo)記基因/特征

最具預(yù)測(cè)性的基因的全球概述

默認(rèn)情況下续担，繪制熱圖時(shí)會(huì)計(jì)算總體上最重要的基因

plot_heatmap(
  model,
  expression_source = dataset$expression,
  grouping = dataset$grouping,
  features_oi = 50
)

Lineage/branch markers

dyno還有提取特定分支的功能，例如：當(dāng)細(xì)胞分化成神經(jīng)元時(shí)改變的基因

branch_feature_importance <- calculate_branch_feature_importance(model, expression_source=dataset$expression)

select_features <- branch_feature_importance %>% 
  filter(to == which(model$milestone_labelling =="CTL")) %>% 
  top_n(30, importance) %>% 
  pull(feature_id)
  
plot_heatmap(
  model, 
  expression_source = dataset$expression, 
  features_oi = select_features
)

分叉點(diǎn)重要的基因

提取在分支點(diǎn)發(fā)生變化的特征/基因

branching_milestone <- model$milestone_network %>% group_by(from) %>% filter(n() > 1) %>% pull(from) %>% first()

branch_feature_importance <- calculate_branching_point_feature_importance(model, expression_source=dataset$expression, milestones_oi = branching_milestone)

branching_point_features <- branch_feature_importance %>% top_n(20, importance) %>% pull(feature_id)

plot_heatmap(
  model,
  expression_source = dataset$expression,
  features_oi = branching_point_features
)

space <- dyndimred::dimred_mds(dataset$expression)
map(branching_point_features[1:12], function(feature_oi) {
  plot_dimred(model, dimred = space, expression_source = dataset$expression, feature_oi = feature_oi, label_milestones = FALSE) +
    theme(legend.position = "none") +
    ggtitle(feature_oi)
}) %>% patchwork::wrap_plots()

快速免安裝體驗(yàn)方式

對(duì)于MACOS或者有root權(quán)限的用戶活孩，如果想快速體驗(yàn)dyno物遇，可以在安裝并啟動(dòng)docker之后，執(zhí)行以下命令行進(jìn)入docker環(huán)境命令行憾儒，然后打開(kāi)R終端询兴。

sudo docker run -it --memory=20G --memory-swap=20G --oom-kill-disable -v /home/zanyuan:/zanyuan -v /tmp:/tmp -v /usr/bin/:/dockerbin -v /var/run/docker.sock:/var/run/docker.sock seqyuan/seqyuan-r:v0.0.1 /bin/bash

理解上面的代碼可能需要一些docker方面的知識(shí)，其中-v /home/zanyuan:/zanyuan 這一行指的是本地的存儲(chǔ)掛載為docker容器里的/zanyuan地址起趾；-v /tmp:/tmp是臨時(shí)存儲(chǔ)目錄掛載诗舰；-v /usr/bin/:/dockerbin指的是docker所在本地路徑掛載為容器里的/dockerbin；seqyuan/seqyuan-r:v0.0.1這個(gè)是docker啟動(dòng)的鏡像名稱训裆，我構(gòu)建的dockerfile請(qǐng)見(jiàn)docker hub官網(wǎng)始衅。

如果運(yùn)行時(shí)間比較長(zhǎng)，可以把代碼寫(xiě)到一個(gè)名為test_dyno.R腳本里例如：

library(dyno)
library(tidyverse)
library(Matrix)
library(Seurat)

pso <- readRDS(file = "/zanyuan/data/cd8_t1.rds")

dataset <- wrap_expression(
  counts = t(pso@assays$RNA@counts),
  expression = t(pso@assays$RNA@data)
)


dataset <- add_prior_information(
  dataset,
  start_id = "AGCGTCGTCCCTAACC.5"
)

model <- infer_trajectory(dataset, "paga")

saveRDS(model, file = "/zanyuan/data/cd8_paga.model", ascii = FALSE, version = NULL, compress = TRUE, refhook = NULL)

然后執(zhí)行

sudo docker run -it --memory=20G --memory-swap=20G --oom-kill-disable -v /home/zanyuan:/zanyuan -v /tmp:/tmp -v /usr/bin/:/dockerbin -v /var/run/docker.sock:/var/run/docker.sock seqyuan/seqyuan-r:v0.0.1 Rscript /zanyuan/test_ti.R

后續(xù)

后續(xù)會(huì)推出在linux集群無(wú)root權(quán)限的情況下調(diào)用singularity來(lái)執(zhí)行dyno