探索Bioconductor 中单细胞数据分析工具

在过去的几年里，分析scRNA-seq的方法数量激增，现在有超过200个可用的软件工具。每个工具都需要选择如何存储和表示分析期间使用的数据。在Keegan Korthauer的帮助下，Davide Risso和Aaron Lun创建了singlecellexper包，试图将所使用的数据结构标准化。该软件包作为Bioconductor 3.6版本的一部分于2017年10月发布。由于我们最近有了另一个Bioconductor 版本，我想我应该看看社区已经开发了哪些围绕着单细胞实验的工具。

What is a SingleCellExperiment?

在我们看什么包装使用SingleCellExperiment 之前，简要地讨论一下它是什么可能是有用的。SingleCellExperiment 对象是较老的SummarizedExperiment对象的扩展。这是一个S4类开发用于Bioconductor 包，其主要部分是一组中央矩阵“assays”，以及提供关于行和列的额外信息的表。还有一个元数据槽，它是一个包含与实验相关的任何其他信息的列表。

Diagram of a SummarizedExperiment object

使用SummarizedExperiment这样的结构的一个关键好处是，与分析相关的所有数据都保存在一个地方。这使得在函数或输出结果之间传递信息变得更容易，并减少了不匹配的可能性。该SingleCellExperiment 增加了一些额外的功能，这是有用的scRNA-seq分析，包括:

• Slots for holding:
• Dimenstionality reductions
• Spike-in information
• Size factors for normalisation
• Convenient access for named assays - counts, normcounts, cpm etc.

这个想法是这样的，SingleCellExperiment 可以被一系列的程序包用来在分析过程中存储数据，并在需要的时候进行扩展。让我们看看这些包目前是什么，看看什么取决于SingleCellExperiment 。

Getting package information

首先，让我们加载分析所需的包。BiocPkgTools需要从GitHub上安装，但其余都在CRAN上。

library("BiocPkgTools") # https://github.com/seandavi/BiocPkgTools
library("tidygraph")
library("ggraph")
library("tidyverse")

我们可以使用Sean Davis的BiocPkgTools包来获得关于Bioconductor包的信息，但是我们需要做一些过滤来获得我们想要的信息。我需要做很多次所以这是我写的一个函数让它更简单。它获取关于所有Bioconductor包的信息的数据库、我们感兴趣的包的名称以及一个表示我们是否需要正常或反向依赖关系的标志。

get_bioc_deps <- function(bpi, pkg, reverse) {
    deps <- bpi %>%
        filter(Package == pkg)
        
    if (reverse) {
        deps <- deps %>%
            select(depends = dependsOnMe, imports = importsMe,
                   suggests = suggestsMe)
    } else {
        deps <- deps %>%
            select(depends = Depends, imports = Imports,
                   suggests = Suggests)
    }
    
    deps <- deps %>%
        gather(key = "type", value = "package") %>%
        separate_rows() %>%
        filter(!is.na(package))
    
    if (reverse) {
        deps <- deps %>%
            mutate(package2 = pkg) %>%
            rename(package1 = package)
    } else {
        deps <- deps %>%
            mutate(package1 = pkg) %>%
            rename(package2 = package)
    }
    
    deps <- deps %>% select(package1, uses = type, package2)
}

如果我们用它来搜索SingleCellExperiment 的反向依赖关系，我们可以看到它返回一个数据帧，其中Bioconductor 包使用SingleCellExperiment 和它们之间的关系。

bpi = biocPkgList()
bioc_revdeps <- get_bioc_deps(bpi, "SingleCellExperiment", reverse = TRUE)
head(bioc_revdeps)

# A tibble: 3 x 3
  package1   uses     package2            
  <list>     <chr>    <chr>               
1 <chr [15]> depends  SingleCellExperiment
2 <chr [15]> imports  SingleCellExperiment
3 <chr [5]>  suggests SingleCellExperiment

我比较好奇bpi 是什么：

 bpi
# A tibble: 1,649 x 45
   Package Version Depends Suggests License MD5sum NeedsCompilation Title Description biocViews Author Maintainer git_url git_branch
   <chr>   <chr>   <list>  <list>   <chr>   <chr>  <chr>            <chr> <chr>       <list>    <list> <list>     <chr>   <chr>     
 1 a4      1.30.0  <chr [~ <chr [4~ GPL-3   4e44a~ no               Auto~ Automated ~ <chr [1]> <chr ~ <chr [2]>  https:~ RELEASE_3~
 2 a4Base  1.30.0  <chr [~ <chr [2~ GPL-3   a8238~ no               Auto~ Automated ~ <chr [1]> <chr ~ <chr [2]>  https:~ RELEASE_3~
 3 a4Clas~ 1.30.0  <chr [~ <chr [1~ GPL-3   f3c57~ no               Auto~ Automated ~ <chr [1]> <chr ~ <chr [2]>  https:~ RELEASE_3~
 4 a4Core  1.30.0  <chr [~ <chr [1~ GPL-3   1c18e~ no               Auto~ Automated ~ <chr [1]> <chr ~ <chr [2]>  https:~ RELEASE_3~
 5 a4Prep~ 1.30.0  <chr [~ <chr [2~ GPL-3   79f08~ no               Auto~ Automated ~ <chr [1]> <chr ~ <chr [2]>  https:~ RELEASE_3~
 6 a4Repo~ 1.30.0  <chr [~ <chr [1~ GPL-3   4a7eb~ no               Auto~ Automated ~ <chr [1]> <chr ~ <chr [2]>  https:~ RELEASE_3~
 7 ABAEnr~ 1.12.0  <chr [~ <chr [3~ GPL (>~ f5934~ yes              Gene~ "The packa~ <chr [2]> <chr ~ <chr [1]>  https:~ RELEASE_3~
 8 ABarray 1.50.0  <chr [~ <chr [2~ GPL     84761~ no               "Mic~ "Automated~ <chr [3]> <chr ~ <chr [1]>  https:~ RELEASE_3~
 9 abseqR  1.0.0   <chr [~ <chr [1~ GPL-3 ~ de9e8~ no               "Rep~ "AbSeq is ~ <chr [5]> <chr ~ <chr [1]>  https:~ RELEASE_3~
10 ABSSeq  1.36.0  <chr [~ <chr [1~ GPL (>~ 03058~ no               "ABS~ "Inferring~ <chr [1]> <chr ~ <chr [1]>  https:~ RELEASE_3~
# ... with 1,639 more rows, and 31 more variables: git_last_commit <chr>, git_last_commit_date <chr>, `Date/Publication` <chr>,
#   source.ver <chr>, win.binary.ver <chr>, `mac.binary.el-capitan.ver` <chr>, vignettes <list>, vignetteTitles <list>, hasREADME <chr>,
#   hasNEWS <chr>, hasINSTALL <chr>, hasLICENSE <chr>, Rfiles <chr>, Enhances <list>, dependsOnMe <list>, Imports <list>, importsMe <list>,
#   suggestsMe <list>, LinkingTo <chr>, Archs <chr>, VignetteBuilder <chr>, URL <chr>, SystemRequirements <chr>, BugReports <chr>,
#   PackageStatus <chr>, Video <chr>, linksToMe <chr>, OS_type <chr>, License_restricts_use <chr>, License_is_FOSS <chr>, organism <chr>

我们可以使用工具::package_dependencies函数对CRAN包执行类似的操作。

get_cran_deps <- function(pkg, db, reverse) {
    
    types <- c("Depends", "Imports", "Suggests")
    
    deps <- sapply(types, function(type) {
        deps <- tools::package_dependencies(pkg, db, which = type,
                                            reverse = reverse)
        c(type = type, package = paste(deps[[1]], collapse = ", "))
    })
    
    deps <- deps %>%
        t() %>%
        as_data_frame() %>%
        mutate(type = tolower(type)) %>%
        filter(package != "") %>%
        separate_rows(package)
    
    if (nrow(deps) == 0) {
        return(tibble(package1 = character(), uses = character(),
                      package2 = character()))
    }
    
    if (reverse) {
        deps <- deps %>%
            mutate(package2 = pkg) %>%
            rename(package1 = package)
    } else {
        deps <- deps %>%
            mutate(package1 = pkg) %>%
            rename(package2 = package)
    }
    
    deps <- deps %>% select(package1, uses = type, package2)
}

db <- available.packages(repos = "http://cran.r-project.org")
cran_revdeps <- get_cran_deps("SingleCellExperiment", db, reverse = TRUE)
head(cran_revdeps)

# A tibble: 2 x 3
  package1 uses     package2            
  <chr>    <chr>    <chr>               
1 clustree suggests SingleCellExperiment
2 Seurat   suggests SingleCellExperiment

What uses SingleCellExperiment?

现在我们有两张表告诉我们哪些生物导体和CRAN包装使用了单细胞搽剂。不过，查看表可能相当乏味，因此让我们使用tidygraph使用这些关系来构建一个图。

package1=unlist(bioc_revdeps$package1)
uses=c(rep("depends",15),rep("imports",15),rep("suggests",5))
package2<-rep("SingleCellExperiment",35)
mybioc_revdeps<-data.frame(package1,uses,package2)

nodes <- mybioc_revdeps %>%
  bind_rows(cran_revdeps) %>%
  select(-uses) %>%
  gather(key = id, value = package) %>%
  select(-id) %>%
  distinct() %>%
  mutate(repo = if_else(package %in% bpi$Package, "Bioconductor", "CRAN"))

edges <- mybioc_revdeps %>%
  bind_rows(cran_revdeps) %>%
  rename(from = package1, to = package2)

graph <- tbl_graph(nodes = nodes, edges = edges)

我们现在可以用ggraph可视化这些关系:

ggraph(graph, layout = "fr") +
    geom_edge_fan(aes(colour = uses),
                  arrow = arrow(length = unit(4, 'mm')), 
                  end_cap = circle(3, 'mm')) +
    geom_node_point(aes(colour = repo)) +
    geom_node_text(aes(label = package, colour = repo), repel = TRUE) +
    scale_color_brewer(palette = "Set1") +
    scale_edge_color_brewer(palette = "Dark2") +
    theme_graph()

这并没有告诉我们很多我们不知道的东西，但它确实让我们可以在一个地方看到所有的东西。我们可以看到只有几个CRAN包，这并不奇怪，因为SingleCellExperiment是Bioconductor的一部分，而且大多数包要么“import”要么“depend”SingleCellExperiment。

依赖于SingleCellExperiment包之间的关系如何?是否有相关的scRNA-seq分析工具?

Adding an extra hop

我们可以重用前面编写的函数来获得我们的scRNA-seq包列表的依赖项(和反向依赖项)。我们还将做一些额外的处理来整理一些结果。

它还使您的使用变得容易得多，因为您不需要学习新的数据结构来使用包，并且可以使用一系列包，而不必在对象之间进行转换。如果您不想被限制在Bioconductor生态系统中，可以考虑使用Seurat对象，或者如果您使用Python，可以考虑使用anndata对象。还有一种loom
格式，它同时具有R和Python接口。如果社区能够利用少量的数据结构，而不是每个包都使用自己的数据结构，那么无论哪种标准适合您，每个人都会受益。

---

---

exploring-the-sce-verse