Cicero:一个单细胞染色质可及性实验可视化R包
本文可在http://xuzhougeng.top/免费阅读原文
Cicero是一个单细胞染色质可及性实验可视化R包。Cicero的主要功能就是使用单细胞染色质可及性数据通过分析共开放去预测基因组上顺式调节作用(cis-regulatory interactions)，例如增强子和启动子。Cicero能够利用染色质开放性进行单细胞聚类，排序和差异可及性分析。关于Cicero的算法原理，参考他们发表的文章.
简介
Cicero的主要目标是使用单细胞染色质可及性数据去预测基因组上那些在细胞核中存在物理邻近的区域。这能够用于鉴定潜在的增强子-启动子对，对基因组顺式作用有一个总体了解。
由于单细胞数据的稀疏性，细胞必须根据相似度进行聚合，才能够对数据中的多种技术因素进行可靠纠正。
最终，Cicero给出了"Cicero 共开放"得分，分数在-1和1之间，用于评估用于给定距离中每个可及peak对之间的共开放性，分数越高，越可能是共开放。
此外，Monocle3的框架让Cicero能对单细胞ATAC-seq实验完成Monocle3上的聚类，拟时间等分析。
Cicero主要提供了两种核心分析:
- 构建和分析顺式调控网络.
- 常规单细胞染色质可及性分析:
Cicero安装和加载
Cicero运行在R语言分析环境中。尽管能够从Bioconductor安装，但是推荐从GitHub上安装。
install.packages("BiocManager")
BiocManager::install(c("Gviz", "GenomicRanges", "rtracklayer"))
install.packages("devtools")
devtools::install_github("cole-trapnell-lab/cicero-release", ref = "monocle3")
加载Cicero
library(cicero)
数据加载
Cicero以cell_data_set(CDS)对象存放数据，该对象继承自Bioconductor的SingleCellExperiment对象。我们可以用下面三个函数来操纵数据
- fData: 获取feature的元信息, 这里的featureI指的是peak
- pData: 获取cell/sample的元信息
- assay: 获取cell-by-peak的count矩阵
除了CDS外，还需要基因坐标信息和基因注释信息。
关于坐标信息，以mouse.mm9.genome为例，它是一个数据框，只有两列，一列是染色体编号，另一列对应的长度
data("mouse.mm9.genome")
基因组注释信息可以用rtracklayer::readGFF加载
temp <- tempfile()
download.file("ftp://ftp.ensembl.org/pub/release-65/gtf/mus_musculus/Mus_musculus.NCBIM37.65.gtf.gz", temp)
gene_anno <- rtracklayer::readGFF(temp)
unlink(temp)
Bioconductor也有现成的TxDb，可以从中提取注释信息。
简单稀疏矩阵格式
以一个测试数据集介绍如何加载简单稀疏矩阵格式
temp <- textConnection(readLines(gzcon(url("http://staff.washington.edu/hpliner/data/kidney_data.txt.gz"))))
cicero_data <- read.table(temp)
这里的简单稀疏矩阵格式指的是数据以三列进行存放，第一列是peak位置信息，第二列是样本信息，第三列则是样本中有多少fragment和该peak重叠。
chr1_4756552_4757256 GAATTCGTACCAGGCGCATTGGTAGTCGGGCTCTGA 1
对于这种格式，Cicero提供了一个函数make_atac_cds, 用于构建一个有效的cell_data_set对象，用于下游分析，输入既可以是一个数据框，也可以是一个文件路径。
input_cds <- make_atac_cds(cicero_data, binarize = TRUE)
10X scATAC-seq 数据
如果数据已经经过10X scATAC-seq处理，那么结果中的filtered_peak_bc_matrix里就有我们需要的信息
加载cell-by-peak count矩阵，
# read in matrix data using the Matrix package
indata <- Matrix::readMM("filtered_peak_bc_matrix/matrix.mtx")
# binarize the matrix
indata@x[indata@x > 0] <- 1
加载cell元数据
# format cell info
cellinfo <- read.table("filtered_peak_bc_matrix/barcodes.tsv")
row.names(cellinfo) <- cellinfo$V1
names(cellinfo) <- "cells"
加载peak元数据
# format peak info
peakinfo <- read.table("filtered_peak_bc_matrix/peaks.bed")
names(peakinfo) <- c("chr", "bp1", "bp2")
peakinfo$site_name <- paste(peakinfo$chr, peakinfo$bp1, peakinfo$bp2, sep="_")
row.names(peakinfo) <- peakinfo$site_name
为cell-by-peak的count矩阵增加行名(peak元数据)和列名(cell元数据)
row.names(indata) <- row.names(peakinfo)
colnames(indata) <- row.names(cellinfo)
然后用new_cell_data_set根据peak元数据，cell元数据构建出一个cell_data_set对象
# make CDS
input_cds <- suppressWarnings(new_cell_data_set(indata,
cell_metadata = cellinfo,
gene_metadata = peakinfo))
之后用detect_genes统计，对于每个基因在多少细胞中的表达量超过了给定阈值。
input_cds <- monocle3::detect_genes(input_cds)
过滤没有细胞开放的peak
input_cds <- input_cds[Matrix::rowSums(exprs(input_cds)) != 0,]
构建顺式作用网络
运行Cicero
构建Cicero CDS对象
因为单细胞染色质开放数据特别稀疏，为了能够比较准确的估计共开放得分，需要将一些比较相近的细胞进行聚合，得到一个相对比较致密的count数据。Cicero使用KNN算法构建细胞的相互重叠集。而细胞之间距离关系则是根据降维后坐标信息。降维方法有很多种，这里以UMAP为例
input_cds <- detect_genes(input_cds)
input_cds <- estimate_size_factors(input_cds)
# PCA或LSI线性降维
input_cds <- preprocess_cds(input_cds, method = "LSI")
# UMAP非线性降维
input_cds <- reduce_dimension(input_cds, reduction_method = 'UMAP',
preprocess_method = "LSI")
此时用reducedDimNames(input_cds)就会发现它多了LSI和UMAP两个信息，我们可以用reducedDims来提取UMAP坐标信息。
umap_coords <- reducedDims(input_cds)$UMAP
make_cicero_cds函数就需要其中UMAP坐标信息。
注1: 假如你已经有了UMAP的坐标信息，那么就不需要运行preprocess_cds和reduce_dimension， 直接提供UMAP坐标信息给make_cicero_cds就可以了。
注2: umap_coords的行名需要和pData得到表格中的行名一样, 即all(row.names(umap_coords) == row.names(pData(input_cds)))结果为TRUE。
cicero_cds <- make_cicero_cds(input_cds, reduced_coordinates = umap_coords)
此时的cicero_cds是没有UMAP降维信息。
运行Cicero
Cicero包的主要功能就是预测基因组上位点之间的共开放. 有两种方式可以获取该信息
-
run_cicero: 一步到位。默认参数适用于人类和小鼠，但是其他物种，参考此处 - 单独运行每个函数，这样子获取中间的信息。
#测试部分数据
data("mouse.mm9.genome")
sample_genome <- subset(mouse.mm9.genome, V1 == "chr2")
sample_genome$V2[1] <- 10000000
conns <- run_cicero(cicero_cds, sample_genome, sample_num = 2)
nrow(conns) # 212302
## 全部数据, 时间真久
conns <- run_cicero(cicero_cds, mouse.mm9.genome, sample_num = 2)
nrows(conns) #59741738
其中conns存放的是peak之间共开放得分。