2.1加载数据

我们使用TENxBrainData函数来下载相关文件Bioconductor的ExperimentHub网络资源。这包括包含计数的HDF5文件，以及行(基因)和列(单元格)上的元数据。输出为单个SingleCellExperiment对象的SingleCellExperiment包中。它等价于aSummarizedExperiment类，但具有许多特定于单单元数据的特性。

library(TENxBrainData) tenx <- TENxBrainData() tenx

##类:singlecel实验## dim: 27998 1306127 ##元数据(0):## assays(1):计数## rownames: NULL ## rowData names(2):集合符号## colnames(1306127): AAACCTGAGATAGGAG-1 AAACCTGAGCGGCTTC-1…## TTTGTCAGTTAAAGTG-133 tttgtcatctgaaga -133 ## colData names(4):条形码序列库鼠标## reducedDimNames(0): ## mainExpName: NULL ## altExpNames(0):

第一次调用TENxBrainData ()由于需要下载一些中等大小的文件，需要一些时间。然后将这些文件存储在本地，以便相同或新会话中的后续调用能够快速进行。

计数矩阵本身表示为aDelayedMatrix从DelayedArray包中。这将底层HDF5文件包装在一个可以用r进行操作的容器中，每个计数表示分配给特定细胞中特定基因的唯一分子标识符(UMIs)的数量。

计数(tenx)

## <27998 x 1306127>矩阵类DelayedMatrix和类型“integer”:## aaacctgagatagag -1…Tttgtcatctgaaaga-133 ##[1，] 0。0 ##[2，] 0。0 ##[3，] 0。0 ##[4，] 0。0 ##[5，] 0。0 ## ... ...##[27994，] 0。0 ## [27995，]0 ##[27996，] 0。 0 ## [27997,] 0 . 0 ## [27998,] 0 . 0

2.2探索数据

为了快速浏览数据集，我们在计数矩阵上计算一些汇总统计信息。我们增加DelayedArray块大小，表明我们可以使用最多2 GB的内存从磁盘加载数据到内存。

选项(DelayedArray.block.size = 2 e9)

我们对图书馆的大小感兴趣colSums(计数(tenx))，每个细胞表达的基因数量colsum (counts(tenx) != 0)，以及单元格的rowMeans(counts(tenx))的平均表达。天真的实现可能是

自由。n.exprs <- colsum (counts(tenx) != 0L) ave.exprs <- rowMeans(counts(tenx))

但是，数据是从磁盘读取的，磁盘访问相对较慢，并且幼稚的实现读取数据三次。相反，我们将把数据分成大约10,000个单元格的列“块”;我们在数据的子集上执行此操作，以减少探索阶段的计算时间。

tenx20k <- tenx[， seq_len(20000)] chunksize <- 5000 cidx <- snow::splitIndices(ncol(tenx20k)， ncol(tenx20k) / chunksize)

并遍历文件，读取数据并在每次迭代中积累统计信息。

自由。size <- n.exprs <- numeric(ncol(tenx20k))exprs <- numeric(nrow(tenx20k)) for (i in head(cidx, 2)) {message("."， appendLF=FALSE) m <- as.matrix(counts(tenx20k)[，i, drop=FALSE])size [i] <- colsum (m) n.exprs[i] <- colsum (m != 0) tot. exprs[i] <- colsum (m != 0)Exprs <- tot。exprs + rowsum (m)} ave.exprs <- tot. exprs。Exprs / ncol(tenx20k)

由于计算的开销很大，并且在将来可能很有用，所以我们注释tenx20k对象

colData (tenx20k)美元的自由。大小<- lib。大小colData (tenx20k) $ n。exprs <- n.exprs rowData(tenx20k)$ave. exprs rowData(tenx20k)Exprs <- ave.exprs

库的大小遵循近似对数正态分布，并且惊人地小。

hist(log10(colData(tenx20k)$lib.sizes)， xlab=expression(Log[10] ~ "库大小")，col="grey80")

在每个样本中只检测到几千个基因的表达。

嘘(colData (tenx20k) $ n。exprs, xlab="检测到的基因数量"，col="grey80")

平均表达式值(读计数)很小。

hist(log10(rowData(tenx20k)$ave.exprs)， xlab=expression(Log[10] ~ "Average count")， col="grey80")

我们还检查了这个数据集中最高表达的基因。

o <- order(rowData(tenx20k)$ave。exprs，递减=TRUE)头(rowData(tenx20k)[o，])

## 6行3列的数据帧##集成符号ave.exprs ## <字符> <数组> <数字> ## # 1 ENSMUSG00000092341 Malat1 49.1875 ## 2 ENSMUSG00000049775 Tmsb4x 24.3624 ## 3 ENSMUSG00000072235 Tuba1a 23.2797 ## 4 ENSMUSG00000064357 mt-Atp6 17.1341 ## 5 ENSMUSG00000064358 mt-Co3 14.3276 ## 6 ENSMUSG00000028832 Stmn1 13.9624

更先进的分析程序在各种Bioconductor软件包-请参阅SingleCellbiocViews获取更多详细信息。

2.3节约计算

保存tenx对象以标准的方式，例如，

destination <- tempfile() savverds (tenx, file = destination)

保存行数据、列数据和元数据为R对象，并记住从该对象派生的HDF5文件的位置和子集。对象可以读入newR会话readRDS(目的地)，只要HDF5文件保留在原始位置。

3．1并行计算

行和列的汇总统计信息可以并行计算，例如使用bpiterate ()在BiocParallel包中。我们加载包并启动5个“snow”工人(单独的进程)。

库(BiocParallel)注册(bpstart (SnowParam (5)))

此函数需要一个迭代器生成数据块。迭代器返回一个函数，该函数本身返回每个块的开始列和结束列索引，直到不再有块为止。

迭代器<- function(tenx, cols_per_chunk = 5000, n = Inf) {start <- seq(1, ncol(tenx)， by = cols_per_chunk) end <- c(tail(start， -1) - 1L, ncol(tenx)) n <- min(n, length(start)) i <- 0L function() {if (i == n) return(NULL) i <<- i + 1L c(start[i]， end[i])}}

下面是正在运行的迭代器

Iter <- iterator(tenx) Iter ()

## [1] 1 5000

iter ()

## [1] 5001 10000

iter ()

## [1] 10001 15000

bpiterate ()需要一个作用于每个数据块的函数。它接收迭代器的输出，以及它可能需要的任何其他参数，并返回该块的汇总统计信息

fun <- function(crng, counts，…){## ' fun() '需要自包含一些并行后端suppressPackageStartupMessages({library(TENxBrainData)}) m <- as。matrix(counts[， seq(crng[1]， crng[2])]) list(row = list(n = rownumbers (m != 0)， sum = rownumbers (m)， sumsq = rownumbers (m * m))， column = list(n = colnumbers (m != 0)， sum = colnumbers (m)， sumsq = colnumbers (m * m)))}

我们可以用

Res <- fun(iter()， unname(counts(tenx))) str(Res)

##列表2 ## $行:列表3 ## ..$ n: num[1:27998] 114 0 0 4 0…# # . .$ sum: num[1:27998] 120 0 0 4 0…# # . .$ sumsq: num[1:27998] 134 0 0 4 0…## $列:列表3 ## ..$ n: num[1:5000] 2077 2053 2503 1617 1402…# # . .$ sum: num[1:5000] 4740 4309 6652 2930 2408…# # . .$ sumsq: num [1:5000] 52772 32577 90380 19598 16622 ...

最后,bpiterate ()需要一个函数来减少返回的连续值有趣的()

reduce <- function(x, y) {list(row = Map(' + '， x$row, y$row)， column = Map(' c '， x$column, y$column))}

测试是

Str (reduce(res, res))

##列表2 ## $行:列表3 ## ..$ n: num[1:27998] 228 0 0 8 0…# # . .$ sum: num[1:27998] 240 0 0 8 0…# # . .$ sumsq: num[1:27998] 268 0 0 8 0…## $列:列表3 ## ..$ n: num[1:10000] 2077 2053 2503 1617 1402…# # . .$ sum: num[1:10000] 4740 4309 6652 2930 2408…# # . .$ sumsq: num [1:10000] 52772 32577 90380 19598 16622 ...

把这些片段放在一起，计算前25000列，我们有

res <- bpiterator (iterator(tenx, n = 5)， fun, counts = unname(counts(tenx))， REDUCE = REDUCE, REDUCE .in.order = TRUE) str(res)

##列表2 ## $行:列表3 ## ..$ n: num[1:27998] 579 1 0 29 0…# # . .$ sum: num[1:27998] 602 1 0 29 0…# # . .$ sumsq: num[1:27998] 652 1 0 29 0…## $列:列表3 ## ..$ n: num[1:25000] 1807 1249 2206 1655 3326…# # . .$ sum: num[1:25000] 4046 2087 4654 3193 8444…# # . .$ sumsq: num [1:25000] 35338 14913 31136 21619 106780 ...

3.2使用rle压缩的HDF5数据

10x Genomics的数据也以压缩格式发布，可从ExperimentHub获得

library(ExperimentHub) hub <- ExperimentHub()查询(hub， "TENxBrainData")

## ExperimentHub与8条记录## # snapshotDate(): 2022-04-19 ## $dataprovider: 10X Genomics ## # $species: Mus musculus ## # $rdataclass: character ## #附加mcols():例如，'object[["EH1039"]]]' ## ## title ## EH1039 | Brain scRNA-seq数据，'基于hdf5的10X基因组学'格式## EH1040 | Brain scRNA-seq数据，'密集矩阵'格式## EH1041 | Brain scRNA-seq数据，样本(列)注释## EH1042 | Brain scRNA-seq数据，基因(行)注释## EH1689 | Brain scRNA-seq数据20k子集，“基于hdf5的10x基因组学”格式## EH1690 |脑scRNA-seq数据20k子集，“密集矩阵”格式## EH1691 |脑scRNA-seq数据20k子集，样本(列)注释## EH1692 |脑scRNA-seq数据20k子集，基因(行)注释

fname <- hub[["EH1039"]]

文件的结构可以使用h5ls ()命令从rhdf5．

h5ls(帧)

## group name otype dclass dim ## 0 /mm10 H5I_GROUP ## 1 /mm10 barcode H5I_DATASET STRING 1306127 ## 2 /mm10 data H5I_DATASET INTEGER 2624828308 ## 3 /mm10 gene_names H5I_DATASET STRING 27998 ## 4 /mm10 genes H5I_DATASET STRING 27998 ## 5 /mm10 indices H5I_DATASET INTEGER 2624828308 ## 6 /mm10 indptr H5I_DATASET INTEGER 1306128 ## 7 /mm10 shape H5I_DATASET INTEGER 2

非零计数在/ mm10 /数据路径。/ mm10 /指标表示对应于每个非零计数的行索引。/ mm10 indptr将数据和索引分成连续的列。例如

Start <- h5read(fname， "/mm10/indptr"， Start =1, count=25001) head(Start)

## [1] 0 1807 3056 5262 6917 10243

检索偏移量到/ mm10 /数据数据的前25001列。偏移量是基于0的，因为HDF5使用基于0的索引;我们有时需要添加1以方便在中使用R．

这里我们将前25000列数据读入R,使用data.table为了有效地计算这些大数据。

Library (data.table) dt <- data.table。Table (row = h5read(fname， "/mm10/ indexes "， start = 1, count = tail(start, 1)) + 1, column = rep(seq_len(length(start) - 1)， diff(start)))， count = h5read(fname， "/mm10/data"， start = 1, count = tail(start, 1))) dt

##行列计数## 1:27995 1 1 ## 2:27921 1 19 ## 3:27918 1 14 ## 4:27915 1 40 ## 5:27914 1 29 ##——## 51028822:63 25000 9 ## 51028823:39 25000 2 ## 51028824:38 25000 1 ## 51028825:17 25000 2 ## 51028826:8 25000 1

然后是行和列摘要

dt[， list(n = .N, sum = sum(count)， sumsq = sum(count * count))， keyby=row]

##行n求和求和## 1:1 579 602 652 ## 2:21 1 1 ## 3:4 29 29 29 ## 4:6 215 615 3501 ## 5:7 5 5 5 ##——## 20191:27980 77 77 ## 20192:27994 22 22 22 ## 20193:27995 14773 28213 77081 ## 20194:27996 1946 2085 2393 ## 20195:27998 16 16 16

dt[， list(n = .N, sum = sum(count)， sumsq = sum(count * count))， keyby=column]

##列n sum sumsq ## 1: 1 1807 4046 35338 ## 2: 2 1249 2087 14913 ## 3: 3 2206 4654 31136 ## 4: 4 1655 3193 21619 ## 5: 5 3326 8444 106780 ##——# 24996:24996 1142 2091 16827 ## 24997:24997 2610 5772 131212 ## 24998:24998 2209 4492 36600 ## 24999:24999 1049 1791 15795 ## 25000:25000 2015 4043 28809

遍历25000列密集数据大约需要3分钟的计算时间(使用6个核大约需要30秒)，而对于稀疏数据只需要几秒钟。处理整个稀疏数据集仍然需要块处理，除非在大内存机器上，并且将受益于并行计算。在后一种情况下，每个块处理少于25000列将减少每个块的内存消耗，从而允许更多的处理核心运行，从而提高整体处理速度。