往期精選
使用Signac包進行單細胞ATAC-seq數(shù)據(jù)分析（一）：Analyzing PBMC scATAC-seq
使用Signac包進行單細胞ATAC-seq數(shù)據(jù)分析（二）：Motif analysis with Signac
使用Signac包進行單細胞ATAC-seq數(shù)據(jù)分析（三）：scATAC-seq data integration
使用Signac包進行單細胞ATAC-seq數(shù)據(jù)分析（四）：Merging objects
Cicero是一個用于分析
單細胞染色質可及性實驗的R工具包。Cicero的主要功能是使用單細胞染色質可及性數(shù)據(jù)，通過分析共同可及性來預測基因組中的順式調控相互作用（如增強子和啟動子之間的相互作用）。此外，Cicero包還擴展了Monocle的功能，能夠利用染色質可及性數(shù)據(jù)對單細胞進行聚類，排序和差異可及性分析。
Cicero包的簡介
Cicero包的主要功能是使用單細胞染色質可及性數(shù)據(jù)來預測基因組中更可能位于細胞核附近的區(qū)域。這可用于鑒定潛在的增強子-啟動子互作對，并了解全基因組區(qū)域內的順式相互作用的整體結構。
由于單細胞數(shù)據(jù)的稀疏性，細胞必須根據(jù)相似度進行聚合，以對數(shù)據(jù)中的各種技術因素進行可靠的校正。最終，Cicero根據(jù)用戶給定的距離，計算在指定距離內每個可及性peaks對之間的開放性，并給出"Cicero co-accessibility"得分，其分數(shù)介于-1和1之間，得分越高，表示更高的共可及性（co-accessibility）。
此外，Cicero包還提供了擴展工具包，可使用Monocle提供的框架來分析單細胞ATAC-seq數(shù)據(jù)。
Cicero包提供了兩種主要的分析功能：
- 構建和分析順式調控網(wǎng)絡。 Cicero通過分析共可及性（co-accessibility），來識別鑒定潛在的順式調控相互作用，并使用各種技術對其進行可視化和分析。
- 常規(guī)單細胞染色質可及性分析。 Cicero還擴展了Monocle包，以使用單細胞染色質可及性數(shù)據(jù)來進行差異可及性（differential accessibility）分析，細胞聚類與可視化，和發(fā)育軌跡重建。
Cicero包的安裝
通過Bioconductor進行安裝
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
# 安裝依賴包
BiocManager::install(c("Gviz", "GenomicRanges", "rtracklayer"))
BiocManager::install("cicero")
# 加載cicero包
library(cicero)
通過Github進行安裝
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
# 安裝依賴包
BiocManager::install(c("Gviz", "GenomicRanges", "rtracklayer"))
install.packages("devtools")
devtools::install_github("cole-trapnell-lab/cicero-release")
# 加載cicero包
library(cicero)
數(shù)據(jù)集的加載
Cicero將數(shù)據(jù)存儲在CellDataSet（CDS）類的對象中，該類繼承自Bioconductor的ExpressionSet類。我們可以使用以下三個函數(shù)來操作該對象：
- fData: 獲取feature的元信息
- pData: 獲取cell/sample的元信息
- exprs: 獲取cell-by-peak的count矩陣
為了修改CDS對象以保留染色質可及性數(shù)據(jù)而不是表達數(shù)據(jù)，Cicero使用peaks作為feature數(shù)據(jù)而不是基因或轉錄本。具體來說，許多Cicero函數(shù)需要形式為chr1_10390134_10391134的peaks值信息，如下所示：
Loading data from a simple sparse matrix format
Cicero可以讀取以簡單稀疏矩陣格式存儲的數(shù)據(jù)，這里以Cicero包自帶的一個小數(shù)據(jù)集cicero_data為例。
# 加載示例數(shù)據(jù)集
data(cicero_data)
# 查看示例數(shù)據(jù)
head(cicero_data)
Peak Cell Count
140 chr18_30209631_30210783 AGCGATAGGCGCTATGGTGGAATTCAGTCAGGACGT 4
150 chr18_45820294_45821666 AGCGATAGGTAGCAGCTATGGTAATCCTAGGCGAAG 2
185 chr18_32820116_32820994 TAATGCGCCGCTTATCGTTGGCAGCTCGGTACTGAC 2
266 chr18_41888433_41890138 AGCGATAGGCGCTATGGTGGAATTCAGTCAGGACGT 2
273 chr18_33038287_33039444 AGCGATAGGGTTATCGAACTCCATCGAGGTACTGAC 2
285 chr18_25533921_25534483 ATTACTCGAACGCGCAGAGGCGGAGGTCGTACTGAC 1
為了方便起見，Cicero提供了一個名為make_atac_cds的函數(shù)。該函數(shù)以稀疏矩陣格式將data.frame或文件的路徑作為輸入。具體來說，此文件應該是由三列組成的制表符分隔的文本文件。 第一列是peak的坐標，格式為“ chr10_100013372_100013596”，第二列是細胞的名稱，第三列是整數(shù)，表示細胞在該peak重疊的讀取次數(shù)，且此文件不應包含標題行，如下所示：
# 使用make_atac_cds函數(shù)將稀疏矩陣轉換為CDS對象
input_cds <- make_atac_cds(cicero_data, binarize = TRUE)
# 查看CDS對象
input_cds
CellDataSet (storageMode: environment)
assayData: 6146 features, 200 samples
element names: exprs
protocolData: none
phenoData
sampleNames: AGCGATAGAACGAATTCGGCGCAATGACCCTATCCT
AGCGATAGAAGTACGCGATCCGCGGACTGTACTGAC ...
TCTCGCGCTCTTGAGGTTTTATGACCAAATAGAGGC (200 total)
varLabels: cells Size_Factor num_genes_expressed
varMetadata: labelDescription
featureData
featureNames: chr18_10025_10225 chr18_10603_11103 ...
chr18_78015362_78016311 (6146 total)
fvarLabels: site_name chr ... num_cells_expressed (5 total)
fvarMetadata: labelDescription
experimentData: use 'experimentData(object)'
Annotation:
# 查看feature的metadata
head(fData(input_cds))
site_name chr bp1 bp2
chr18_10025_10225 chr18_10025_10225 18 10025 10225
chr18_10603_11103 chr18_10603_11103 18 10603 11103
chr18_11604_13986 chr18_11604_13986 18 11604 13986
chr18_49557_50057 chr18_49557_50057 18 49557 50057
chr18_50240_50740 chr18_50240_50740 18 50240 50740
chr18_104385_104585 chr18_104385_104585 18 104385 104585
num_cells_expressed
chr18_10025_10225 5
chr18_10603_11103 1
chr18_11604_13986 9
chr18_49557_50057 2
chr18_50240_50740 2
chr18_104385_104585 1
# 查看cell的metadata
head(pData(input_cds))
cells
AGCGATAGAACGAATTCGGCGCAATGACCCTATCCT AGCGATAGAACGAATTCGGCGCAATGACCCTATCCT
AGCGATAGAAGTACGCGATCCGCGGACTGTACTGAC AGCGATAGAAGTACGCGATCCGCGGACTGTACTGAC
AGCGATAGAATACGATAAGGCCGTCAACTAATCTTA AGCGATAGAATACGATAAGGCCGTCAACTAATCTTA
AGCGATAGATTATGCAAGCCAGTACTTGCCTATCCT AGCGATAGATTATGCAAGCCAGTACTTGCCTATCCT
AGCGATAGCAGACTAAGGGGAATTCAGTGGCTCTGA AGCGATAGCAGACTAAGGGGAATTCAGTGGCTCTGA
AGCGATAGCCGTATGATTAGATCTTGGTCAGGACGT AGCGATAGCCGTATGATTAGATCTTGGTCAGGACGT
Size_Factor num_genes_expressed
AGCGATAGAACGAATTCGGCGCAATGACCCTATCCT NA 290
AGCGATAGAAGTACGCGATCCGCGGACTGTACTGAC NA 490
AGCGATAGAATACGATAAGGCCGTCAACTAATCTTA NA 253
AGCGATAGATTATGCAAGCCAGTACTTGCCTATCCT NA 181
AGCGATAGCAGACTAAGGGGAATTCAGTGGCTCTGA NA 85
AGCGATAGCCGTATGATTAGATCTTGGTCAGGACGT NA 251
# 查看cell-by-peak的count矩陣
head(exprs(input_cds))
6 x 200 sparse Matrix of class "dgCMatrix"
[[ suppressing 32 column names 'AGCGATAGAACGAATTCGGCGCAATGACCCTATCCT', 'AGCGATAGAAGTACGCGATCCGCGGACTGTACTGAC', 'AGCGATAGAATACGATAAGGCCGTCAACTAATCTTA' ... ]]
chr18_10025_10225 . . . . . . . . . . . . . . . . . . . . . . . . . . . .
chr18_10603_11103 . . . . . . . . . . . . . . . . . . . . . . . . . . . .
chr18_11604_13986 . 1 . . . . . . . . . . . . . . . . . . . . . . . . . .
chr18_49557_50057 . . 1 . . . . . . . . . . . . . . . . . . . . . . . . .
chr18_50240_50740 . . . . . . . . . . . . . . . . . . . . . . . . . . . .
chr18_104385_104585 . . . . . . . . . . . . . . . . . . . . . . . . . . . .
chr18_10025_10225 . . . . ......
chr18_10603_11103 . . . . ......
chr18_11604_13986 . . . . ......
chr18_49557_50057 . . . . ......
chr18_50240_50740 . . . . ......
chr18_104385_104585 . 1 . . ......
.....suppressing 168 columns in show(); maybe adjust 'options(max.print= *, width = *)'
..............................
Loading 10X scATAC-seq data
如果我們的scATAC-seq數(shù)據(jù)來自10x Genomics平臺，使用Cell Ranger ATAC軟件處理將數(shù)據(jù)輸出到一個名為filtered_peak_bc_matrix的文件夾中，可以通過以下方式將數(shù)據(jù)轉換為CDS對象。
加載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的metadata
# format cell info
cellinfo <- read.table("filtered_peak_bc_matrix/barcodes.tsv")
row.names(cellinfo) <- cellinfo$V1
names(cellinfo) <- "cells"
# 加載peak的metadata
# 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
row.names(indata) <- row.names(peakinfo)
colnames(indata) <- row.names(cellinfo)
# 使用newCellDataSet函數(shù)構建CDS對象
# make CDS
fd <- methods::new("AnnotatedDataFrame", data = peakinfo)
pd <- methods::new("AnnotatedDataFrame", data = cellinfo)
input_cds <- suppressWarnings(newCellDataSet(indata,
phenoData = pd,
featureData = fd,
expressionFamily=VGAM::binomialff(),
lowerDetectionLimit=0))
input_cds@expressionFamily@vfamily <- "binomialff"
input_cds <- monocle::detectGenes(input_cds)
# 數(shù)據(jù)初步過濾
#Ensure there are no peaks included with zero reads
input_cds <- input_cds[Matrix::rowSums(exprs(input_cds)) != 0,]
參考來源：https://cole-trapnell-lab.github.io/cicero-release/docs/
