Analyzing PBMC scATAC-seq
Compiled: April 28, 2021
Source:vignettes/pbmc_vignette.Rmd
pbmc_vignette.RmdFor this tutorial, we will be analyzing a single-cell ATAC-seq dataset of human peripheral blood mononuclear cells (PBMCs) provided by 10x Genomics. The following files are used in this vignette, all available through the 10x Genomics website:
- The Raw data
- The Metadata
- The fragments file
- The fragments file index
To download all the required files, you can run the following lines in a shell:
wget https://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5
wget https://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_singlecell.csv
wget https://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_fragments.tsv.gz
wget https://cf.10xgenomics.com/samples/cell-atac/1.0.1/atac_v1_pbmc_10k/atac_v1_pbmc_10k_fragments.tsv.gz.tbiFirst load in Signac, Seurat, and some other packages we will be using for analyzing human data.
library(Signac)
library(Seurat)
library(GenomeInfoDb)
library(EnsDb.Hsapiens.v75)
library(ggplot2)
library(patchwork)
set.seed(1234)Pre-processing workflow
When pre-processing chromatin data, Signac uses information from two related input files, both of which can be created using CellRanger:
Peak/Cell matrix. This is analogous to the gene expression count matrix used to analyze single-cell RNA-seq. However, instead of genes, each row of the matrix represents a region of the genome (a peak), that is predicted to represent a region of open chromatin. Each value in the matrix represents the number of Tn5 integration sites for each single barcode (i.e. a cell) that map within each peak. You can find more detail on the 10X Website.
Fragment file. This represents a full list of all unique fragments across all single cells. It is a substantially larger file, is slower to work with, and is stored on-disk (instead of in memory). However, the advantage of retaining this file is that it contains all fragments associated with each single cell, as opposed to only fragments that map to peaks. More information about the fragment file can be found on the 10x Genomics website or on the sinto website.
We start by creating a Seurat object using the peak/cell matrix and cell metadata generated by cellranger-atac, and store the path to the fragment file on disk in the Seurat object:
counts <- Read10X_h5(filename = "/home/stuartt/github/chrom/vignette_data/atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5")
metadata <- read.csv(
file = "/home/stuartt/github/chrom/vignette_data/atac_v1_pbmc_10k_singlecell.csv",
header = TRUE,
row.names = 1
)
chrom_assay <- CreateChromatinAssay(
counts = counts,
sep = c(":", "-"),
genome = 'hg19',
fragments = '/home/stuartt/github/chrom/vignette_data/atac_v1_pbmc_10k_fragments.tsv.gz',
min.cells = 10,
min.features = 200
)
pbmc <- CreateSeuratObject(
counts = chrom_assay,
assay = "peaks",
meta.data = metadata
)
pbmc## An object of class Seurat
## 87561 features across 8728 samples within 1 assay
## Active assay: peaks (87561 features, 0 variable features)The ATAC-seq data is stored using a custom assay, the ChromatinAssay. This enables some specialized functions for analysing genomic single-cell assays such as scATAC-seq. By printing the assay we can see some of the additional information that can be contained in the ChromatinAssay, including motif information, gene annotations, and genome information.
pbmc[['peaks']]## ChromatinAssay data with 87561 features for 8728 cells
## Variable features: 0
## Genome: hg19
## Annotation present: FALSE
## Motifs present: FALSE
## Fragment files: 1For example, we can call granges on a Seurat object with a ChromatinAssay set as the active assay (or on a ChromatinAssay) to see the genomic ranges associated with each feature in the object. See the object interaction vignette for more information about the ChromatinAssay class.
granges(pbmc)## GRanges object with 87561 ranges and 0 metadata columns:
## seqnames ranges strand
## <Rle> <IRanges> <Rle>
## [1] chr1 565107-565550 *
## [2] chr1 569174-569639 *
## [3] chr1 713460-714823 *
## [4] chr1 752422-753038 *
## [5] chr1 762106-763359 *
## ... ... ... ...
## [87557] chrY 58993392-58993760 *
## [87558] chrY 58994571-58994823 *
## [87559] chrY 58996352-58997331 *
## [87560] chrY 59001782-59002175 *
## [87561] chrY 59017143-59017246 *
## -------
## seqinfo: 24 sequences from an unspecified genome; no seqlengthsWe can also add gene annotations to the pbmc object for the human genome. This will allow downstream functions to pull the gene annotation information directly from the object.
# extract gene annotations from EnsDb
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v75)
# change to UCSC style since the data was mapped to hg19
seqlevelsStyle(annotations) <- 'UCSC'
genome(annotations) <- "hg19"
# add the gene information to the object
Annotation(pbmc) <- annotationsComputing QC Metrics
We can now compute some QC metrics for the scATAC-seq experiment. We currently suggest the following metrics below to assess data quality. As with scRNA-seq, the expected range of values for these parameters will vary depending on your biological system, cell viability, and other factors.
Nucleosome banding pattern: The histogram of DNA fragment sizes (determined from the paired-end sequencing reads) should exhibit a strong nucleosome banding pattern corresponding to the length of DNA wrapped around a single nucleosome. We calculate this per single cell, and quantify the approximate ratio of mononucleosomal to nucleosome-free fragments (stored as
nucleosome_signal)Transcriptional start site (TSS) enrichment score. The ENCODE project has defined an ATAC-seq targeting score based on the ratio of fragments centered at the TSS to fragments in TSS-flanking regions (see https://www.encodeproject.org/data-standards/terms/). Poor ATAC-seq experiments typically will have a low TSS enrichment score. We can compute this metric for each cell with the
TSSEnrichment()function, and the results are stored in metadata under the column nameTSS.enrichment.Total number of fragments in peaks: A measure of cellular sequencing depth / complexity. Cells with very few reads may need to be excluded due to low sequencing depth. Cells with extremely high levels may represent doublets, nuclei clumps, or other artefacts.
Fraction of fragments in peaks: Represents the fraction of all fragments that fall within ATAC-seq peaks. Cells with low values (i.e. <15-20%) often represent low-quality cells or technical artifacts that should be removed. Note that this value can be sensitive to the set of peaks used.
Ratio reads in genomic blacklist regions The ENCODE project has provided a list of blacklist regions, representing reads which are often associated with artefactual signal. Cells with a high proportion of reads mapping to these areas (compared to reads mapping to peaks) often represent technical artifacts and should be removed. ENCODE blacklist regions for human (hg19 and GRCh38), mouse (mm10), Drosophila (dm3), and C. elegans (ce10) are included in the Signac package.
Note that the last three metrics can be obtained from the output of CellRanger (which is stored in the object metadata), but can also be calculated for non-10x datasets using Signac (more information at the end of this document).
# compute nucleosome signal score per cell
pbmc <- NucleosomeSignal(object = pbmc)
# compute TSS enrichment score per cell
pbmc <- TSSEnrichment(object = pbmc, fast = FALSE)
# add blacklist ratio and fraction of reads in peaks
pbmc$pct_reads_in_peaks <- pbmc$peak_region_fragments / pbmc$passed_filters * 100
pbmc$blacklist_ratio <- pbmc$blacklist_region_fragments / pbmc$peak_region_fragmentsWe can inspect the TSS enrichment scores by grouping the cells based on the score and plotting the accessibility signal over all TSS sites. Setting the fast=TRUE option in TSSEnrichment() will only compute the TSS enrichment score without storing the entire cell by position matrix of Tn5 insertion frequency for each cell, and can save memory. However, setting fast=TRUE will not allow downstream plotting of the TSS enrichment signal for different groups of cells using the TSSPlot() function, shown here:
pbmc$high.tss <- ifelse(pbmc$TSS.enrichment > 2, 'High', 'Low')
TSSPlot(pbmc, group.by = 'high.tss') + NoLegend()
We can also look at the fragment length periodicity for all the cells, and group by cells with high or low nucleosomal signal strength. You can see that cells that are outliers for the mononucleosomal / nucleosome-free ratio (based on the plots above) have different nucleosomal banding patterns. The remaining cells exhibit a pattern that is typical for a successful ATAC-seq experiment.
pbmc$nucleosome_group <- ifelse(pbmc$nucleosome_signal > 4, 'NS > 4', 'NS < 4')
FragmentHistogram(object = pbmc, group.by = 'nucleosome_group')
VlnPlot(
object = pbmc,
features = c('pct_reads_in_peaks', 'peak_region_fragments',
'TSS.enrichment', 'blacklist_ratio', 'nucleosome_signal'),
pt.size = 0.1,
ncol = 5
)
Finally we remove cells that are outliers for these QC metrics.
pbmc <- subset(
x = pbmc,
subset = peak_region_fragments > 3000 &
peak_region_fragments < 20000 &
pct_reads_in_peaks > 15 &
blacklist_ratio < 0.05 &
nucleosome_signal < 4 &
TSS.enrichment > 2
)
pbmc## An object of class Seurat
## 87561 features across 7060 samples within 1 assay
## Active assay: peaks (87561 features, 0 variable features)Normalization and linear dimensional reduction
Normalization: Signac performs term frequency-inverse document frequency (TF-IDF) normalization. This is a two-step normalization procedure, that both normalizes across cells to correct for differences in cellular sequencing depth, and across peaks to give higher values to more rare peaks.
Feature selection: The low dynamic range of scATAC-seq data makes it challenging to perform variable feature selection, as we do for scRNA-seq. Instead, we can choose to use only the top n% of features (peaks) for dimensional reduction, or remove features present in less than n cells with the
FindTopFeatures()function. Here, we will all features, though we note that we see very similar results when using only a subset of features (try setting min.cutoff to ‘q75’ to use the top 25% all peaks), with faster runtimes. Features used for dimensional reduction are automatically set asVariableFeatures()for the Seurat object by this function.Dimension reduction: We next run singular value decomposition (SVD) on the TD-IDF matrix, using the features (peaks) selected above. This returns a reduced dimension representation of the object (for users who are more familiar with scRNA-seq, you can think of this as analogous to the output of PCA).
The combined steps of TF-IDF followed by SVD are known as latent semantic indexing (LSI), and were first introduced for the analysis of scATAC-seq data by Cusanovich et al. 2015.
pbmc <- RunTFIDF(pbmc)
pbmc <- FindTopFeatures(pbmc, min.cutoff = 'q0')
pbmc <- RunSVD(pbmc)The first LSI component often captures sequencing depth (technical variation) rather than biological variation. If this is the case, the component should be removed from downstream analysis. We can assess the correlation between each LSI component and sequencing depth using the DepthCor() function:
DepthCor(pbmc)
Here we see there is a very strong correlation between the first LSI component and the total number of counts for the cell, so we will perform downstream steps without this component.
Non-linear dimension reduction and clustering
Now that the cells are embedded in a low-dimensional space, we can use methods commonly applied for the analysis of scRNA-seq data to perform graph-based clustering and non-linear dimension reduction for visualization. The functions RunUMAP(), FindNeighbors(), and FindClusters() all come from the Seurat package.
pbmc <- RunUMAP(object = pbmc, reduction = 'lsi', dims = 2:30)
pbmc <- FindNeighbors(object = pbmc, reduction = 'lsi', dims = 2:30)
pbmc <- FindClusters(object = pbmc, verbose = FALSE, algorithm = 3)
DimPlot(object = pbmc, label = TRUE) + NoLegend()
Create a gene activity matrix
The UMAP visualization reveals the presence of multiple cell groups in human blood. If you are familiar with our scRNA-seq analyses of PBMC, you may even recognize the presence of certain myeloid and lymphoid populations in the scATAC-seq data. However, annotating and interpreting clusters is more challenging in scATAC-seq data as much less is known about the functional roles of noncoding genomic regions than is known about protein coding regions (genes).
However, we can try to quantify the activity of each gene in the genome by assessing the chromatin accessibility associated with each gene, and create a new gene activity assay derived from the scATAC-seq data. Here we will use a simple approach of summing the fragments intersecting the gene body and promoter region (we also recommend exploring the Cicero tool, which can accomplish a similar goal, and we provide a vignette showing how to run Cicero within a Signac workflow here).
To create a gene activity matrix, we extract gene coordinates and extend them to include the 2 kb upstream region (as promoter accessibility is often correlated with gene expression). We then count the number of fragments for each cell that map to each of these regions, using the using the FeatureMatrix() function. These steps are automatically performed by the GeneActivity() function:
gene.activities <- GeneActivity(pbmc)
# add the gene activity matrix to the Seurat object as a new assay and normalize it
pbmc[['RNA']] <- CreateAssayObject(counts = gene.activities)
pbmc <- NormalizeData(
object = pbmc,
assay = 'RNA',
normalization.method = 'LogNormalize',
scale.factor = median(pbmc$nCount_RNA)
)Now we can visualize the activities of canonical marker genes to help interpret our ATAC-seq clusters. Note that the activities will be much noisier than scRNA-seq measurements. This is because they represent measurements from sparse chromatin data, and because they assume a general correspondence between gene body/promoter accessibility and gene expression which may not always be the case. Nonetheless, we can begin to discern populations of monocytes, B, T, and NK cells based on these gene activity profiles. However, further subdivision of these cell types is challenging based on supervised analysis alone.
DefaultAssay(pbmc) <- 'RNA'
FeaturePlot(
object = pbmc,
features = c('MS4A1', 'CD3D', 'LEF1', 'NKG7', 'TREM1', 'LYZ'),
pt.size = 0.1,
max.cutoff = 'q95',
ncol = 3
)
Integrating with scRNA-seq data
To help interpret the scATAC-seq data, we can classify cells based on an scRNA-seq experiment from the same biological system (human PBMC). We utilize methods for cross-modality integration and label transfer, described here, with a more in-depth tutorial here. We aim to identify shared correlation patterns in the gene activity matrix and scRNA-seq dataset to identify matched biological states across the two modalities. This procedure returns a classification score for each cell for each scRNA-seq-defined cluster label.
Here we load a pre-processed scRNA-seq dataset for human PBMCs, also provided by 10x Genomics. You can download the raw data for this experiment from the 10x website, and view the code used to construct this object on GitHub. Alternatively, you can download the pre-processed Seurat object here.
# Load the pre-processed scRNA-seq data for PBMCs
pbmc_rna <- readRDS("/home/stuartt/github/chrom/vignette_data/pbmc_10k_v3.rds")
transfer.anchors <- FindTransferAnchors(
reference = pbmc_rna,
query = pbmc,
reduction = 'cca'
)
predicted.labels <- TransferData(
anchorset = transfer.anchors,
refdata = pbmc_rna$celltype,
weight.reduction = pbmc[['lsi']],
dims = 2:30
)
pbmc <- AddMetaData(object = pbmc, metadata = predicted.labels)
plot1 <- DimPlot(
object = pbmc_rna,
group.by = 'celltype',
label = TRUE,
repel = TRUE) + NoLegend() + ggtitle('scRNA-seq')
plot2 <- DimPlot(
object = pbmc,
group.by = 'predicted.id',
label = TRUE,
repel = TRUE) + NoLegend() + ggtitle('scATAC-seq')
plot1 + plot2
You can see that the scRNA-based classifications are entirely consistent with the UMAP visualization, computed only on the scATAC-seq data. We can now easily annotate our scATAC-seq-derived clusters (alternatively, we could use the RNA classifications themselves). We note that cluster 14 maps to CD4 Memory T cells, but is a very small cluster with lower QC metrics. As this group is likely representing low-quality cells, we remove it from downstream analysis.
pbmc <- subset(pbmc, idents = 14, invert = TRUE)
pbmc <- RenameIdents(
object = pbmc,
'0' = 'CD14 Mono',
'1' = 'CD4 Memory',
'2' = 'CD8 Effector',
'3' = 'CD4 Naive',
'4' = 'CD14 Mono',
'5' = 'DN T',
'6' = 'CD8 Naive',
'7' = 'NK CD56Dim',
'8' = 'pre-B',
'9' = 'CD16 Mono',
'10' = 'pro-B',
'11' = 'DC',
'12' = 'NK CD56bright',
'13' = 'pDC'
)Find differentially accessible peaks between clusters
To find differentially accessible regions between clusters of cells, we can perform a differential accessibility (DA) test. We utilize logistic regression for DA, as suggested by Ntranos et al. 2018 for scRNA-seq data, and add the total number of fragments as a latent variable to mitigate the effect of differential sequencing depth on the result. Here we will focus on comparing Naive CD4 cells and CD14 monocytes, but any groups of cells can be compared using these methods. We can also visualize these marker peaks on a violin plot, feature plot, dot plot, heat map, or any visualization tool in Seurat.
# change back to working with peaks instead of gene activities
DefaultAssay(pbmc) <- 'peaks'
da_peaks <- FindMarkers(
object = pbmc,
ident.1 = "CD4 Naive",
ident.2 = "CD14 Mono",
min.pct = 0.2,
test.use = 'LR',
latent.vars = 'peak_region_fragments'
)
head(da_peaks)## p_val avg_log2FC pct.1 pct.2 p_val_adj
## chr14-99721608-99741934 1.622701e-307 1.1465091 0.879 0.023 1.420853e-302
## chr14-99695477-99720910 3.713144e-241 1.1351216 0.825 0.023 3.251266e-236
## chr17-80084198-80086094 1.789358e-231 1.6124656 0.695 0.006 1.566780e-226
## chr7-142501666-142511108 1.284659e-230 0.9959048 0.767 0.032 1.124860e-225
## chr2-113581628-113594911 2.820671e-212 -1.2820550 0.055 0.693 2.469808e-207
## chr6-44025105-44028184 3.989239e-206 -1.2508915 0.063 0.654 3.493017e-201
plot1 <- VlnPlot(
object = pbmc,
features = rownames(da_peaks)[1],
pt.size = 0.1,
idents = c("CD4 Memory","CD14 Mono")
)
plot2 <- FeaturePlot(
object = pbmc,
features = rownames(da_peaks)[1],
pt.size = 0.1
)
plot1 | plot2
Another way to find DA regions between two groups of cells is to look at the fold change accessibility between two groups of cells. This can be much faster than running more sophisticated DA tests, but is not able to account for latent variables such as differences in the total sequencing depth between cells, and does not perform any statistical test. However, this can still be a useful way to quickly explore data, and can be performed using the FoldChange() function in Signac.
fc <- FoldChange(pbmc, ident.1 = "CD4 Naive", ident.2 = "CD14 Mono")
head(fc)## avg_log2FC pct.1 pct.2
## chr1-565107-565550 0.1947304 0.012 0.004
## chr1-569174-569639 0.3883631 0.032 0.008
## chr1-713460-714823 0.1935680 0.421 0.208
## chr1-752422-753038 -1.1668007 0.010 0.094
## chr1-762106-763359 0.1033536 0.290 0.162
## chr1-779589-780271 1.2677633 0.051 0.002Peak coordinates can be difficult to interpret alone. We can find the closest gene to each of these peaks using the ClosestFeature() function. If you explore the gene lists, you will see that peaks open in Naive T cells are close to genes such as BCL11B and GATA3 (key regulators of T cell differentiation ), while peaks open in monocytes are close to genes such as CEBPB (a key regulator of monocyte differentiation). We could follow up this result further by doing gene ontology enrichment analysis on the gene sets returned by ClosestFeature(), and there are many R packages that can do this (see the GOstats package for example).
open_cd4naive <- rownames(da_peaks[da_peaks$avg_log2FC > 0.5, ])
open_cd14mono <- rownames(da_peaks[da_peaks$avg_log2FC < -0.5, ])
closest_genes_cd4naive <- ClosestFeature(pbmc, regions = open_cd4naive)
closest_genes_cd14mono <- ClosestFeature(pbmc, regions = open_cd14mono)
head(closest_genes_cd4naive)## tx_id gene_name gene_id gene_biotype type
## ENST00000443726 ENST00000443726 BCL11B ENSG00000127152 protein_coding cds
## ENST00000357195 ENST00000357195 BCL11B ENSG00000127152 protein_coding cds
## ENST00000583593 ENST00000583593 CCDC57 ENSG00000176155 protein_coding cds
## ENSE00002456092 ENST00000463701 PRSS1 ENSG00000204983 protein_coding exon
## ENST00000455990 ENST00000455990 HOOK1 ENSG00000134709 protein_coding cds
## ENST00000557595 ENST00000557595 AE000662.92 ENSG00000259003 protein_coding cds
## closest_region query_region distance
## ENST00000443726 chr14-99737498-99737555 chr14-99721608-99741934 0
## ENST00000357195 chr14-99697682-99697894 chr14-99695477-99720910 0
## ENST00000583593 chr17-80085568-80085694 chr17-80084198-80086094 0
## ENSE00002456092 chr7-142460719-142460923 chr7-142501666-142511108 40742
## ENST00000455990 chr1-60280790-60280852 chr1-60279767-60281364 0
## ENST00000557595 chr14-23027629-23027664 chr14-23012859-23028219 0
head(closest_genes_cd14mono)## tx_id gene_name gene_id gene_biotype
## ENST00000432018 ENST00000432018 IL1B ENSG00000125538 protein_coding
## ENSE00001638912 ENST00000455005 RP5-1120P11.3 ENSG00000231881 lincRNA
## ENST00000445003 ENST00000445003 RP11-290F20.3 ENSG00000224397 lincRNA
## ENST00000568649 ENST00000568649 PPCDC ENSG00000138621 protein_coding
## ENST00000484822 ENST00000484822 RXRA ENSG00000186350 protein_coding
## ENST00000560869 ENST00000560869 TNFAIP2 ENSG00000185215 protein_coding
## type closest_region query_region
## ENST00000432018 cds chr2-113593760-113593806 chr2-113581628-113594911
## ENSE00001638912 exon chr6-44041650-44042535 chr6-44025105-44028184
## ENST00000445003 gap chr20-48884201-48894027 chr20-48889794-48893313
## ENST00000568649 cds chr15-75335782-75335877 chr15-75334903-75336779
## ENST00000484822 gap chr9-137211331-137293477 chr9-137263243-137268534
## ENST00000560869 utr chr14-103589798-103590288 chr14-103587259-103591113
## distance
## ENST00000432018 0
## ENSE00001638912 13465
## ENST00000445003 0
## ENST00000568649 0
## ENST00000484822 0
## ENST00000560869 0Plotting genomic regions
We can plot the frequency of Tn5 integration across regions of the genome for cells grouped by cluster, cell type, or any other metadata stored in the object for any genomic region using the CoveragePlot() function. These represent pseudo-bulk accessibility tracks, where signal from all cells within a group have been averaged together to visualize the DNA accessibility in a region (thanks to Andrew Hill for giving the inspiration for this function in his excellent blog post). Alongside these accessibility tracks we can visualize other important information including gene annotation, peak coordinates, and genomic links (if they’re stored in the object). See the visualization vignette for more information.
# set plotting order
levels(pbmc) <- c("CD4 Naive","CD4 Memory","CD8 Naive","CD8 Effector","DN T","NK CD56bright","NK CD56Dim","pre-B",'pro-B',"pDC","DC","CD14 Mono",'CD16 Mono')
CoveragePlot(
object = pbmc,
region = rownames(da_peaks)[1],
extend.upstream = 40000,
extend.downstream = 20000
)
We can also create an interactive version of these plots using the CoverageBrowser() function. Here is a recorded demonstration showing how we can use CoverageBrowser() to browser the genome and adjust plotting parameters interactively. The “Save plot” button in CoverageBrowser() will add the current plot to a list of ggplot objects that is returned when the browser session is ended by pressing the “Done” button, allowing interesting views to be saved during an interactive session.
Working with datasets that were not quantified using CellRanger
The CellRanger software from 10x Genomics generates several useful QC metrics per-cell, as well as a peak/cell matrix and an indexed fragments file. In the above vignette, we utilize the CellRanger outputs, but provide alternative functions in Signac for many of the same purposes here.
Generating a peak/cell or bin/cell matrix
The FeatureMatrix function can be used to generate a count matrix containing any set of genomic ranges in its rows. These regions could be a set of peaks, or bins that span the entire genome.
# not run
# peak_ranges should be a set of genomic ranges spanning the set of peaks to be quantified per cell
peak_matrix <- FeatureMatrix(
fragments = Fragments(pbmc),
features = peak_ranges
)For convenience, we also include a GenomeBinMatrix() function that will generate a set of genomic ranges spanning the entire genome for you, and run FeatureMatrix() internally to produce a genome bin/cell matrix.
# not run
bin_matrix <- GenomeBinMatrix(
fragments = Fragments(pbmc),
genome = seqlengths(pbmc),
binsize = 5000
)Counting fraction of reads in peaks
The function FRiP() will count the fraction of reads in peaks for each cell, given a peak/cell assay and a bin/cell assay. Note that this can be run on a subset of the genome, so that a bin/cell assay does not need to be computed for the whole genome. This will return a Seurat object will metadata added corresponding to the fraction of reads in peaks for each cell.
# not run
total_fragments <- CountFragments("'/home/stuartt/github/chrom/vignette_data/atac_v1_pbmc_10k_fragments.tsv.gz'")
pbmc$fragments <- total_fragments[colnames(pbmc), "frequency_count"]
pbmc <- FRiP(
object = pbmc,
assay = 'peaks',
total.fragments = 'fragments'
)Counting fragments in genome blacklist regions
The ratio of reads in genomic blacklist regions, that are known to artifactually accumulate reads in genome sequencing assays, can be diagnostic of low-quality cells. We provide blacklist region coordinates for several genomes (hg19, hg38, mm9, mm10, ce10, ce11, dm3, dm6) in the Signac package for convenience. These regions were provided by the ENCODE consortium, and we encourage users to cite their paper if you use the regions in your analysis. The FractionCountsInRegion() function can be used to calculate the fraction of all counts within a given set of regions per cell. We can use this function and the blacklist regions to find the fraction of blacklist counts per cell.
# not run
pbmc$blacklist_fraction <- FractionCountsInRegion(
object = pbmc,
assay = 'peaks',
regions = blacklist_hg19
)Session Info
## R version 4.0.1 (2020-06-06)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats4 parallel stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] patchwork_1.1.1 ggplot2_3.3.3
## [3] EnsDb.Hsapiens.v75_2.99.0 ensembldb_2.14.0
## [5] AnnotationFilter_1.14.0 GenomicFeatures_1.42.3
## [7] AnnotationDbi_1.52.0 Biobase_2.50.0
## [9] GenomicRanges_1.42.0 GenomeInfoDb_1.26.5
## [11] IRanges_2.24.1 S4Vectors_0.28.1
## [13] BiocGenerics_0.36.0 SeuratObject_4.0.0
## [15] Seurat_4.0.1.9005 Signac_1.2.0
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.1 reticulate_1.19
## [3] tidyselect_1.1.0 RSQLite_2.2.7
## [5] htmlwidgets_1.5.3 grid_4.0.1
## [7] docopt_0.7.1 BiocParallel_1.24.1
## [9] Rtsne_0.15 munsell_0.5.0
## [11] codetools_0.2-18 ragg_1.1.2
## [13] ica_1.0-2 future_1.21.0
## [15] miniUI_0.1.1.1 withr_2.4.2
## [17] colorspace_2.0-0 highr_0.9
## [19] knitr_1.33 rstudioapi_0.13
## [21] ROCR_1.0-11 tensor_1.5
## [23] listenv_0.8.0 labeling_0.4.2
## [25] MatrixGenerics_1.2.1 slam_0.1-48
## [27] GenomeInfoDbData_1.2.4 polyclip_1.10-0
## [29] bit64_4.0.5 farver_2.1.0
## [31] rprojroot_2.0.2 parallelly_1.24.0
## [33] vctrs_0.3.7 generics_0.1.0
## [35] xfun_0.22 biovizBase_1.38.0
## [37] BiocFileCache_1.14.0 lsa_0.73.2
## [39] ggseqlogo_0.1 R6_2.5.0
## [41] hdf5r_1.3.3 bitops_1.0-7
## [43] spatstat.utils_2.1-0 cachem_1.0.4
## [45] DelayedArray_0.16.3 assertthat_0.2.1
## [47] promises_1.2.0.1 scales_1.1.1
## [49] nnet_7.3-15 gtable_0.3.0
## [51] globals_0.14.0 goftest_1.2-2
## [53] rlang_0.4.10 systemfonts_1.0.1
## [55] RcppRoll_0.3.0 splines_4.0.1
## [57] rtracklayer_1.50.0 lazyeval_0.2.2
## [59] dichromat_2.0-0 checkmate_2.0.0
## [61] spatstat.geom_2.1-0 yaml_2.2.1
## [63] reshape2_1.4.4 abind_1.4-5
## [65] backports_1.2.1 httpuv_1.6.0
## [67] Hmisc_4.5-0 tools_4.0.1
## [69] ellipsis_0.3.1 spatstat.core_2.1-2
## [71] jquerylib_0.1.4 RColorBrewer_1.1-2
## [73] ggridges_0.5.3 Rcpp_1.0.6
## [75] plyr_1.8.6 base64enc_0.1-3
## [77] progress_1.2.2 zlibbioc_1.36.0
## [79] purrr_0.3.4 RCurl_1.98-1.3
## [81] prettyunits_1.1.1 rpart_4.1-15
## [83] openssl_1.4.3 deldir_0.2-10
## [85] pbapply_1.4-3 cowplot_1.1.1
## [87] zoo_1.8-9 SummarizedExperiment_1.20.0
## [89] ggrepel_0.9.1 cluster_2.1.2
## [91] fs_1.5.0 magrittr_2.0.1
## [93] RSpectra_0.16-0 data.table_1.14.0
## [95] scattermore_0.7 lmtest_0.9-38
## [97] RANN_2.6.1 SnowballC_0.7.0
## [99] ProtGenerics_1.22.0 fitdistrplus_1.1-3
## [101] matrixStats_0.58.0 hms_1.0.0
## [103] mime_0.10 evaluate_0.14
## [105] xtable_1.8-4 XML_3.99-0.6
## [107] jpeg_0.1-8.1 sparsesvd_0.2
## [109] gridExtra_2.3 compiler_4.0.1
## [111] biomaRt_2.46.3 tibble_3.1.1
## [113] KernSmooth_2.23-18 crayon_1.4.1
## [115] htmltools_0.5.1.1 mgcv_1.8-33
## [117] later_1.2.0 Formula_1.2-4
## [119] tidyr_1.1.3 DBI_1.1.1
## [121] tweenr_1.0.2 dbplyr_2.1.1
## [123] MASS_7.3-53.1 rappdirs_0.3.3
## [125] Matrix_1.3-2 igraph_1.2.6
## [127] pkgconfig_2.0.3 pkgdown_1.6.1
## [129] GenomicAlignments_1.26.0 foreign_0.8-81
## [131] plotly_4.9.3 spatstat.sparse_2.0-0
## [133] xml2_1.3.2 bslib_0.2.4
## [135] XVector_0.30.0 VariantAnnotation_1.36.0
## [137] stringr_1.4.0 digest_0.6.27
## [139] sctransform_0.3.2 RcppAnnoy_0.0.18
## [141] spatstat.data_2.1-0 Biostrings_2.58.0
## [143] rmarkdown_2.7 leiden_0.3.7
## [145] fastmatch_1.1-0 htmlTable_2.1.0
## [147] uwot_0.1.10 curl_4.3
## [149] shiny_1.6.0 Rsamtools_2.6.0
## [151] lifecycle_1.0.0 nlme_3.1-152
## [153] jsonlite_1.7.2 BSgenome_1.58.0
## [155] desc_1.3.0 viridisLite_0.4.0
## [157] askpass_1.1 fansi_0.4.2
## [159] pillar_1.6.0 lattice_0.20-41
## [161] fastmap_1.1.0 httr_1.4.2
## [163] survival_3.2-11 glue_1.4.2
## [165] qlcMatrix_0.9.7 png_0.1-7
## [167] bit_4.0.4 ggforce_0.3.3
## [169] stringi_1.5.3 sass_0.3.1
## [171] blob_1.2.1 textshaping_0.3.3
## [173] latticeExtra_0.6-29 memoise_2.0.0
## [175] dplyr_1.0.5 irlba_2.3.3
## [177] future.apply_1.7.0