本文介绍了如何利用R语言中的liana包，将流行的细胞通讯分析软件如CellPhoneDB、CellChat等整合起来，对单细胞数据进行分析。作者展示了如何调用这些方法，以及如何通过liana获得多软件和数据库的共识排名。
hello，今天我们分享一个简单的，将最近流行的细胞通讯分析软件整合在一起，分析我们的单细胞数据，直接看看代码部分
注意一个名词， cell-cell communication (CCC)，细胞通讯的软件方法很多了，these methods and resources are usually in a fixed combination of a tool and its corresponding resource, but in principle any resource could be combined with any method.（随意结合）。
图片.png
加载包
require(liana)
require(tibble)
require(purrr)
首先看看数据库来源
# Resource currently included in OmniPathR (and hence `liana`) include:
get_resources()
#> [1] "Reshuffled" "Default" "CellChatDB" "CellPhoneDB"
#> [5] "Ramilowski2015" "Baccin2019" "LRdb" "Kirouac2010"
#> [9] "ICELLNET" "iTALK" "EMBRACE" "HPMR"
#> [13] "Guide2Pharma" "connectomeDB2020" "talklr" "CellTalkDB"
#> [17] "OmniPath"
# A list of resources can be obtained using the `select_resource()` function:
# See `?select_resource()` documentation for further information.
# select_resource(c('OmniPath')) %>% glimpse()
来源是真地多。
CCC Methods
- CellPhoneDB algorithm (via Squidpy)
- CellChat
- NATMI
- Connectome
- SingleCellSignalR (SCA)
- iTALK
目前就支持这么多方法，但是后续本人会再往上添加分析方法。
调用包装函数 liana
liana_path <- system.file(package = "liana")
testdata <-
readRDS(file.path(liana_path , "testdata", "input", "testdata.rds"))
testdata %>% glimpse()
#> Formal class 'Seurat' [package "Seurat"] with 13 slots
#> ..@ assays :List of 1
#> .. ..$ RNA:Formal class 'Assay' [package "Seurat"] with 8 slots
#> ..@ meta.data :'data.frame': 90 obs. of 4 variables:
#> .. ..$ orig.ident : Factor w/ 1 level "pbmc3k": 1 1 1 1 1 1 1 1 1 1 ...
#> .. ..$ nCount_RNA : num [1:90] 4903 3914 4973 3281 2641 ...
#> .. ..$ nFeature_RNA : int [1:90] 1352 1112 1445 1015 928 937 899 1713 960 888 ...
#> .. ..$ seurat_annotations: Factor w/ 3 levels "B","CD8 T","NK": 1 1 1 1 3 3 1 1 1 1 ...
#> ..@ active.assay: chr "RNA"
#> ..@ active.ident: Factor w/ 3 levels "B","CD8 T","NK": 1 1 1 1 3 3 1 1 1 1 ...
#> .. ..- attr(*, "names")= chr [1:90] "AAACATTGAGCTAC" "AAACTTGAAAAACG" "AAAGGCCTGTCTAG" "AAAGTTTGGGGTGA" ...
#> ..@ graphs : list()
#> ..@ neighbors : list()
#> ..@ reductions : list()
#> ..@ images : list()
#> ..@ project.name: chr "SeuratProject"
#> ..@ misc : list()
#> ..@ version :Classes 'package_version', 'numeric_version' hidden list of 1
#> .. ..$ : int [1:3] 3 2 3
#> ..@ commands :List of 1
#> .. ..$ FindVariableFeatures.RNA:Formal class 'SeuratCommand' [package "Seurat"] with 5 slots
#> ..@ tools : list()
liana_wrap 调用了许多方法，并且每个方法都使用提供的资源运行。
在这里，我们为了展示调用 Connectome 和 Squidpy，但是任何 c(“sca”, “natmi”, “italk”, “cellchat”) 也可以通过简单地将它们的名称添加到方法参数来运行。
同样，这里我们仅使用“OmniPath”CCC 资源，但可以将上述任何资源（可通过 get_resources() 获得）添加到资源参数中 。
# Run liana
liana_test <- liana_wrap(testdata,
method = c('squidpy', 'connectome'),
resource = c('OmniPath'))
#> Warning in data.frame(source = sources[i], target = targets[j], ligand =
#> ligands, : row names were found from a short variable and have been discarded
#> Warning in data.frame(source = sources[i], target = targets[j], ligand =
#> ligands, : row names were found from a short variable and have been discarded
#> Warning in data.frame(source = sources[i], target = targets[j], ligand =
#> ligands, : row names were found from a short variable and have been discarded
#> Warning in data.frame(source = sources[i], target = targets[j], ligand =
#> ligands, : row names were found from a short variable and have been discarded
#> Warning in data.frame(source = sources[i], target = targets[j], ligand =
#> ligands, : row names were found from a short variable and have been discarded
#> Warning in data.frame(source = sources[i], target = targets[j], ligand =
#> ligands, : row names were found from a short variable and have been discarded
#> Warning in data.frame(source = sources[i], target = targets[j], ligand =
#> ligands, : row names were found from a short variable and have been discarded
#> Warning in data.frame(source = sources[i], target = targets[j], ligand =
#> ligands, : row names were found from a short variable and have been discarded
#> Warning in data.frame(source = sources[i], target = targets[j], ligand =
#> ligands, : row names were found from a short variable and have been discarded
获得共识排名(多软件，多数据库共同支持的显著配受体对)
liana 还为使用不同方法获得的结果提供了共识等级。 默认情况下，liana将提供均值、中位数和聚合*共识等级。
# Liana returns a list of results, each element of which corresponds to a method
liana_test
#> $squidpy
#> # A tibble: 9,396 x 8
#> ligand receptor source target means pvalue uniprot_source unprot_target
#> <chr> <chr> <chr> <chr> <dbl> <dbl> <chr> <chr>
#> 1 NRP2 SEMA4F B B 0 NaN O60462 O95754
#> 2 NRP2 SEMA4F B CD8 T 0 NaN O60462 O95754
#> 3 NRP2 SEMA4F B NK 0 NaN O60462 O95754
#> 4 NRP2 SEMA4F CD8 T B 0 NaN O60462 O95754
#> 5 NRP2 SEMA4F CD8 T CD8 T 0 NaN O60462 O95754
#> 6 NRP2 SEMA4F CD8 T NK 0 NaN O60462 O95754
#> 7 NRP2 SEMA4F NK B 0 NaN O60462 O95754
#> 8 NRP2 SEMA4F NK CD8 T 0 NaN O60462 O95754
#> 9 NRP2 SEMA4F NK NK 0 NaN O60462 O95754
#> 10 BAMBI BMPR1A B B 0 NaN Q13145 P36894
#> # … with 9,386 more rows
#>
#> $connectome
#> # A tibble: 9,396 x 8
#> source target ligand receptor weight_norm weight_sc p_val_adj.lig
#> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 CD8 T B BMPR2 ACVR1 0 0.0375 1
#> 2 B B BMPR2 ACVR1 0 0.0375 1
#> 3 NK B BMPR2 ACVR1 0 0.0375 1
#> 4 CD8 T CD8 T BMPR2 ACVR1 0 -0.0750 1
#> 5 B CD8 T BMPR2 ACVR1 0 -0.0750 1
#> 6 NK CD8 T BMPR2 ACVR1 0 -0.0750 1
#> 7 CD8 T NK BMPR2 ACVR1 0 0.0375 1
#> 8 NK NK BMPR2 ACVR1 0 0.0375 1
#> 9 B NK BMPR2 ACVR1 0 0.0375 1
#> 10 B B GDF11 ACVR1B 0 0.105 1
#> # … with 9,386 more rows, and 1 more variable: p_val_adj.rec <dbl>
# We can aggregate these results into a dataframe with consensus ranks
liana_test %>%
liana_aggregate()
#> # A tibble: 9,396 x 11
#> source ligand target receptor aggregate_rank mean_rank median_rank
#> <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
#> 1 NK HLA-B CD8 T CD8A 0.00000999 11 11
#> 2 NK B2M NK CD247 0.00000999 4 4
#> 3 NK B2M NK KLRD1 0.0000191 5 5
#> 4 NK B2M NK KLRC1 0.0000522 24.5 24.5
#> 5 NK B2M NK KIR2DL3 0.000102 34 34
#> 6 NK HLA-B NK KIR2DL3 0.000108 58.5 58.5
#> 7 B LGALS9 NK HAVCR2 0.000179 45 45
#> 8 NK HLA-B NK KIR3DL2 0.000192 70 70
#> 9 NK HLA-G NK KLRC1 0.000196 70.5 70.5
#> 10 CD8 T PTPRC B CD22 0.000236 75 75
#> # … with 9,386 more rows, and 4 more variables: squidpy.pvalue <dbl>,
#> # squidpy.rank <dbl>, connectome.weight_sc <dbl>, connectome.rank <dbl>
使用 RobustRankAggreg 包获得聚合共识等级。 RobustRankAggreg 包提供了不同的函数来获取聚合排名，默认情况下使用 RRA 方法。 RRA 提供可以解释为 p 值的等级分数。 排名始终高于随机的交互被分配低分数p 值。
使用覆盖的默认设置调用liana
默认情况下，liana 使用每个方法的默认值运行，可以通过 liana_default() 获得
或者，也可以通过简单地将它们传递给 liana 包装函数来覆盖默认设置
# Overwrite CellChat default parameters by providing a list of parameters
liana_test <- liana_wrap(testdata,
method = c('cellchat', 'squidpy'),
resource = 'Default', # Run with Default Resources
cellchat.params =
list(
nboot = 10,
exclude_anns = NULL,
thresh = 1,
assay = "RNA",
.normalize = TRUE,
.do_parallel = FALSE,
.raw_use = TRUE
)
)
#> The cell groups used for CellChat analysis are B CD8 T NK
# This returns a list of results for each method
liana_test %>% glimpse
#> List of 2
#> $ cellchat: tibble [56 × 6] (S3: tbl_df/tbl/data.frame)
#> ..$ source : Factor w/ 3 levels "B","CD8 T","NK": 1 2 3 1 2 3 1 2 3 1 ...
#> ..$ target : Factor w/ 3 levels "B","CD8 T","NK": 1 1 1 2 2 2 3 3 3 1 ...
#> ..$ ligand : chr [1:56] "MIF" "MIF" "MIF" "MIF" ...
#> ..$ receptor: chr [1:56] "CD74_CXCR4" "CD74_CXCR4" "CD74_CXCR4" "CD74_CXCR4" ...
#> ..$ prob : num [1:56] 0.123 0.151 0.166 0.111 0.137 ...
#> ..$ pval : num [1:56] 0.7 0.1 0 0.8 0.2 0 0.9 0.9 0.8 0.8 ...
#> $ squidpy : Named list()
单独使用一个方法
# RUN cellchat alone with OmniPath
cc_res <- call_cellchat(
seurat_object = testdata,
op_resource = select_resource('OmniPath')[[1]],
.normalize = TRUE,
nboot = 10)
#> Create a CellChat object from a data matrix
#> Set cell identities for the new CellChat object
#> The cell groups used for CellChat analysis are B CD8 T NK
#> Issue identified!! Please check the official Gene Symbol of the following genes:
#> CYR61 WISP1 CTGF FAM19A4 CD97
# Show CellChat Results
cc_res
#> # A tibble: 70 x 6
#> source target ligand receptor prob pval
#> <fct> <fct> <chr> <chr> <dbl> <dbl>
#> 1 NK NK TYROBP KLRD1 0.144 0
#> 2 B CD8 T HLA-B CD8A 0.154 0
#> 3 CD8 T CD8 T HLA-B CD8A 0.169 0
#> 4 NK CD8 T HLA-B CD8A 0.183 0
#> 5 B CD8 T HLA-B CD3D 0.266 0
#> 6 CD8 T CD8 T HLA-B CD3D 0.288 0
#> 7 NK CD8 T HLA-B CD3D 0.307 0
#> 8 B NK HLA-B KLRD1 0.171 0
#> 9 CD8 T NK HLA-B KLRD1 0.187 0
#> 10 NK NK HLA-B KLRD1 0.201 0
#> # … with 60 more rows
本人呢，其实能力还有限，不过还是希望大家支持我一下，今后肯定会联合更多的软件来开发一些好用的包和函数，甚至于，开发出我自己创新性的R包供大家使用，希望大家支持我一下，如果有志同道合的人，多多交流，多多联系。
生活很好，有你更好~~~
扫一扫
1429
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
