CyTOF分析流程（一）：基于R语言的CyTOF数据预处理

小鱼爸

2023-08-03

简介

质谱流式(Mass cytometry，CyTOF)得到越来越多的应用，适合用于深入表征健康和疾病组织异质性。目前，原始文件的预处理通常是在独立软件工具中完成的，这使得重复再现比较困难。在这里，我们介绍一个基于CATALYST的R处理流程，它以完全可重复的方式涵盖了质谱流式分析下游所需的所有预处理步骤。该流程包括文件串联file concatenation、磁珠归一化bead-based normalization、单细胞反卷积single-cell deconvolution、溢出补偿spillover compensation以及碎片和双峰去除后的活细胞门控。

数据集介绍

本分析流程中使用的数据是在Tumor Profiler project中产生的multi-center observational study研究，研究不同创新技术的相关性，包括CyTOF、成像质谱流式细胞术、单细胞DNA和RNA测序以及离体药物测试，以提高晚期癌症患者的诊断。 分析的样本包括2020年在苏黎世大学医院收集的肿瘤活检和血液样本。这些样本以一系列正常背景（包括商业可用的细胞系、健康捐献者的PBMCs和PHA激活的PBMCs）作为参考对照，进行CyTOF评估。 本研究使用的数据集来自一次CyTOF实验，也称为批处理batch，其中9个参考样本、2个血液样本和2个肿瘤样本使用20孔条形码板进行barcode编码。选择这种设计是为了在基于分位数缩放的独立实验中实现质量控制和批量校正。用40-Ab panel对混合细胞进行染色，以对样品的免疫系统进行深入的表征。 在整个分析流程中，我们将利用一组元数据metadata进行标准预处理步骤（标准化、去barcode和补偿），以及之前在七个独立实验中获得的各种质量控制。原始数据下载地址：https://tpreports.nexus.ethz.ch/download/catalyst/。
表1给出了所使用的元数据的概述。

表1：整个流程中使用的元数据文件的概述，包括每个文件的描述、以及预处理或质量控制的目的。

Description	Purpose
normalization_beads.fcs
Beads identified using CATALYST during the normalization step of a previous CyTOF experiment.	Used as reference beads to correct for changes in signal sensitivity over time.
ref_bead_counts.csv
A table of mean dual counts for the 6 different bead channels (columns) obtained from 7 previous CyTOF experiments (rows).	Used as a reference to assess the measurement sensitivity.
debarcoding_scheme.csv
A binary barcoding scheme of 6-choose-3 = 20 barcodes with columns corresponding to barcode channel masses (101, 104, 105, 106, 108, 110) and rows corresponding to barcodes (7 empty, 9 references, 2 PBMC and 2 tumor samples)	Used for single-cell deconvolution of multiplexed samples.
spillover_matrix.csv
A spillover matrix calculated with CATALYST from beads single-stained with each of the 40 antibodies included in the panel used in this study. The matrix contains, for each measurement channel (rows), the percentage of signal received by all other channels (columns).	Used for correction of spillover.
ref_cell_counts.csv
A table of the number of cells measured in 7 previous experiments, each including 4 cell line, 3 PBMC and 2 tumor references samples (63 samples in total).	Used to assess reference sample cell yields in the current in comparison to previous experiments
sample_cell_counts.csv
A table of the number of cells measured in 7 previous experiments, each including 2 PBMC and 2 tumor samples each (28 samples in total).	Used to assess sample cell yields in the current in comparison to previous experiments
ref_marker_levels.csv
A table of the 98th expression percentiles for each target (columns) across 7 previous experiments (rows).	Used to assess the staining efficiency of the current experiment

整个分析流程中使用和返回的大多数数据都保存在SingleCellExperiment包中的SingleCellExperiment(SCE)类对象中。这种数据结构可以存储所有单细胞相关数据（测量数据及其转换；细胞、features和实验范围的元数据；降维结果），从而减少容易出错的数据操作。 SCE的关键组成部分是类矩阵的assay，其中行是特征features（目标），列是观察值（cells），用于存储测量数据及其派生的任何数据。与细胞相关的元数据存储在colData下，特征元数据存储在rowData下，任何实验范围的元数据都可以存储在元数据槽中。最后，SCE可以在reduceDims下存储任意数量的降维。有关SCE的用法和结构的更详细描述，可以参考SingleCellExperiment包的文档。 在本分析示例中，原始测量数据（FCS文件）以及所有元数据（用于去barcode、标准化、补偿和质量控制）需要位于 data/子目录内（相对于代码运行的位置）；否则，所提供的文件路径需要修改。

结果

这里介绍的分析流程描述了将原始质谱流式数据处理到用户可以继续进行下游分析（例如，降维、差异分析、轨迹推断）的状态所需的所有步骤。该过程包括串联各个采集文件、排除信号不稳定的部分采集、通过微珠归一化校正时间依赖性信号漂移、通过补偿校正信号溢出、通过自动门控选择感兴趣的细胞 ，以及批次效应的校正。 该工作流程以来自单个CyTOF实验的数据为例，该实验通过连续3次采集（单个 FCS文件）收集与校准珠混合的15个barcoded样本。 我们使用CATALYST执行关键的预处理步骤，包括：串联、标准化、去barcode标记和补偿； 用于门控的openCyto和flowWorkspace；ggplot2、ggcyto和patchwork用于可视化；flowCore、reshape2和dplyr用于数据访问和操作；和mvtnorm来计算多边门控。 因此，我们的工作流程仅限于以下依赖项：

library(CATALYST)
library(dplyr)
library(flowCore)
library(flowWorkspace)
library(ggcyto)
library(ggplot2)
library(mvtnorm)
library(openCyto)
library(patchwork)
library(reshape2)
library(pander)

除了标准预处理步骤外，我们还包括质量控制(QC)步骤来评估CyTOF灵敏度、染色功效和细胞产量； 这些依赖于先前实验（n = 7）的结果作为参考。 为了实现一致的可视化，我们为箱线图定义了一个通用的绘图主题，用于将当前实验与之前的实验进行比较：

qc_theme <- list(
 theme_bw(base_size = 8), theme(
  panel.grid.minor = element_blank(),
  panel.grid.major.x = element_blank(),
  plot.title = element_text(face = "bold"),
  axis.text = element_text(color = "black"),
  axis.text.x = element_text(angle = 45, hjust = 1)))

构建一个SingleCellExperiment对象

默认情况下，flowCore的read.FCS()函数是read.flowSet()的基础，用于读取FCS文件，转换通道强度并删除离群值。为了避免这些情况，我们建议读取时参数设置为transformation = FALSE和truncate_max_range = FALSE；默认情况下，CATALYST 的prepData()函数将使用这些设置读入文件。

如上所述，SCE类允许在单个对象中保留多个数据转换transformations。因此，当应用转换以获得类似表达的数据时，我们可以将转换后的数据存储在单独的测定中，而不会覆盖原始计数数据。通过这种方式，在整个预处理过程中生成和使用的任何数据（例如，标准化、补偿或批量校正计数及其 arcsinh 转换的对应项）原则上都可以保留，并写入中间FCS文件，以便在R之外进行备份或质量控制。但值得注意的是，此过程可能会导致大型数据集内存不足，在这种情况下，建议在每个步骤中覆盖数据； 如果没有另外指定，CATALYST默认会覆盖。

可以使用CATALYST的prepData()函数构建SCE，该函数接受包含一个或多个FCS文件的目录路径、FCS 文件名的字符向量、单个flowFrame或flowFrame列表或flowSet。默认情况下（transform = TRUE），将cofactor = 5的arcsinh变换应用于输入（count）数据，并将生成的表达矩阵存储在SCE的exprs中：

# construct ’SingleCellExperiment’
fcs <- list.files("data", "acquisition", full.names = TRUE)
(sce <- prepData(fcs, transform = TRUE, cofactor = 5))

## class: SingleCellExperiment 
## dim: 63 368152 
## metadata(2): experiment_info chs_by_fcs
## assays(2): counts exprs
## rownames(63): 75As CD15 ... 208Pb CD45
## rowData names(4): channel_name marker_name marker_class use_channel
## colnames: NULL
## colData names(1): sample_id
## reducedDimNames(0):
## mainExpName: NULL
## altExpNames(0):

SCE有两种检测，包含双离子计数（counts）和cofactor-5 arcsinh转换计数（exprs）。 用于转换的辅因子存储在对象的内部元数据(int_metadata(sce)$cofactor)中，以及原始FCS文件中每个细胞的来源存储在元数据列sample_id下（可通过colData(sce)$sample_id或sce$sample_id访问）。在我们的数据集中，FCS文件对应于采集数据而不是生物样本。因此，我们将元数据变量sample_id重命名为file_id以避免歧义：

i  <- match("sample_id", names(colData(sce)))
names(colData(sce))[i]  <- "file_id"

所有采集的细胞总数对应于SCE中的列数(ncol(sce)：368152)。我们可以通过列出file_ids来总结每个文件中的细胞数量：

data.frame(
   file_id = levels(sce$file_id),
   n_cells = tabulate(sce$file_id))

##   file_id n_cells
## 1      V1   48675
## 2      V2  125607
## 3      V3  193870

在质谱和流式细胞术中，每个feature都有与其对应的通道和目标。从上面输出sce变量可以看出，prepData()默认使用目标作为行名。 我们可以使用channels()访问器检索每个特征feature的测量通道，并使用通道金属mentals和质量masses来提取与不同预处理步骤相关的特征索引。也就是说，我们将测量DNA的通道分配给变量dna（在这里是Ir191和Ir193），并将用于活细胞门控的通道（此处为DNA的Ir191和顺铂的Pt194）分配给live：

# store character vector or channels names
chs <- channels(sce)
# store DNA & live channel indices
dna <- grep("^Ir", chs)
live <- grep("191|194", chs)

过滤稳定信号

高质量的数据生成需要在整个采集过程中保持稳定的信号。随着时间的推移，各种问题都可能导致信号发生变化，包括不稳定的流速、系统中引入空气或引入金属污染物。这些信号强度的变化在持续时间和强度方面可能有所不同，并且可能同时影响所有通道。为了检测受信号不稳定影响的采集区域，我们在散点图中将选定通道的信号显示为时间函数（图 1）。

# plot channels of interest vs. time
coi <- chs[c(dna[1], which(rowData(sce)$use_channel))]
plotScatter(sce, chs = c("Time", coi), label = "both") +
   labs(y = "expression") +
   scale_x_continuous(
     expression("Time ("*10^6~"ms)"),
     labels = function(u) u/1e6) +
   theme_bw(base_size = 8) + theme(
     aspect.ratio = 2/3,
     panel.grid = element_blank(),
     axis.text = element_text(color = "black"),
     strip.background = element_rect(fill = NA))

图1：DNA通道Ir191和测量抗原的41个通道对时间的散点图。

在本实验中，我们没有观察到与时间相关的信号不稳定性。在应剔除部分采集的情况下，可以通过在信号稳定的区域上手动选择完成，随后的留下来的子集只保留满足稳定条件(参数pop="+")内的细胞来完成。反之亦然，可以选择一个信号不稳定的区域，将其从SCE对象(pop="-")中移除。为了完整起见，我们展示了如何通过手动门控排除不稳定信号的区域：

# construct ’GatingSet’
ff <- sce2fcs(sce[dna, ], assay = "exprs")
gs <- GatingSet(flowSet(ff))

# apply rectangular gate to exclude unstable signal
min_t <- ...
max_t <- ...
gs_add_gating_method(
  gs, alias = "stable",
   pop = "-", parent = "root",
   dims = paste0("Time,", chs[dna[1]]),
   gating_method = "boundary",
   gating_args = sprintf("min=c(%s,0),max=c(%s,10)", min_t, max_t))

# plot scatter of DNA vs. Time
ggcyto(gs,
   aes_string("Time", chs[dna[1]])) +
   geom_hex(bins = 128) +
   geom_gate("stable") +
   facet_null() + theme_bw() +
   ggtitle(NULL) + theme(
     legend.position = "none",
     panel.grid = element_blank())

# subset selected events
sce <- sce[, gh_pop_get_indices(gs, "stable")]

标准化

对于CyTOF，必须标准化采集过程中由于灵敏度逐渐丧失而导致的信号漂移。一种最常用的策略是将样品与嵌入金属镧系元素的聚苯乙烯珠混合，从而可以在整个数据采集过程中监控仪器性能。这些珠子依次用于估计和校正信号的时间漂移。当必须在相同背景下分析独立实验时，还必须考虑由于仪器性能随时间的变化以及计划维护之间的间隔而产生的变化。在这种情况下，磁珠信号应以早期实验中的磁珠作为参考进行标准化。这可确保不同的实验标准化至同一水平，而与CyTOF的灵敏度无关。 目前的R实现可以通过CATALYST和premessa完成。CATALYST提供了基于磁珠的标准化的方法，可以自动识别磁珠单态(用于归一化)，以及磁珠-磁珠和细胞-磁珠二联体(剔除)，从而无需手动设门。具体实现如下：

• 最开始，确定磁珠在哪些events中具有最高信号
• 通过在最低磁珠和最高非磁珠通道信号之间利用分离截止来去除细胞-磁珠二联体
• 所有低于磁珠信号下界的events被移除(这些包括磁珠-磁珠和细胞-磁珠二联体)
• 通过对步骤2中确定的events应用默认中值median ± 5 mad来删除磁珠-磁珠二联体；剩余的磁珠events用于标准化

上述过程是通过函数normCytof()来实现的，它以SCE和一组控制归一化和输出格式的参数作为输入。这里，我们指定dna = 191(Ir191)和beads = "dvs"，对应DVS Science beads(镧系元素Ce140、Eu151、Eu153、Ho165、Lu175)。其次，我们提供了一组用于计算基线强度进行归一化的参考磁珠(参数norm_to)路径。最后，我们设置overwrite=FALSE来保留原始数据和归一化数据，并且remove_beads=TRUE来排除beads和二联体events：

# specify path to reference beads
ref_beads <- file.path("data", "normalization_beads.fcs")
# apply bead-based normalization
system.time(res <- normCytof(sce, beads = "dvs", dna = 191,
   norm_to = ref_beads, remove_beads = TRUE, overwrite = FALSE))

## Identifying beads...

## Computing normalization factors...

##  用户  系统  流逝 
## 21.39  0.79 22.21

当remove_beads = TRUE(默认值)时，normCytof()将返回一个包含过滤、beads和remove events的三个SCE列表，以及两个ggplot对象：

names(res)

## [1] "data"    "beads"   "removed" "scatter" "lines"

第一个SCE(res$data)包含过滤后的数据，其中“normed”保存着标准化后的表达值。剩下的两个SCE分别是包含任何被识别为beads(beads)的events和所有移除events(包括磁珠、磁珠-磁珠和细胞-磁珠二联体;removed)的数据子集；因此，磁珠beads本身就是被移除eventsremoved的子集。在这里，我们比较了每个子集包含的细胞数量和百分比：

# view no. of remaining, bead & removed events
ns <- sapply(res[1:3], ncol)
ps <- sprintf("%1.2f", ns/ncol(sce)*100)
data.frame(t(cbind("# events" = ns, "% of total" = ps)))

##              data beads removed
## # events   337525 27544   30627
## % of total  91.68  7.48    8.32

作为第一个质量控制图，res$scatter(图2)显示了磁珠通道(x轴)与DNA(y轴)的散点图，其中被识别为磁珠的events及其表达范围在颜色上突出显示；磁珠events应该具有较低的DNA强度(因为它们不是细胞)和所有通道都有的高信号强度。

res$scatter

图2：磁珠通道的散射图vsDNA。

其次，res$lines(图3)显示了归一化前后平滑的磁珠强度中位值；在归一化之前，它们通常随着时间的推移而减少，而在归一化之后，它们应该是近似不变的，并且以基线为中心。在我们的数据集中，归一化是基于先前获得的参考磁珠进行的。因此，基线值对应于参考磁珠的平均磁珠通道强度。如图3所示，归一化后的磁珠通道水平显著降低，表明在当前实验中具有更高的灵敏度。重要的是，信号随时间的轻微下降在归一化后不再存在。

res$lines

图3：标准化前后磁珠强度中位值与时间的关系。

为了评估CyTOF在采集过程中的灵敏度，并识别在仪器调优过程中未检测到的潜在问题，我们计算了识别为beads events的平均beads通道计数(res$beads子集)。其中通道对应于beads的逻辑向量存储在rowData列bead_ch下，我们可以使用该逻辑向量对counts进行子集划分，使其仅包含beads通道。

# compute mean bead channel counts for current experiment
is_bead <- rowData(res$beads)$bead_ch  # get bead channels
bead_cs <- counts(res$beads)[is_bead, ] # subset counts
rownames(bead_cs) <- chs[is_bead]    # use channels as names
(bead_ms <- rowMeans(bead_cs))      # compute means

##  Ce140Di  Eu151Di  Eu153Di  Ho165Di  Lu175Di 
## 2842.462 2111.367 2660.618 2538.095 2323.409

为了评估当前实验过程中的测量灵敏度，我们将上述计算得到的平均磁珠通道计数与元数据表ref_beads_count.csv中的7个先前获得的实验结果进行了比较。得到的箱线图(图4)表明当前实验的灵敏度相对较高，但在之前实验的范围内表现良好。

# read in reference mean bead channel counts
ref <- read.csv(file.path("data", "ref_bead_counts.csv"))

# join into single tidy data.frame
df <- bind_rows(ref, bead_ms, .id = "group")
df <- melt(df, id.var = "group")

# boxplot of reference vs. current experiment's mean bead channel counts
ggplot(df, aes(variable, value)) +
  geom_boxplot(data = df[df$group == 1, ]) +
  geom_point(data = df[df$group == 2, ],
    col = "red", pch = 4, stroke = 1) +
  labs(x = "bead channel", y = "mean count") +
  qc_theme + ggtitle(
    "QC on bead channel counts",
    "[-] = reference | x = current experiment")

图4：磁珠通道计数质量控制。

归一化后，我们用过滤后的子集覆盖输入数据集，该子集不再包含磁珠events，或者磁珠-磁珠和磁珠-细胞二联体：

sce <- res$data

结语

由于内容太多，第一部分只介绍了文件读取串联、过滤稳定信号和磁珠归一化这几个部分，后面还有单细胞反卷积single-cell deconvolution、溢出补偿spillover compensation以及碎片和双峰去除后的活细胞门控的内容，马上更新，敬请期待！