在《生信小记：映射和注释查询数据集》中，我们介绍了将scRNA-seq的数据映射到参考scRNA-seq数据集的方法。scRNA-seq是细胞类型注释的金标准，但有许多实验技术（如scACTA-seq、CyTOF等）并不直接测量基因的表达量，无法直接利用scRNA-seq参考数据集进行细胞类型的标注。为之奈何？

一种理论可行的方法是，找到一个已经使用其他模态数据完成细胞类型标注的参考数据集，并将查询数据集映射到这个参考数据集上。然而，其他模态的参考数据集通常规模较小、注释不够完善，且缺乏统一的标准，直接映射这些参考数据集产生的结果准确性并不能让人满意。但是，scRNA-seq技术由于发展早，当前已经形成了大量高质量、多组织、多物种的参考数据集，具有精细的多层次细胞类型注释。另外，scRNA-seq测量的是mRNA的转录水平，基于此的细胞类型定义在功能上更为直接，准确性也更高。把scRNA-seq的参考数据集“转译”为其他模态，以“转译”后的数据集作为参考数据集进行映射，或许是更好的选择。

那么，如何进行“转译”呢？假设现在有两个国家的人儿，他们的语言互不相通，如果他们之间想要相互交流，应当怎么办？答案是找到一位同时精通两国语言的翻译专家，或者是找到两个国家语言互译的“词典”，以此来跨过语言上的鸿沟、实现对话和交流。对于其他模态的数据集（查询数据集）和scRNA-seq参考数据集亦是如此，它们也需要一本语言互译的词典，使得scRNA-seq参考数据集能理解并整合查询数据集。这里的词典即为 桥接数据集（bridge datasets） ，其本质为包含同一细胞的“配对”多模态信息（包括查询数据集和参考数据集两种类型的标准数据）的数据集，使得算法可以学习不同模态数据之间的对应关系，从而将scRNA-seq参考数据集的细胞类型注释能力拓展到其他模态的数据集上，形成 扩展参考（extended reference）。使用扩展参考直接映射查询数据集，即可实现查询细胞集的细胞类型标注。

上述过程即为 桥接整合（bridge integration） ，即利用桥接数据集作为分子桥梁，将不直接测量基因表达的单细胞数据映射并整合到已注释的scRNA-seq参考数据集中的方法。在本篇文章中，我们将，以10X Multiome数据集作为桥接数据集，将人PBMC scATCA-seq查询数据集映射到已经注释的scRNA-seq参考数据集上，以实现细胞类型的注释。

一些背景知识介绍
scATAC-seq是一种测量染色质可及性的单细胞测序技术。根据开放状态，细胞内的染色质可以分为开放染色质和闭合染色质。前者结构松散，可以被转录因子等调控因子结合，继而启动基因的转录和表达，相应区域基因处于高表达状态，因此染色质可及性高；而后者结构紧密，调控因子难以结合，相应区域基因处于低表达状态，因此染色质可及性低。
在scATAC-seq的数据分析中，峰值是指基因组上染色质可及性明显升高的连续区域。每个峰值代表一个潜在的基因调控元件（如启动子、增强子等），其开放意味着它正在活跃地调控附近基因的转录。

和之前一样，我们先加载好所需的依赖包。

1
2
3
4
5
6
7
8
BiocManager::install("biovizBase")

library(Seurat)
library(Signac)
library(EnsDb.Hsapiens.v86)
library(dplyr)
library(ggplot2)
library(biovizBase)

加载桥接、查询与参考数据集

这一步中，我们会分别加载所需要的三个数据集，并进行预处理。

加载并预处理桥接数据集

我们所使用的桥接数据集包含了健康人12000个PBMCs scRNA-seq和scATCA-seq的数据，可以通过此处下载。不过，我们所下载到的scRNA-seq和scATAC-seq数据是分开储存的，因此我们需要把数据先“打包”好，构建多模态Seurat对象。

1
2
3
4
5
6
# 加载桥接数据集
inputdata.10x <- Read10X_h5("C:/Users/86150/Downloads/pbmc_granulocyte_sorted_10k_filtered_feature_bc_matrix.h5")

# 提取出其中的RNA和ATAC矩阵信息
rna_counts <- inputdata.10x$`Gene Expression`
atac_counts <- inputdata.10x$Peaks

创建RNA Assay

1
2
3
# 以scRNA-seq为基础创建Seurat对象，并获取线粒体基因的百分比
obj.multi <- CreateSeuratObject(counts = rna_counts)
obj.multi[["percent.mt"]] <- PercentageFeatureSet(obj.multi, pattern = "^MT-")

创建ATAC Assay（ChromatinAssay对象）

ChromatinAssay对象是一种特殊类型的Assay，用于存储与染色质可及性相关的内容，主要包括以下三大内容：

  • 计数矩阵（Counts）： 包含scATAC-seq数据中，每个基因组区间的染色质可及性的大小信息。
  • 基因组区间信息（Fragments）： 包括scATAC-seq数据中，每个基因组区间所在的位置信息。
  • 基因组注释信息（Annotation）： 包括基因组区间的功能性注释信息，标注了某段基因组区间是否属于、属于哪种基因调控元件（如启动子、增强子等）。

我们需要对上述三个内容分别进行预处理，然后对其进行整合为完整的ChromatinAssay对象。

计数矩阵（Counts）的预处理

人细胞中存在22对常染色体和1对性染色体。但除此之外，scATAC-seq还可以测得一些独立于前述染色体的DNA物质，包括线粒体DNA、已被测序但尚未精确定位到某条特定染色体上的DNA片段、来源不明或组装不完整的DNA片段等。但我们所研究的重心是核内染色体的基因表达，上述DNA片段的数据会为我们的分析引入大量噪声，所以应当予以去除。在这里，我们先筛出不属于标准染色体的基因组片段，在后续步骤予以去除。

1
2
3
grange.counts <- StringToGRanges(rownames(atac_counts), sep = c(":", "-"))
grange.use <- seqnames(grange.counts) %in% standardChromosomes(grange.counts)
atac_counts <- atac_counts[as.vector(grange.use), ]
  • rownames(atac_counts)以类似“染色体编号:区间起始位置-区间结束位置”的字符串的形式，存储着峰值所在的基因组位置信息。StringToGRanges函数可以将这些字符串转换为一个标准的基因组区间对象（即GRange对象），便于后续进行基因组操作。具体来说，我们最开始拿到的数据是一长串“染色体编号:区间起始位置-区间结束位置”的文本，而StringToGranges函数则可以帮助我们识别出，这段文本的哪一部分是染色质编号的信息、哪一部分是区间起始位置的信息、哪一部分是区间结束位置的信息、哪一部分是分隔这些信息的符号。
  • seqnames函数可以提取出每个峰值所在染色体编号。standardChromosomes函数则可以查询这个基因组所属类型的标准染色体列表（即chr1-22，chrX，chrY）。通过%in%这一管道符，可以检验每个峰值所在染色体是否属于标准染色体，是则返回TRUE，否则返回FALSE。
  • as.vector函数可以抹除grange.use中除了基因组区间外可能存在的额外信息，使其成为一个最基础、纯净的逻辑向量，避免后续数据分析中出现问题。

基因组注释信息（annotation）的预处理

Ensembl数据库包含了多个物种的基因组数据、注释和相关的生物学信息，其中包括人类基因组的功能性注释，即人类基因组中的基因、外显子、启动子等区域的位置。GetGRangesFromEnsDb函数可以帮助我们获得Ensembl数据库中人类基因组的注释信息，以GRange对象的形式储存。

不同来源的数据对染色体的编号方式不同，如Ensembl数据库将染色体命名为“1”、“2”等，而UCSC则会将染色体命名为“chr1”、“chr2”等。这里我们要“书同文”，使用seqlevesStyle函数将染色体的命名风格统一为“UCSC”格式。

另外，我们所需使用的数据，是基于第38版人类基因组测序的注释信息数据库，我们需要对此加以指定，以避免基因组版本不一致导致分析出错。

1
2
3
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86) # 获取人类基因组的注释信息
seqlevelsStyle(annotations) <- 'UCSC' # 统一注释风格
genome(annotations) <- "hg38" # 明确指定基因组版本为hg38

不过，如果你像我一样，在运行seqleveStyle函数时遇到了转换映射表下载失败的问题，可以试试手动修改染色体命名。

1
2
seqlevels(annotations) <- paste0("chr", seqlevels(annotations))
seqlevels(annotations) <- gsub("chrMT", "chrM", seqlevels(annotations)) #线粒体基因的处理

载入基因组区间信息（Fragments）

基因组区间信息不太需要做预处理，但有一点需要注意：在下载文件时，我们需要同时下载pbmc_granulocyte_sorted_10k_atac_fragments.tsv.gz和它的索引文件pbmc_granulocyte_sorted_10k_atac_fragments.tsv.gz.tbi，否则后面创建ChromatinAssay对象时会报错。

1
frag.file <- "D:/firstp/file/dictionary learning/pbmc_granulocyte_sorted_10k_atac_fragments.tsv.gz"

创建ChromatinAssay对象

咱使用CreateChromatinAssay函数来创建ChromatinAssay对象。

1
2
3
4
5
6
7
8
chrom_assay <- CreateChromatinAssay(
  counts = atac_counts,
  sep = c(":", "-"),
  genome = 'hg38',
  fragments = frag.file,
  min.cells = 10, # 这里是在做初步的质量控制，只保留在至少10个细胞中检测到的峰值。
  annotation = annotations
)

整合Seurat对象，并质量控制

1
2
3
4
5
6
7
8
9
obj.multi[["ATAC"]] <- chrom_assay
obj.multi <- subset(
  x = obj.multi,
  subset = nCount_ATAC < 7e4 &
    nCount_ATAC > 5e3 &
    nCount_RNA < 25000 &
    nCount_RNA > 1000 &
    percent.mt < 20
)

加载并预处理查询数据集

查询数据集的下载地址在这里。我们需要下载的数据包括：Fragments（TSV）、Fragments index（TBI）、Peak by cell matrix HDF5（raw）。

首先载入数据集。和前面处理桥接数据集一样，我们需要使用seqlevesStyle函数将染色体的命名风格统一为“UCSC”格式，再以数据集为基础创建ChromatinAssay对象。

1
2
3
4
5
6
7
8
9
10
11
atac_pbmc_data <- Read10X_h5(filename = "D:/firstp/file/dictionary learning/10k_PBMC_ATAC_nextgem_Chromium_X_filtered_peak_bc_matrix.h5") 
fragpath <- "D:/firstp/file/dictionary learning/10k_PBMC_ATAC_nextgem_Chromium_X_fragments.tsv.gz"

annotation <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
seqlevelsStyle(annotation) <- 'UCSC'
atac_pbmc_assay <- CreateChromatinAssay(
  counts = atac_pbmc_data,
  sep = c(":", "-"),
  fragments = fragpath,
  annotation = annotation
)

接下来，我们将重新量化查询数据集，使其与桥接数据集具有可比性。用人话来解释的话，在不同实验、不同批次、不同计数平台所测得的数据中，同一段染色质开放区域的位置和数量存在不同；因此我们需要利用FeatureMatrix函数使fragments数据集主动对齐features数据集，以便下一步分析。

1
2
3
4
5
requant_multiome_ATAC <- FeatureMatrix(
  fragments = Fragments(atac_pbmc_assay),
  features = granges(obj.multi[['ATAC']]),
  cells = Cells(atac_pbmc_assay)
)

最后，我们需要为量化好的查询数据集创建一个单独的ChromatinAssay对象，并将新的ChromatinAssay对象和原始ChromatinAssay对象放在一起，创建一个查询Seurat对象。（质量控制也不要忘了哦！）

1
2
3
4
5
6
7
8
9
ATAC_assay <- CreateChromatinAssay(
  counts = requant_multiome_ATAC,
  fragments = fragpath,
  annotation = annotation
)

obj.atac  <- CreateSeuratObject(counts = ATAC_assay,assay = 'ATAC')
obj.atac[['peak.orig']] <- atac_pbmc_assay
obj.atac <- subset(obj.atac, subset = nCount_ATAC < 7e4 & nCount_ATAC > 2000)

加载参考数据集

参考数据集在这里下载！下载后直接加载即可。

1