This section introduces how to process Cell Ranger results into an eccDNA object (Seurat object), and complete the processes of dimensionality reduction, clustering, gene activity calculation, and cell annotation. The example dataset combined_atac can be downloaded here, the final combined_atac_annote file can be downloaded here.

Convert Cell Ranger results into a eccDNA(Seurat) object

Each sample needs to include four input files

## eccDNAscope_example/sample_data
## ├── AK1
## │   └── result_cellranger-atac
## │       ├── AK1_fragments.tsv.gz
## │       ├── AK1_fragments.tsv.gz.tbi
## │       ├── peaks.bed
## │       └── singlecell.csv
## └── CSCC1
##     └── result_cellranger-atac
##         ├── CSCC1_fragments.tsv.gz
##         ├── CSCC1_fragments.tsv.gz.tbi
##         ├── peaks.bed
##         └── singlecell.csv

Run the create_eccDNA_obj function

library(eccDNAscope)
## 
##**************************************************************************
##  Welcome to eccDNAscope!
##      Version: 0.1.0
##                   _____________   ________                                 
##____________________  __ \__  | / /__    |_______________________________ 
##_  _ \  ___/  ___/_  / / /_   |/ /__  /| |_  ___/  ___/  __ \__  __ \  _ \
##/  __/ /__ / /__ _  /_/ /_  /|  / _  ___ |(__  )/ /__ / /_/ /_  /_/ /  __/
##\___/\___/ \___/ /_____/ /_/ |_/  /_/  |_/____/ \___/ \____/_  .___/\___/ 
##                                                            /_/    														
##
##**************************************************************************
## 
combined_atac <- create_eccDNA_obj(dir = "eccDNAscope_example/sample_data",
                                               cellranger.result.dir="result_cellranger-atac")
## 1. Reading the peaks.bed file
## 2. Converting the BED file into genomic ranges
## 3. Creating a set of unified peak values to quantify each dataset
## 4. Filtering peaks based on their width
## 5. Reading the singlecell.csv file
## 6. Filtering out cells with low counts
## 7. Generating fragment objects
## 8. Computing the FeatureMatrix
## 9. Creating Chromatin Assay object
## 10. Creating Seurat object
## 11. Adding metadata
## 12. Merging Seurat objects

View the combined_atac object

combined_atac
## An object of class Seurat 
## 96663 features across 140862 samples within 1 assay 
## Active assay: ATAC (96663 features, 0 variable features)
##  2 layers present: counts, data


Add genome annotation information to the eccDNA object

library(EnsDb.Hsapiens.v86)  # Human hg38 2.99.0
library(Signac)  # 1.14.0

# extract gene annotations from EnsDb
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)

# change to UCSC style since the data was mapped to hg38
seqlevels(annotations) <- paste0('chr', seqlevels(annotations))
genome(annotations) <- "hg38"

# add the gene information to the object
Annotation(combined_atac) <- annotations


Perform quality control on the combined_atac

combined_atac <- perform_QC(combined_atac)
## 1. Calculating NucleosomeSignal
## 2. Calculating TSS Enrichment score
## 3. Calculating the proportion of DNA fragments within peaks for each cell
## 4. Calculating the proportion of DNA fragments within genomic blacklist regions for each cell
## 5. Plotting a density scatter plot between n_count_ATAC and TSS.enrichment before QC
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## 
## 6. Plotting a TSSPlot to compare cells with high and low TSS scores
## 7. Plotting a FragmentHistogram to compare cells with high and low nucleosome_signal
## 8. Plot violin plots for various data quality metrics before QC
## 
## 9. performing quality control
## 10. Plotting a density scatter plot between n_count_ATAC and TSS.enrichment after QC
## Scale for colour is already present.
## Adding another scale for colour, which will replace the existing scale.
## 
## 11. Plot violin plots for various data quality metrics after QC

Display changes in the data before and after quality control

# DensityScatter
combined_atac@misc$p_DensityScatter_beforeQC

combined_atac@misc$p_DensityScatter_afterQC

# VlnPlot
combined_atac@misc$p_data_beforeQC

combined_atac@misc$p_data_afterQC

# TSSPlot
combined_atac@misc$p_TSSPlot

# FragmentHistogram
combined_atac@misc$p_FragmentHistogram


Perform dimensionality reduction, clustering, and cell annotation of the eccDNA object using the Signac R package

dimensionality reduction

library(Signac)
library(Seurat)
combined_atac <- RunTFIDF(combined_atac)
combined_atac <- FindTopFeatures(combined_atac, min.cutoff = 10)
combined_atac <- RunSVD(combined_atac)
DepthCor(combined_atac)

combined_atac <- RunUMAP(combined_atac, dims = 2:30, reduction = 'lsi')
DimPlot(combined_atac, group.by = 'dataset', label.size = 4 ,pt.size = 0.5, raster=FALSE)

clustering

DefaultAssay(combined_atac) <- 'ATAC'
combined_atac <- FindNeighbors(object = combined_atac, reduction = 'lsi', dims = 2:30)
combined_atac <- FindClusters(object = combined_atac, verbose = F, algorithm = 3)
DimPlot(object = combined_atac, group.by = 'seurat_clusters', label = TRUE,label.size = 4 ,pt.size = 0.5, raster=FALSE)

Generate the gene activity matrix

gene.activities <- GeneActivity(combined_atac)
combined_atac[['RNA']] <- CreateAssayObject(counts = gene.activities)
combined_atac <- NormalizeData(
  object = combined_atac,
  assay = 'RNA',
  normalization.method = 'LogNormalize',
  scale.factor = median(combined_atac$nCount_RNA),
  verbose = F
)

cell annotation

DefaultAssay(combined_atac) <- 'RNA'
combined_atac_final <- readRDS("eccDNAscope_example/clone_data/combined_atac_final.rds")
combined_atac <- combined_atac[,intersect(colnames(combined_atac),colnames(combined_atac_final))]
combined_atac$celltype_ATAC <- combined_atac_final$celltype_ATAC[match(colnames(combined_atac),colnames(combined_atac_final))]
DimPlot(combined_atac, group.by = 'celltype_ATAC', label = TRUE,label.size = 4 ,pt.size = 0.5, raster=FALSE, repel = TRUE)
saveRDS(combined_atac , file = "eccDNAscope_example/result_merge/combined_atac_annote.rds")