作者：椰子糖
審稿：童蒙
編輯：amethyst

單細胞分析世界里數(shù)據(jù)結(jié)構(gòu)多種多樣煞聪，主流的四種數(shù)據(jù)結(jié)構(gòu)分別是Bioconductor主導(dǎo)的SingleCellExperiment，Seurat中的SeuratObject格式而账，scanpy中的AnnData格式贞让，以及大型數(shù)據(jù)存儲的loom格式涩咖。通常一種數(shù)據(jù)結(jié)構(gòu)對應(yīng)的內(nèi)容可以包含所有的分析项钮，例如seurat就可以一用到底华望，那么我們只要掌握好其中一種數(shù)據(jù)結(jié)構(gòu)就基本夠用哟玷，但也許這樣就可能會錯過其他比較好用的函數(shù)狮辽。為了更深入的了解更多好用的函數(shù)，就可以來看看各種數(shù)據(jù)結(jié)構(gòu)的轉(zhuǎn)換巢寡。

01 各種數(shù)據(jù)結(jié)構(gòu)的介紹

單細胞數(shù)據(jù)中每一個基因可以看做是單細胞的一個特征喉脖，因此單細胞數(shù)據(jù)中每一個細胞就有n個特征，在空間上可以看做是細胞具有的n個特征維度抑月。

1树叽、SingleCellExperiment

單細胞分析中的非常常用的S4對象，里面包羅萬象谦絮，那么它是如何組織的题诵？存儲了什么內(nèi)容？以下這張圖片中就已經(jīng)整體進行了展示說明挨稿。

圖中最核心的部分仇轻，是藍色的data部分；另外還有綠色的基因注釋信息feature metadata奶甘、橙色的細胞注釋信息cell metadata篷店。除了這三大件，還會包含一些下游分析結(jié)果臭家，比如PCA疲陕、tSNE降維結(jié)果就會保存在紫色的Dimension Reductions，現(xiàn)在Bioconductor上的70多個單細胞相關(guān)的R包都使用了這個對象钉赁，是單細胞分析中非常常用的一種數(shù)據(jù)結(jié)構(gòu)蹄殃。
SingleCellExperment對象又簡稱sce，就好比是一輛裝了各種集裝箱的大貨車你踩，每一個集裝箱里裝著不同的物資诅岩，每一個集裝箱都有不同的數(shù)據(jù)接口（有的接口必須是數(shù)值型矩陣結(jié)構(gòu)，有的就需要數(shù)據(jù)框結(jié)構(gòu)）带膜，又都鉚合在一起共同在貨車上前進吩谦。也有人把這些對象結(jié)構(gòu)比作是一個流水線，經(jīng)過不同的步驟膝藕，傳送帶上的對象盒子里就裝上了經(jīng)過不同函數(shù)處理的獨立又相互關(guān)聯(lián)的數(shù)據(jù)內(nèi)容式廷，中間每一步的處理有畫圖函數(shù)和指標(biāo)來進行質(zhì)控。

counts_matrix <- data.frame(cell_1 = rpois(10, 10), 
                    cell_2 = rpois(10, 10), 
                    cell_3 = rpois(10, 30))
rownames(counts_matrix) <- paste0("gene_", 1:10)

counts_matrix <- as.matrix(counts_matrix) 

有了這個矩陣逛绵，就可以用一個list構(gòu)建SingleCellExperiment對象了。

#BiocManager::install('SingleCellExperiment')
sce <- SingleCellExperiment(assays = list(counts = counts_matrix))

>sce
## class: SingleCellExperiment 
## dim: 10 3 
## metadata(0):
## assays(1): counts
## rownames(10): gene_1 gene_2 ... gene_9 gene_10
## rowData names(0):
## colnames(3): cell_1 cell_2 cell_3
## colData names(0):
## reducedDimNames(0):
## spikeNames(0):
## altExpNames(0):

要獲取這個對象中的矩陣倔韭，可以使用assay（sce,"count"）也可以使用counts(sce)，其中第一種方式也可以是其他名稱的矩陣瓢对，而第二種方式只能獲得為counts名稱的矩陣寿酌。

構(gòu)建好assays的核心后，可繼續(xù)進行拓展硕蛹，例如使用標(biāo)準(zhǔn)函數(shù)進行擴展醇疼，擴展到歸一化的矩陣，例如：

sce <- scran::computeSumFactors(sce)
sce <- scater::normalize(sce)

> sce
## class: SingleCellExperiment 
## dim: 10 3 
## metadata(1): log.exprs.offset
## assays(2): counts logcounts
## rownames(10): gene_1 gene_2 ... gene_9 gene_10
## rowData names(0):
## colnames(3): cell_1 cell_2 cell_3
## colData names(0):
## reducedDimNames(0):
## spikeNames(0):
## altExpNames(0):

這個assays就從原來存儲原始矩陣的counts法焰，增加了一個歸一化矩陣logcounts秧荆。這里對兩個矩陣進行比較發(fā)現(xiàn)normalization這一步主要是去除細胞間或者樣本文庫間的差異，是所有的細胞或者樣本具有了可比性埃仪。目前我們的sce對象中就有了2個元素乙濒。如果構(gòu)建了新的矩陣也可以增加到sce對象中去。

assays(sce)
## List of length 2
## names(2): counts logcounts

除了有assays外卵蛉，還有列的注釋信息colData颁股，列的注釋信息通常就是細胞的注釋信息，例如細胞名稱或者樣本名稱傻丝、批次甘有、作者等等，對應(yīng)的上圖中的橙色部分葡缰。首先我們來設(shè)置一個細胞批次的注釋信息：

cell_metadata <- data.frame(batch = c(1, 1, 2))
rownames(cell_metadata) <- paste0("cell_", 1:3)

接著添加到sce對象中亏掀，通過直接構(gòu)建的方式添加：

sce <- SingleCellExperiment(assays = list(counts = counts_matrix),
                           colData = cell_metadata)

也可以后續(xù)添加：

colData(sce) <- DataFrame(cell_metadata)

加入了sce對象后，可以通過colData(sce)來進行獲取泛释，或者sce$batch查看滤愕。

細胞注釋的信息除了手動添加外，還可以通過scater包的calculateQCMetrics()進行計算胁澳，得到幾十項細胞的質(zhì)量信息该互，并添加進去。

sce <- scater::addPerCellQC(sce)
colData(sce)[, 1:5]
## DataFrame with 3 rows and 5 columns
##            batch       sum  detected percent_top_50 percent_top_100
##        <numeric> <integer> <integer>      <numeric>       <numeric>
## cell_1         1        80        10            100             100
## cell_2         1        88        10            100             100
## cell_3         2       309        10            100             100

rowData是行的注釋信息韭畸，同細胞或者樣本一樣宇智，基因也有自己的注釋信息蔓搞，它是一個數(shù)據(jù)框的結(jié)構(gòu)。

 # 一開始rowData(sce)是空的随橘，可以添加
 sce <- scater::addPerFeatureQC(sce)
 rowData(sce)
 
 ## DataFrame with 10 rows and 2 columns
 ##              mean  detected
 ##         <numeric> <numeric>
 ## gene_1    16.0000       100
 ## gene_2    14.3333       100
## gene_3    16.0000       100
## gene_4    18.6667       100
## gene_5    15.3333       100
## gene_6    16.6667       100
## gene_7    17.6667       100
## gene_8    13.0000       100
## gene_9    16.6667       100
## gene_10   14.6667       100

除了以上這些基本的喂分，由于基因的注釋信息是數(shù)據(jù)框的結(jié)構(gòu)，因此可以按照位置或者名字來取子集机蔗，sce[c("gene_1","gene_2")]蒲祈。另外還有reduceDims，用來存儲原始矩陣的降維結(jié)構(gòu)萝嘁，可以通過PCA梆掸，tSNE和Umap獲得。

SingleCellExperiment對象的兼容性強牙言，可以用在多種單細胞分析的R包中作為輸入或者中間文件酸钦，方便了對數(shù)據(jù)的傳輸和協(xié)作。

2咱枉、SeuratObject

單細胞數(shù)據(jù)中的另一常用對象是seurat對象卑硫，主要用在seurat包、monocle包等單細胞分析的R包中蚕断。其結(jié)構(gòu)與SingleCellExperiment對象非常相似欢伏，在R中，這兩個對象也能夠通過as函數(shù)進行轉(zhuǎn)換亿乳。seurat對象結(jié)構(gòu)包含以下這些slots硝拧。

有人把seurat對象比作是流水線上的具有不同盒子的容器，經(jīng)過不同的工序葛假，會在不同的盒子里（slots）增加內(nèi)容河爹，而這部分工序就叫convertor,包括NormallizeData()、FindVariableFeeature()等桐款，而在每一步工序中對數(shù)據(jù)進行監(jiān)管咸这，控制的工具就叫Inspector,包括VlnPlot()，pbmc[['RNA']]@data等魔眨。

Assays：assays里有一個元素“RNA”侄刽，訪問assays對象的內(nèi)部結(jié)構(gòu)，例如pbmc[['RNA']]朋凉。

RNA：是assay眨补，其可以包含多個matrix：

• Counts：原始的表達量count矩陣混稽；
• Data：原始數(shù)據(jù)經(jīng)過normalized的數(shù)據(jù)珍昨；
• Scale.data：數(shù)據(jù)經(jīng)過scaling后，存放位置吓揪；
• Key：每個assay對象都有一個key，例如‘rna_’所计；
• Var.features：普通向量柠辞，高表達變異的基因。

Meta.data：是細胞的注釋信息的數(shù)據(jù)框主胧，行是細胞叭首，列是細胞的屬性。

Active.assay：當(dāng)前激活的assay對象踪栋。

標(biāo)準(zhǔn)的seurat流程可以這樣進行焙格，包括創(chuàng)建對象，標(biāo)準(zhǔn)化夷都，尋找高突變特征间螟，歸一化，聚類等：

pbmc.counts <- Read10X(data.dir ="~/Downloads/pbmc3k/filtered_gene_bc_matrices/hg19/")
pbmc <- CreateSeuratObject(counts = pbmc.counts)
pbmc <- NormalizeData(object = pbmc)
pbmc <- FindVariableFeatures(object = pbmc)
pbmc <- ScaleData(object = pbmc)
pbmc <- RunPCA(object = pbmc)
pbmc <- FindNeighbors(object = pbmc)
pbmc <- FindClusters(object = pbmc)
pbmc <- RunTSNE(object = pbmc)
DimPlot(object = pbmc, reduction = "tsne")

seurat對象操作如下：

# Get cell and feature names, and total numbers
colnames(x = pbmc) ##細胞
Cells(object = pbmc) ##同上獲取細胞名稱
rownames(x = pbmc) ##獲取基因名稱
ncol(x = pbmc) ##獲取細胞數(shù)量
nrow(x = pbmc)  ##獲取基因數(shù)目

可以使用一些函數(shù)對行列等進行獲取损肛，也可以用Idents獲取細胞的分類情況。

# Get cell identity classes
Idents(object = pbmc)
levels(x = pbmc)
table(Idents(pbmc))  ##獲取每個細胞類型數(shù)目的表格

# Stash cell identity classes
pbmc[["old.ident"]] <- Idents(object = pbmc)
pbmc <- StashIdent(object = pbmc, save.name = "old.ident")
##設(shè)置細胞類別
# Set identity classes
Idents(object = pbmc) <- "CD4 T cells"
Idents(object = pbmc, cells = 1:10) <- "CD4 T cells"

# Rename identity classes
pbmc <- RenameIdents(object = pbmc, `CD4 T cells` = "T Helper cells")

另外也可以通過subset來取對象子集荣瑟，或者通過merge來合并多個seurat對象治拿。

由于seurat包的好用易用，所以了解和學(xué)習(xí)seurat對象結(jié)構(gòu)也是很有必要的笆焰。

3劫谅、AnnData

隨著單細胞技術(shù)的發(fā)展，單細胞測序的技術(shù)門檻也越來越低嚷掠，單細胞的數(shù)據(jù)呈幾何倍數(shù)的增加捏检，要分析大百萬級細胞數(shù)量的單細胞樣本，給依賴R分析的軟件帶來了一定的困難（由于R的存儲方式的限制）不皆。于是就有了python版的對象AnnData：

可以通過 sc.pl.highest_expr_genes(adata, n_top=20)獲取表達量top20的基因贯城。

sc.pp.filter_cells(adata, min_genes=200) # 去除表達基因200以下的細胞 sc.pp.filter_genes(adata, min_cells=3) # 去除在3個細胞以下表達的基因

在anndata對象中對線粒體基因進行過濾，以及標(biāo)準(zhǔn)化數(shù)據(jù)和后續(xù)的差異分析及聚類的代碼如下：

mito_genes = adata.var_names.str.startswith('MT-') 
adata.obs['percent_mito'] = np.sum( adata[:, mito_genes].X, axis=1).A1 / np.sum(adata.X, axis=1).A1 
adata.obs['n_counts'] = adata.X.sum(axis=1).A1
sc.pl.scatter(adata, x='n_counts', y='percent_mito') 
sc.pl.scatter(adata, x='n_counts', y='n_genes')

adata = adata[adata.obs['n_genes'] < 4000, :] 
adata = adata[adata.obs['percent_mito'] < 0.3, :]

#數(shù)據(jù)標(biāo)準(zhǔn)化
sc.pp.normalize_per_cell(adata, counts_per_cell_after=1e4)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
sc.pl.highly_variable_genes(adata)

##保留差異基因進行后續(xù)分析
adata = adata[:, adata.var['highly_variable']]
sc.pp.regress_out(adata, ['n_counts', 'percent_mito'])
sc.pp.scale(adata, max_value=10)
sc.tl.pca(adata, svd_solver='arpack') # PCA分析 sc.pl.pca(adata, color='CST3') #繪圖

##碎石圖霹娄，選擇多少PCA進行后續(xù)分析能犯，用于計算細胞間的相鄰距離
sc.pl.pca_variance_ratio(adata, log=True)
adata.write("pca_results.h5ad")
聚類分析
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40) ##neighbor個數(shù)越多，聚類越少
sc.tl.umap(adata) 
sc.pl.umap(adata, color=['CST3', 'NKG7', 'PPBP'])
sc.tl.louvain(adata)
sc.pl.umap(adata, color=['louvain'])

由上也可以看出anndata對象具有切片特性犬耻，熟悉pandas的人踩晶，操作起來就很熟練。

4枕磁、loom

Loom是非常大的omics數(shù)據(jù)集的有效文件格式渡蜻，由主矩陣、可選的附加層、可變數(shù)量的行和列注釋以及稀疏的圖形對象組成茸苇。我們使用loom文件存儲單細胞基因表達數(shù)據(jù)：主矩陣包含實際的表達值（每個細胞一列排苍，每個基因一行）；行注釋和列注釋包含基因和細胞的元數(shù)據(jù)税弃，例如名稱纪岁、染色體、位置（對于基因）和品系则果、性別幔翰、年齡（對于細胞）。圖對象用于存儲用于基于圖的聚類的最近鄰圖西壮。

在這里附上loom處理的相關(guān)代碼遗增，以供大家參考：

##創(chuàng)建loom對象
library(Seurat) 
library(loomR) 
library(dplyr) 
pbmc_small_loom<-create(filename = "pbmc.small.loom",data = pbmc_small@assays$RNA@counts,overwrite = T)
 #這里也可以直接as.loom()
##查看loom文件中的信息
pbmc_small_loom pbmc_small_loom$matrix[1:6,1:6] pbmc_small_loom$col.attrs$CellID[1:6] pbmc_small_loom$row.attrs$Gene[1:6]
##提取loom中的信息：
pbmc_small_loom$get.attribute.df(MARGIN = 1,attributes = "Gene")[1:6,] pbmc_small_loom$get.attribute.df(MARGIN = 2,attributes = "CellID")[1:6,]
##添加信息
# Generate random ENSEMBL IDs for demonstration purposes 
ensembl.ids <- paste0("ENSG0000", 1:length(x = pbmc_small_loom$row.attrs$Gene[])) 
pbmc_small_loom$add.row.attribute(list(ensembl.id = ensembl.ids), overwrite = TRUE) pbmc_small_loom$get.attribute.df(MARGIN = 1)[1:6,]
##進行seurat操作
pbmc_small_seurat<-as.Seurat(pbmc_small_loom) pbmc_small_seurat<-NormalizeData(pbmc_small_seurat)%>%ScaleData()
##關(guān)閉loom
pbmc_small_loom$close_all()

操作loom結(jié)構(gòu)時，一定注意關(guān)閉loom文件款青，否則再次讀取的時候做修，就會有異常。

02 數(shù)據(jù)結(jié)構(gòu)中結(jié)構(gòu)轉(zhuǎn)換的軟件

雖然各軟件都有對應(yīng)的整套的分析流程抡草，但是也可能碰巧有你想用的軟件或者包饰及，不是你熟悉的數(shù)據(jù)格式，那么我們就可以通過以下一些方法進行轉(zhuǎn)換康震。

方法一：seurat中的函數(shù)直接轉(zhuǎn)換

導(dǎo)入對應(yīng)的包：

1# install scater https://bioconductor.org/packages/release/bioc/html/scater.html
2library(scater)
3library(Seurat)
4# install SeuratDisk from GitHub using the remotes package remotes::install_github(repo =
5# 'mojaveazure/seurat-disk', ref = 'develop')
6library(SeuratDisk)
7library(patchwork)

進行轉(zhuǎn)換：
converting to/from SingleCellExperiment：

1pbmc <- readRDS(file = "../data/pbmc3k_final.rds")
2pbmc.sce <- as.SingleCellExperiment(pbmc)
3p1 <- plotExpression(pbmc.sce, features = "MS4A1", x = "ident") + theme(axis.text.x = element_text(angle = 45, 
4    hjust = 1))
5p2 <- plotPCA(pbmc.sce, colour_by = "ident")
6p1 + p2
7
 1# download from hemberg lab
 2# https://scrnaseq-public-datasets.s3.amazonaws.com/scater-objects/manno_human.rds
 3manno <- readRDS(file = "../data/manno_human.rds")
 4manno <- runPCA(manno)
 5manno.seurat <- as.Seurat(manno, counts = "counts", data = "logcounts")
 6# gives the same results; but omits defaults provided in the last line
 7manno.seurat <- as.Seurat(manno)
 8Idents(manno.seurat) <- "cell_type1"
 9p1 <- DimPlot(manno.seurat, reduction = "PCA", group.by = "Source") + NoLegend()
10p2 <- RidgePlot(manno.seurat, features = "ACTB", group.by = "Source")
11p1 + p2

Converting to/from loom:

1# download from linnarsson lab
2# https://storage.googleapis.com/linnarsson-lab-loom/l6_r1_immune_cells.loom
3l6.immune <- Connect(filename = "../data/l6_r1_immune_cells.loom", mode = "r")
4l6.seurat <- as.Seurat(l6.immune)
5Idents(l6.seurat) <- "ClusterName"
6VlnPlot(l6.seurat, features = c("Sparc", "Ftl1", "Junb", "Ccl4"), ncol = 2)

方法二：sceasy進行轉(zhuǎn)換

1##安裝
2BiocManager::install(c("LoomExperiment", "SingleCellExperiment"), update=F)
3install.packages('reticulate')
4devtools::install_github("cellgeni/sceasy")
5
6##python安裝anndata包(版本 < 0.6.20), oompy包(版本 < 3.0.0)使用library(sceasy)
7library(reticulate)
8use_condaenv('EnvironmentName')
9loompy <- reticulate::import('loompy')

安裝好軟件后燎含，就可以進行隨意的轉(zhuǎn)換啦。

 1###Seurat to AnnData
 2sceasy::convertFormat(seurat_object, from="seurat", to="anndata", outFile='filename.h5ad')
 3
 4###AnnData to Seurat
 5sceasy::convertFormat(h5ad_file, from="anndata", to="seurat", outFile='filename.rds')
 6
 7###Seurat to SingleCellExperiment
 8sceasy::convertFormat(seurat_object, from="seurat", to="sce", outFile='filename.rds')
 9
10###SingleCellExperiment to AnnData
11sceasy::convertFormat(sce_object, from="sce", to="anndata", outFile='filename.h5ad')
12
13###SingleCellExperiment to Loom
14sceasy::convertFormat(sce_object, from="sce", to="loom", outFile='filename.loom')
15
16###Loom to AnnData
17sceasy::convertFormat('filename.loom', from="loom", to="anndata", outFile='filename.h5ad')
18
19###Loom to SingleCellExperiment
20sceasy::convertFormat('filename.loom', from="loom", to="sce", outFile='filename.rds')

有了這些工具腿短，就可以讓我們在單細胞數(shù)據(jù)的海洋的愜意遨游屏箍。

參考資料：

1、https://osca.bioconductor.org/data-infrastructure.html
2橘忱、https://bioc.ism.ac.jp/packages/3.6/bioc/vignettes/SingleCellExperiment/inst/doc/intro.html
3赴魁、https://satijalab.org/seurat/articles/essential_commands.html
4、http://linnarssonlab.org/loompy/format/index.html
5钝诚、https://satijalab.org/seurat/mca_loom.html
6颖御、https://satijalab.org/seurat/articles/conversion_vignette.html
7、https://github.com/cellgeni/sceasy
8凝颇、https://www.bilibili.com/read/cv8675277