最近看到Fly Cell Atlas (FCA) 公布果蠅的多個組織的單細胞核轉(zhuǎn)錄組測序數(shù)據(jù)猜谚，雖然FCA組織提供了多個在線平臺（例如SCope）進行非常方便的可視化和多種分析。然而源武，由于網(wǎng)速的原因，還是希望將數(shù)據(jù)下載到本地進行分析称诗。

Figure Source: SCope (https://scope.aertslab.org/#/FlyCellAtlas/)

FCA提供了.loom 和 .h5ad格式的數(shù)據(jù)致讥，而我比較習(xí)慣使用R的Seurat處理單細胞的數(shù)據(jù)富岳，便想看看有沒有方法將這兩種格式的文件轉(zhuǎn)換成Seurat對象。

實際上陈瘦，Seurat也開發(fā)了SeuratDisk這一軟件包支持將.loom轉(zhuǎn)換為Seurat對象（還支持多種單細胞數(shù)據(jù)格式的轉(zhuǎn)換）幌甘。

不巧的是，F(xiàn)CA的數(shù)據(jù)存儲格式似乎與Seurat設(shè)置的格式不一致痊项，導(dǎo)致簡單的數(shù)據(jù)轉(zhuǎn)換變得十分復(fù)雜...因此锅风，本文記錄我在轉(zhuǎn)換.loom/.h5ad為Seurat對象時踩過的坑。

數(shù)據(jù)下載

FCA的測序數(shù)據(jù)都可以在該網(wǎng)站下載：https://flycellatlas.org/

我這篇文章先嘗試了使用果蠅頭部的數(shù)據(jù)進行轉(zhuǎn)換

Figure Source: Fly Cell Atlas (https://flycellatlas.org/)

loom 文件

.loom格式是一種用于處理大型單細胞數(shù)據(jù)的文件格式鞍泉，它并不直接將單細胞數(shù)據(jù)讀入（R中）皱埠，而是與文件創(chuàng)建鏈接。這樣可以只在使用到數(shù)據(jù)時咖驮，讀入其中的一部分边器，避免內(nèi)存的過度占用（個人理解）训枢。

Figure Source: Loompy (http://linnarssonlab.org/loompy/index.html)

FCA 中使用的是3.0.0版本的loom，以下我只對這個版本的格式進行介紹忘巧。

.loom實質(zhì)上是基于HDF5 (Hierarchical Data Format)格式的恒界。所以，我們可以使用HDFView 這種查看HDF格式的軟件觀察其中的數(shù)據(jù)格式砚嘴。

HDFView可在https://www.hdfgroup.org/downloads/hdf5/ 下載

郵箱注冊之后選擇相應(yīng)系統(tǒng)版本下載就行十酣。

loom v3的文件包含了以下幾個層級的數(shù)據(jù)：

• /matrix(must): 基因*細胞的矩陣
• /attrs (must): 全局的屬性信息，例如使用的參考基因組际长、數(shù)據(jù)創(chuàng)建日期等信息
• /col_attrs (must): 列（細胞）的屬性耸采，例如barcode, cell id, cell annotation

• /row_attrs (must): 行（基因）的屬性，包括基因ID也颤，細胞marker之類的

• /layers (optional): 另外的基因*細胞矩陣洋幻，例如可以保存Seurat的data以及scale.data

• /col_graphs (optional): 細胞的圖信息，例如最近鄰圖（nearest-neighbor graphs）

• /row_graphs (optional): 基因的圖信息

loom to Seurat

簡單了解loom文件格式后翅娶，我們就可以在R中著手將其轉(zhuǎn)換為Seurat對象

SeuratDisk提供了在R中處理 loom v3文件所需的函數(shù)文留。

目前SeuratDisk還只能安裝其Github上的版本

remotes::install_github(repo = 'mojaveazure/seurat-disk', ref = 'develop')

創(chuàng)建loom鏈接

library(SeuratDisk)
library(SeuratObject)
library(Seurat)
# Connect to the file
# Loom v3 standards
ds <- Connect(filename = "data/ref/s_fca_biohub_head_10x.loom", mode = "r") # r for read-only and r+ for read/write

> ds
Class: loom
Filename: ~\data\ref\s_fca_biohub_head_10x.loom
Access type: H5F_ACC_RDONLY
Attributes: last_modified
Listing:
       name    obj_type   dataset.dims dataset.type_class
      attrs   H5I_GROUP           <NA>               <NA>
  col_attrs   H5I_GROUP           <NA>               <NA>
 col_graphs   H5I_GROUP           <NA>               <NA>
     layers   H5I_GROUP           <NA>               <NA>
     matrix H5I_DATASET 100527 x 13056          H5T_FLOAT
  row_attrs   H5I_GROUP           <NA>               <NA>
 row_graphs   H5I_GROUP           <NA>               <NA>

訪問loom中的數(shù)據(jù)

可以通過[[的方法訪問loom文件中的數(shù)據(jù)，注意到我們是通過一種路徑的方式來訪問其中的數(shù)據(jù)

# gene * cell matrix 
> ds[["/matrix"]]
Class: H5D
Dataset: /matrix
Filename: ~\data\ref\s_fca_biohub_head_10x.loom
Access type: H5F_ACC_RDONLY
Datatype: H5T_IEEE_F32LE
Space: Type=Simple     Dims=100527 x 13056     Maxdims=Inf x 13056
Chunk: 64 x 64

> ds[["/matrix"]][1:3,1:3]
     [,1] [,2] [,3]
[1,]    0    7   19
[2,]    0    1    2
[3,]    0    1    4

> ds[["/matrix"]]$dims
[1] 100527  13056

# cell annotation
> ds[["/col_attrs/annotation"]][1:3]
[1] "unannotated"                "unannotated"                "lamina wide-field 2 neuron"

在R中竭沫，HDF5格式的文件被轉(zhuǎn)換為H5D類的對象燥翅，類似于R6對象，H5D也定義了一系列的方法

> ds[['/matrix']]$methods()
<H5D>
  Inherits from: <H5RefClass>
  Public:
    attr_delete: function (attr_name) 
    attr_delete_by_idx: function (n, obj_name, idx_type = h5const$H5_INDEX_NAME, order = h5const$H5_ITER_NATIVE, 
...
  Private:
    closeFun: function (id) 
    get_robj_dim: function (reg_eval_res) 
    pid: environment
    pmessage: 
    stopOnInvalid: function (error.msg = "The object is invalid")

由于我沒有深入去探究H5D類的各種方法是如何調(diào)用蜕提，這里就不做展開了森书，有興趣的朋友可以查看Rhdf5r包的文檔：

https://www.rdocumentation.org/packages/hdf5r/versions/1.3.3/topics/H5D-class

Converting from loom to seurat

接下來，我嘗試通過Seurat::as.Seurat()將loom格式轉(zhuǎn)換為Seurat對象谎势。但在運行幾個步驟后就報錯了...

> fca_head_seurat <- as.Seurat(ds)
# Reading in /matrix
# Storing /matrix as counts
# Saving /matrix to assay 'RNA'
# Error in `[<-.data.frame`(`*tmp*`, names.intersect, i[index], value = c(`128up` = 0,  : 
#    missing values are not allowed in subscripted assignments of data frames

搜索該報錯凛膏，可以找到seurat github上的一個issue有討論這種錯誤。

https://github.com/satijalab/seurat/issues/5124

大意是FCA存儲loom文件的格式跟Seurat識別的格式不太一樣脏榆，所以就報錯了猖毫。目前這種格式對應(yīng)的問題還沒解決，只能通過使用H5AD格式的文件轉(zhuǎn)換為Seurat對象來繞過這個問題须喂。

H5AD

H5AD 也是一種HDF5格式的文件吁断，其中包含額外的AnnData

Figure Source: AnnData (https://anndata.readthedocs.io/en/latest/)

這里不再過多介紹H5AD，我們大概知道它也是用于存儲單細胞數(shù)據(jù)就可以了坞生，有興趣的朋友可以自行了解：https://anndata.readthedocs.io/en/latest/

Converting h5ad to seurat

接下來仔役，我們將H5AD轉(zhuǎn)換為h5seurat格式，再讀入R中就是Seurat對象了

library(Seurat)
library(SeuratData)
library(SeuratDisk)
# Converting from AnnData to Seurat via h5Seurat
Convert("data/ref/s_fca_biohub_head_10x.h5ad", dest="h5seurat", overwrite=T)
fca_head <- LoadH5Seurat("data/ref/s_fca_biohub_head_10x.h5seurat")

> fca_head
An object of class Seurat 
1838 features across 100527 samples within 1 assay 
Active assay: RNA (1838 features, 0 variable features)
 3 dimensional reductions calculated: pca, tsne, umap

這次轉(zhuǎn)換成功了是己，我們獲得Seurat對象 ??

但是又兵，注意到h5ad轉(zhuǎn)換過來的對象只有1838個基因，而不是loom文件中的13056個卒废。有可能是因為這個數(shù)據(jù)中的基因是經(jīng)過一定過濾的沛厨。

讀入的這個對象居然有10.1 GB這么大乘盼，難怪他們會用loom/h5ad的文件格式來進行分析

在讀入數(shù)據(jù)后，事情也還沒結(jié)束俄烁。因為我們需要的是它的細胞注釋信息绸栅。但是讀入的數(shù)據(jù)結(jié)構(gòu)中，有關(guān)細胞信息的meta.data居然是空的页屠。

> fca_head@meta.data
data frame with 0 columns and 100527 rows

仔細觀察其數(shù)據(jù)結(jié)構(gòu)發(fā)現(xiàn)粹胯，相關(guān)的注釋信息都存儲在fca_head@misc當(dāng)中，但問題在于找不到一個注釋數(shù)據(jù)與細胞ID是匹配的辰企。換言之风纠，我也不能借助h5ad這組數(shù)據(jù)獲得細胞注釋。

Create seurat object

在轉(zhuǎn)換的辦法行不通后牢贸，自然想到直接使用loom文件和h5ad文件中的數(shù)據(jù)構(gòu)建seurat對象.

首先竹观，提取loom文件中的基因 X 細胞矩陣，注意將其轉(zhuǎn)換為稀疏矩陣（節(jié)省空間）潜索。由于提取出的矩陣是細胞 X 基因的形式臭增，故對其轉(zhuǎn)置

library(tidyverse)
library(SeuratDisk)
library(SeuratObject)
library(Seurat)
# gene * cell matrix 
fca_head_mat <- ds[["/matrix"]][,]
fca_head_mat <- Matrix::Matrix(fca_head_mat, sparse=T)
fca_head_mat <- Matrix::t(fca_head_mat)

接著，提取cell id and gene id竹习，并設(shè)置為矩陣的列名和行名

# cell id -- 100,527
fca_head_cellid <- ds[["/col_attrs/CellID"]][]
# gene -- 13,056
fca_head_geneid <- ds[["/row_attrs/Gene"]][]

colnames(fca_head_mat) <- fca_head_cellid
rownames(fca_head_mat) <- fca_head_geneid

提取我們需要的meta data. 由于我找不到一個H5D類的方法可以直接通過字符向量提取相應(yīng)meta data誊抛，這里用了一個循環(huán)的方法。

注意將meta.data的行名設(shè)置為與基因 X 細胞矩陣匹配的細胞ID

# pull metadata from loom
attrs <- c('Barcode', 'CellID', 'ClusterID', 'n_counts', 'n_genes', 'percent_mito',
           'annotation', 'annotation_broad', 'annotation_broad_extrapolated','sex')
# can't find method from loom to directly extract attributes
attrs_df <- map_dfc(attrs, ~ ds[[paste0("col_attrs/", .)]][]) %>% as.data.frame()
colnames(attrs_df) <- attrs
rownames(attrs_df) <- fca_head_cellid

構(gòu)建Seurat對象

# create seurat object from loom data
fca_head_seurat <- CreateSeuratObject(counts = fca_head_mat,
                                      meta.data = attrs_df)

又由于我找不到loom文件中降維的信息在哪一個section整陌，這里用h5ad的降維信息 ??

同時拗窃，注意到h5ad的細胞ID和loom的細胞ID是不一樣的，將其轉(zhuǎn)為loom的細胞ID

fca_head_h5 <- LoadH5Seurat("data/ref/s_fca_biohub_head_10x.h5seurat")
fca_head_seurat@reductions <- fca_head_h5@reductions
rownames(fca_head_seurat@reductions$pca@cell.embeddings) <- fca_head_cellid
rownames(fca_head_seurat@reductions$tsne@cell.embeddings) <- fca_head_cellid
rownames(fca_head_seurat@reductions$umap@cell.embeddings) <- fca_head_cellid

畫一個tsne plot 看看我們構(gòu)建的數(shù)據(jù)是否正確

DimPlot(fca_head_seurat, 
        reduction = 'tsne', 
        group.by = 'annotation', 
        label = T, 
        label.size = 4) + NoLegend()

label 太多看不清了泌辫，但整體結(jié)構(gòu)是一致的随夸，說明我們構(gòu)建成功了 ??

補充一點： Seurat團隊有提到Seurat相關(guān)包可以直接對loom文件進行分析。但他這里提到的版本是Seurat V2震放，而現(xiàn)在Seurat都是V4了宾毒，而V4的函數(shù)也不能直接使用loom對象作為輸入。

https://github.com/mojaveazure/loomR/issues/35

總結(jié)這次數(shù)據(jù)轉(zhuǎn)換的折騰過程澜搅，問題主要是由于我對單細胞數(shù)據(jù)伍俘、文件格式和分析流程的熟悉度不足夠所導(dǎo)致的邪锌。這讓我感覺到當(dāng)接觸一個新的領(lǐng)域時勉躺，我們都只能參考軟件文檔中的例子進行分析，超出例子范圍的操作就只能一個error觅丰、一個error的去搜索解決方案饵溅。這種debug的過程很痛苦，因為面對海量的信息妇萄，我們需要仔細地甄別哪一種是能為我們所用的蜕企，但這種debug的能力也是做生信必不可缺的咬荷。

最后，要說明我這里分享的是我探索的路徑轻掩，它并不一定是最優(yōu)的解決方案幸乒。

> sessionInfo()
R version 4.0.3 (2020-10-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19041)

Matrix products: default

locale:
[1] LC_COLLATE=Chinese (Simplified)_China.936  LC_CTYPE=Chinese (Simplified)_China.936   
[3] LC_MONETARY=Chinese (Simplified)_China.936 LC_NUMERIC=C                              
[5] LC_TIME=Chinese (Simplified)_China.936    

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] Seurat_4.0.5          SeuratObject_4.0.2    SeuratDisk_0.0.0.9019 forcats_0.5.0        
 [5] stringr_1.4.0         dplyr_1.0.2           purrr_0.3.4           readr_1.4.0          
 [9] tidyr_1.1.2           tibble_3.0.4          ggplot2_3.3.5         tidyverse_1.3.0      

loaded via a namespace (and not attached):
  [1] Rtsne_0.15            colorspace_2.0-0      deldir_1.0-2          ellipsis_0.3.1       
  [5] ggridges_0.5.2        fs_1.5.0              spatstat.data_2.1-0   rstudioapi_0.13      
  [9] leiden_0.3.6          listenv_0.8.0         bit64_4.0.5           ggrepel_0.9.0        
 [13] lubridate_1.7.9.2     xml2_1.3.2            codetools_0.2-16      splines_4.0.3        
 [17] polyclip_1.10-0       jsonlite_1.7.2        broom_0.7.5           ica_1.0-2            
 [21] cluster_2.1.0         dbplyr_2.0.0          png_0.1-7             uwot_0.1.10          
 [25] spatstat.sparse_2.0-0 shiny_1.6.0           sctransform_0.3.2     compiler_4.0.3       
 [29] httr_1.4.2            backports_1.2.0       assertthat_0.2.1      Matrix_1.3-4         
 [33] fastmap_1.0.1         lazyeval_0.2.2        cli_2.5.0             later_1.1.0.1        
 [37] htmltools_0.5.1.1     tools_4.0.3           igraph_1.2.6          gtable_0.3.0         
 [41] glue_1.4.2            reshape2_1.4.4        RANN_2.6.1            tinytex_0.28         
 [45] Rcpp_1.0.7            scattermore_0.7       cellranger_1.1.0      vctrs_0.3.6          
 [49] nlme_3.1-149          lmtest_0.9-38         xfun_0.19             globals_0.14.0       
 [53] rvest_0.3.6           mime_0.9              miniUI_0.1.1.1        lifecycle_1.0.0      
 [57] irlba_2.3.3           goftest_1.2-2         future_1.21.0         MASS_7.3-53          
 [61] zoo_1.8-8             scales_1.1.1          spatstat.core_2.3-0   spatstat.utils_2.2-0 
 [65] hms_0.5.3             promises_1.1.1        parallel_4.0.3        RColorBrewer_1.1-2   
 [69] gridExtra_2.3         reticulate_1.18       pbapply_1.4-3         rpart_4.1-15         
 [73] stringi_1.5.3         rlang_0.4.11          pkgconfig_2.0.3       matrixStats_0.57.0   
 [77] lattice_0.20-41       tensor_1.5            ROCR_1.0-11           patchwork_1.1.1      
 [81] htmlwidgets_1.5.3     bit_4.0.4             cowplot_1.1.0         tidyselect_1.1.0     
 [85] parallelly_1.22.0     RcppAnnoy_0.0.18      plyr_1.8.6            magrittr_2.0.1       
 [89] R6_2.5.0              generics_0.1.0        DBI_1.1.1             mgcv_1.8-33          
 [93] pillar_1.4.7          haven_2.3.1           withr_2.4.2           fitdistrplus_1.1-3   
 [97] abind_1.4-5           survival_3.2-7        future.apply_1.6.0    hdf5r_1.3.4          
[101] modelr_0.1.8          crayon_1.3.4          KernSmooth_2.23-17    spatstat.geom_2.2-2  
[105] plotly_4.9.4.1        grid_4.0.3            readxl_1.3.1          data.table_1.13.4    
[109] reprex_0.3.0          digest_0.6.27         xtable_1.8-4          httpuv_1.5.4         
[113] munsell_0.5.0         viridisLite_0.3.0   

Ref:

Loom file format specs: https://linnarssonlab.org/loompy/format/index.html

Seurat::as.loom Doc: https://www.rdocumentation.org/packages/Seurat/versions/3.1.4/topics/as.loom

Introduction to loomR: https://satijalab.org/loomr/loomr_tutorial

Interoperability between single-cell object formats: https://satijalab.org/seurat/articles/conversion_vignette.html

Conversions: h5Seurat and AnnData: https://mojaveazure.github.io/seurat-disk/articles/convert-anndata.html#converting-from-anndata-to-seurat-via-h5seurat-1

R convert matrix or data frame to sparseMatrix: https://stackoverflow.com/questions/10555210/r-convert-matrix-or-data-frame-to-sparsematrix

loom文件轉(zhuǎn)換成seurat對象: http://www.reibang.com/p/7067e0ec6ed8