本文可在http://xuzhougeng.top/免費(fèi)閱讀原文

Cicero是一個(gè)單細(xì)胞染色質(zhì)可及性實(shí)驗(yàn)可視化R包她倘。Cicero的主要功能就是使用單細(xì)胞染色質(zhì)可及性數(shù)據(jù)通過分析共開放去預(yù)測基因組上順式調(diào)節(jié)作用(cis-regulatory interactions)窿克，例如增強(qiáng)子和啟動(dòng)子湿痢。Cicero能夠利用染色質(zhì)開放性進(jìn)行單細(xì)胞聚類末盔，排序和差異可及性分析支竹。關(guān)于Cicero的算法原理眷细，參考他們發(fā)表的文章.

簡介

Cicero的主要目標(biāo)是使用單細(xì)胞染色質(zhì)可及性數(shù)據(jù)去預(yù)測基因組上那些在細(xì)胞核中存在物理鄰近的區(qū)域太抓。這能夠用于鑒定潛在的增強(qiáng)子-啟動(dòng)子對(duì)，對(duì)基因組順式作用有一個(gè)總體了解。

由于單細(xì)胞數(shù)據(jù)的稀疏性肺孵，細(xì)胞必須根據(jù)相似度進(jìn)行聚合匀借，才能夠?qū)?shù)據(jù)中的多種技術(shù)因素進(jìn)行可靠糾正。

最終平窘，Cicero給出了"Cicero 共開放"得分吓肋，分?jǐn)?shù)在-1和1之間，用于評(píng)估用于給定距離中每個(gè)可及peak對(duì)之間的共開放性瑰艘，分?jǐn)?shù)越高是鬼，越可能是共開放。

此外紫新，Monocle3的框架讓Cicero能對(duì)單細(xì)胞ATAC-seq實(shí)驗(yàn)完成Monocle3上的聚類均蜜，擬時(shí)間等分析。

Cicero主要提供了兩種核心分析:

• 構(gòu)建和分析順式調(diào)控網(wǎng)絡(luò).
• 常規(guī)單細(xì)胞染色質(zhì)可及性分析:

Cicero安裝和加載

Cicero運(yùn)行在R語言分析環(huán)境中芒率。盡管能夠從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)

數(shù)據(jù)加載

Cicero以cell_data_set(CDS)對(duì)象存放數(shù)據(jù)偶芍，該對(duì)象繼承自Bioconductor的SingleCellExperiment對(duì)象充择。我們可以用下面三個(gè)函數(shù)來操縱數(shù)據(jù)

• fData: 獲取feature的元信息, 這里的featureI指的是peak
• pData: 獲取cell/sample的元信息
• assay: 獲取cell-by-peak的count矩陣

除了CDS外，還需要基因坐標(biāo)信息和基因注釋信息匪蟀。

關(guān)于坐標(biāo)信息椎麦，以mouse.mm9.genome為例，它是一個(gè)數(shù)據(jù)框材彪，只有兩列观挎，一列是染色體編號(hào)，另一列對(duì)應(yīng)的長度

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也有現(xiàn)成的TxDb段化，可以從中提取注釋信息键兜。

簡單稀疏矩陣格式

以一個(gè)測試數(shù)據(jù)集介紹如何加載簡單稀疏矩陣格式

temp <- textConnection(readLines(gzcon(url("http://staff.washington.edu/hpliner/data/kidney_data.txt.gz"))))
cicero_data <- read.table(temp)

這里的簡單稀疏矩陣格式指的是數(shù)據(jù)以三列進(jìn)行存放，第一列是peak位置信息穗泵，第二列是樣本信息普气，第三列則是樣本中有多少fragment和該peak重疊。

chr1_4756552_4757256    GAATTCGTACCAGGCGCATTGGTAGTCGGGCTCTGA    1

對(duì)于這種格式佃延，Cicero提供了一個(gè)函數(shù)make_atac_cds, 用于構(gòu)建一個(gè)有效的cell_data_set對(duì)象现诀，用于下游分析，輸入既可以是一個(gè)數(shù)據(jù)框履肃，也可以是一個(gè)文件路徑仔沿。

input_cds <- make_atac_cds(cicero_data, binarize = TRUE)

10X scATAC-seq 數(shù)據(jù)

如果數(shù)據(jù)已經(jīng)經(jīng)過10X scATAC-seq處理，那么結(jié)果中的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元數(shù)據(jù)

# format cell info
cellinfo <- read.table("filtered_peak_bc_matrix/barcodes.tsv")
row.names(cellinfo) <- cellinfo$V1
names(cellinfo) <- "cells"

加載peak元數(shù)據(jù)

# 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元數(shù)據(jù))和列名(cell元數(shù)據(jù))

row.names(indata) <- row.names(peakinfo)
colnames(indata) <- row.names(cellinfo)

然后用new_cell_data_set根據(jù)peak元數(shù)據(jù)封锉，cell元數(shù)據(jù)構(gòu)建出一個(gè)cell_data_set對(duì)象

# make CDS
input_cds <-  suppressWarnings(new_cell_data_set(indata,
    cell_metadata = cellinfo,
    gene_metadata = peakinfo))

之后用detect_genes統(tǒng)計(jì)，對(duì)于每個(gè)基因在多少細(xì)胞中的表達(dá)量超過了給定閾值。

input_cds <- monocle3::detect_genes(input_cds)

過濾沒有細(xì)胞開放的peak

input_cds <- input_cds[Matrix::rowSums(exprs(input_cds)) != 0,] 

構(gòu)建順式作用網(wǎng)絡(luò)

運(yùn)行Cicero

構(gòu)建Cicero CDS對(duì)象

因?yàn)閱渭?xì)胞染色質(zhì)開放數(shù)據(jù)特別稀疏成福，為了能夠比較準(zhǔn)確的估計(jì)共開放得分碾局，需要將一些比較相近的細(xì)胞進(jìn)行聚合，得到一個(gè)相對(duì)比較致密的count數(shù)據(jù)奴艾。Cicero使用KNN算法構(gòu)建細(xì)胞的相互重疊集净当。而細(xì)胞之間距離關(guān)系則是根據(jù)降維后坐標(biāo)信息。降維方法有很多種蕴潦，這里以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")

此時(shí)用reducedDimNames(input_cds)就會(huì)發(fā)現(xiàn)它多了LSI和UMAP兩個(gè)信息像啼，我們可以用reducedDims來提取UMAP坐標(biāo)信息。

umap_coords <- reducedDims(input_cds)$UMAP

make_cicero_cds函數(shù)就需要其中UMAP坐標(biāo)信息潭苞。

注1: 假如你已經(jīng)有了UMAP的坐標(biāo)信息忽冻，那么就不需要運(yùn)行preprocess_cds和reduce_dimension， 直接提供UMAP坐標(biāo)信息給make_cicero_cds就可以了此疹。

注2: umap_coords的行名需要和pData得到表格中的行名一樣, 即all(row.names(umap_coords) == row.names(pData(input_cds)))結(jié)果為TRUE僧诚。

cicero_cds <- make_cicero_cds(input_cds, reduced_coordinates = umap_coords)

此時(shí)的cicero_cds是沒有UMAP降維信息。

運(yùn)行Cicero

Cicero包的主要功能就是預(yù)測基因組上位點(diǎn)之間的共開放. 有兩種方式可以獲取該信息

• run_cicero: 一步到位秀菱。默認(rèn)參數(shù)適用于人類和小鼠振诬，但是其他物種蹭睡，參考此處
• 單獨(dú)運(yùn)行每個(gè)函數(shù)衍菱，這樣子獲取中間的信息。

#測試部分?jǐn)?shù)據(jù)
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
## 全部數(shù)據(jù), 時(shí)間真久
conns <- run_cicero(cicero_cds, mouse.mm9.genome, sample_num = 2)
nrows(conns) #59741738

其中conns存放的是peak之間共開放得分肩豁。

Cicero Connection可視化