scPipe: a flexible data preprocessing pipeline for scATAC-seq data
Introduction
scPipe is a package initially designed to process
single-cell RNA-sequencing (scRNA-seq) data generated by different
protocols. We have modified it to accommodate pre-processing capability
of single-cell ATAC-Seq (Assay for Transposase-Accessible Chromatin
using sequencing) data pre-processing. scPipe ATAC-Seq
module is designed for protocols without UMIs, but can also adapt to any
UMI protocols if the need arise.
scPipe ATAC-Seq module consists of two major components.
The first is data pre-processing with barcodes as raw fastq as input and
a feature count matrix as output. Second is the data pre-processing with
barcode CSV file as input and a feature count matrix as output.
10X ATAC method currently is the most popular method to
generate scATAC-Seq data with higher sensitivity and lower cost. The
structure of the 10X ATAC library is shown below.
The output fastq files from a 10X ATAC experiment is
paired-ended and data is contained within both reads.
Workflow
The pipeline is visually described via the workflow diagram depicted below. The functions will be explained in greater detail.
Getting started
It is not mandatory to specify an output folder even though it can be
specified. If no output_folder is defined a folder named
scPipe-atac-output will get created in the working
directory.
We begin by loading the library.
Specify the output folder.
Fastq reformatting
To process the data, we need the fastq files (both
compressed and uncompressed versions are accepted) and a cell barcode
annotation. We have supplied very small sample FASTQ files of chr21. The
barcode annotation could be either in a .fastq format or a
.csv file with at least two columns, where the first column
has the cell id and the second column contains the barcode sequence. All
these files can be found in the data folder of the
scPipe package:
# data fastq files
r1 <- file.path(data.folder, "small_chr21_R1.fastq.gz")
r2 <- file.path(data.folder, "small_chr21_R3.fastq.gz")
# barcodes in fastq format
barcode_fastq <- file.path(data.folder, "small_chr21_R2.fastq.gz")
# barcodes in .csv format
barcode_1000 <- file.path(data.folder, "chr21_modified_barcode_1000.csv")The pipeline starts with fastq file reformatting. We move the barcode
and UMI sequences (if available) to the read name and leave the
transcript sequence as is. This outputs a read name that looks like
@[barcode_sequence]_[UMI_sequence]#[readname]. Usually the
scATAC-Seq data is paired-end and a 16bp long barcode is located on both
reads. Here the barcode information is located on a separate fastq file
and the length of the barcode fastq file matches the length of the reads
files. Therefore, you need a minimal example like below to generate the
output.
sc_atac_trim_barcode (r1 = r1,
r2 = r2,
bc_file = barcode_fastq,
rmN = FALSE,
rmlow = FALSE,
output_folder = output_folder)## Output Directory Does Not Exist. Created Directory: /tmp/RtmpPhWjmQ/scPipeATACVignette## Saving the output at location:
## /tmp/RtmpPhWjmQ/scPipeATACVignette## No valid_barcode_file provided; no barcode error correction will occur.
## Total reads: 35000
## Total reads removed due to N's in barcodes: 0
## Total reads removed due to low quality barcodes: 0
## Total barcodes provided in FASTQ file: 16567
## time elapsed: 194 millisecondsThis generates two output fastq files that are appended by the prefix
demux_ to notify that the new files are the reformatted
(a.k.a. demultiplexed) .fastq files. These will get saved
in the scPipe-atac-output directory if the user has not
specified an output_folder.
However, if the barcodes are in the form of a .csv file,
some extra information on 0-indexed barcode start (id1_st,
id2_st) and the barcode length (id1_len,
id2_len) are also required to be entered by the user. The
algorithm is flexible to do a “look-around” to identify whether the
correct parameters are used for a subset of data (hence saving time) and
report back to the console if it believes the barcode position is
incorrect and/or should be shifted.
sc_atac_trim_barcode (r1 = r1,
r2 = r2,
bc_file = barcode_1000,
id1_st = -1,
id1_len = -1,
id2_st = -1,
id2_len = -10,
output_folder = output_folder,
rmN = TRUE)To accommodate combinatorial indexing in some scATAC-seq protocols,
the bc_file parameter will accept a list of barcode files
as well (currently only implemented for the fastq
approach).
Additionally, rmN, rmlow and
min_qual parameters define the quality thresholds for the
reads. If there are Ns in the barcode or UMI positions
those will be discarded by rmN = TRUE.
rmlow =TRUE will remove reads having lower quality than
what is defined by min_qual (default: 20).
Completion of this function will output three different outputs depending on the findings; * complete matches: When the barcode is completely matched and identified in the correct position * partial matches: When the barcode is identified in the location specified but corrected with hamming distance approach * unmatched: no barcode match is found in the given position even after hamming distance corrections are applied
NOTE: we use a zero based index system, so the indexing of the sequence starts at zero.
Aligning reads to a reference genome
Next, we align reads to the genome. This example uses
Rsubread but any aligner that support RNA-seq alignment and
gives standard BAM output can be used here.
demux_r1 <- file.path(output_folder, "demux_completematch_small_chr21_R1.fastq.gz")
demux_r2 <- file.path(output_folder, "demux_completematch_small_chr21_R3.fastq.gz")
reference = file.path(data.folder, "small_chr21.fa")
aligned_bam <- sc_aligning(ref = reference,
R1 = demux_r1,
R2 = demux_r2,
nthreads = 12,
output_folder = output_folder)## ATAC-Seq mode is selected...## Genome index location not specified. Looking for the index in/tmp/RtmpPhWjmQ/scPipeATACVignette## Genome index not found. Creating one in /tmp/RtmpPhWjmQ/scPipeATACVignette ...##
## ========== _____ _ _ ____ _____ ______ _____
## ===== / ____| | | | _ \| __ \| ____| /\ | __ \
## ===== | (___ | | | | |_) | |__) | |__ / \ | | | |
## ==== \___ \| | | | _ <| _ /| __| / /\ \ | | | |
## ==== ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
## ========== |_____/ \____/|____/|_| \_\______/_/ \_\_____/
## Rsubread 2.21.1
##
## //================================= setting ==================================\\
## || ||
## || Index name : genome_index ||
## || Index space : base space ||
## || Index split : no-split ||
## || Repeat threshold : 100 repeats ||
## || Gapped index : no ||
## || ||
## || Free / total memory : 13.3GB / 15.6GB ||
## || ||
## || Input files : 1 file in total ||
## || o small_chr21.fa ||
## || ||
## \\============================================================================//
##
## //================================= Running ==================================\\
## || ||
## || Check the integrity of provided reference sequences ... ||
## || No format issues were found ||
## || Scan uninformative subreads in reference sequences ... ||
## || 3 uninformative subreads were found. ||
## || These subreads were excluded from index building. ||
## || Estimate the index size... ||
## || 8%, 0 mins elapsed, rate=1640.8k bps/s ||
## || 16%, 0 mins elapsed, rate=2788.4k bps/s ||
## || 24%, 0 mins elapsed, rate=3874.0k bps/s ||
## || 33%, 0 mins elapsed, rate=4556.1k bps/s ||
## || 41%, 0 mins elapsed, rate=5342.3k bps/s ||
## || 49%, 0 mins elapsed, rate=5788.5k bps/s ||
## || 58%, 0 mins elapsed, rate=6127.8k bps/s ||
## || 66%, 0 mins elapsed, rate=6674.5k bps/s ||
## || 74%, 0 mins elapsed, rate=6896.3k bps/s ||
## || 83%, 0 mins elapsed, rate=7324.6k bps/s ||
## || 91%, 0 mins elapsed, rate=7946.4k bps/s ||
## || 3.0 GB of memory is needed for index building. ||
## || Build the index... ||
## || 8%, 0 mins elapsed, rate=159.7k bps/s ||
## || 16%, 0 mins elapsed, rate=308.0k bps/s ||
## || 24%, 0 mins elapsed, rate=453.2k bps/s ||
## || 33%, 0 mins elapsed, rate=586.4k bps/s ||
## || 41%, 0 mins elapsed, rate=720.2k bps/s ||
## || 49%, 0 mins elapsed, rate=840.3k bps/s ||
## || 58%, 0 mins elapsed, rate=954.3k bps/s ||
## || 66%, 0 mins elapsed, rate=1073.6k bps/s ||
## || 74%, 0 mins elapsed, rate=1176.4k bps/s ||
## || 83%, 0 mins elapsed, rate=1287.1k bps/s ||
## || 91%, 0 mins elapsed, rate=1405.1k bps/s ||
## || Save current index block... ||
## || [ 0.0% finished ] ||
## || [ 10.0% finished ] ||
## || [ 20.0% finished ] ||
## || [ 30.0% finished ] ||
## || [ 40.0% finished ] ||
## || [ 50.0% finished ] ||
## || [ 60.0% finished ] ||
## || [ 70.0% finished ] ||
## || [ 80.0% finished ] ||
## || [ 90.0% finished ] ||
## || [ 100.0% finished ] ||
## || ||
## || Total running time: 0.1 minutes. ||
## |Index /tmp/RtmpPhWjmQ/scPipeATACVignette/genome_index was successfully built. ||
## || ||
## \\============================================================================//## No partial matches, checking for reads with non-matched barcodes.## No reads found with non-matching barcodes.## Outputted demultiplexing stats file to/tmp/RtmpPhWjmQ/scPipeATACVignette/scPipe_atac_stats//demultiplexing_stats.csv## Output file name is not provided. Aligned reads are saved in /tmp/RtmpPhWjmQ/scPipeATACVignette/demux_completematch_small_chr21_R1_aligned.bam##
## ========== _____ _ _ ____ _____ ______ _____
## ===== / ____| | | | _ \| __ \| ____| /\ | __ \
## ===== | (___ | | | | |_) | |__) | |__ / \ | | | |
## ==== \___ \| | | | _ <| _ /| __| / /\ \ | | | |
## ==== ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
## ========== |_____/ \____/|____/|_| \_\______/_/ \_\_____/
## Rsubread 2.21.1
##
## //================================= setting ==================================\\
## || ||
## || Function : Read alignment (DNA-Seq) ||
## || Input file 1 : demux_completematch_small_chr21_R1.fastq.gz ||
## || Input file 2 : demux_completematch_small_chr21_R3.fastq.gz ||
## || Output file : demux_completematch_small_chr21_R1_aligned.bam (BAM), ... ||
## || Index name : genome_index ||
## || ||
## || ------------------------------------ ||
## || ||
## || Threads : 12 ||
## || Phred offset : 33 ||
## || # of extracted subreads : 10 ||
## || Min read1 vote : 3 ||
## || Min read2 vote : 1 ||
## || Max fragment size : 600 ||
## || Min fragment size : 50 ||
## || Max mismatches : 3 ||
## || Max indel length : 5 ||
## || Report multi-mapping reads : yes ||
## || Max alignments per multi-mapping read : 1 ||
## || ||
## \\============================================================================//
##
## //================ Running (18-Dec-2024 08:39:03, pid=12833) =================\\
## || ||
## || Check the input reads. ||
## || The input file contains base space reads. ||
## || Initialise the memory objects. ||
## || Estimate the mean read length. ||
## || The range of Phred scores observed in the data is [2,37] ||
## || Create the output BAM file. ||
## || Check the index. ||
## || Init the voting space. ||
## || Global environment is initialised. ||
## || Load the 1-th index block... ||
## || The index block has been loaded. ||
## || Start read mapping in chunk. ||
## || Estimated fragment length : 123 bp ||
## || ||
## || Completed successfully. ||
## || ||
## \\==================================== ====================================//
##
## //================================ Summary =================================\\
## || ||
## || Total fragments : 35,000 ||
## || Mapped : 18,851 (53.9%) ||
## || Uniquely mapped : 18,057 ||
## || Multi-mapping : 794 ||
## || ||
## || Unmapped : 16,149 ||
## || ||
## || Properly paired : 14,524 ||
## || Not properly paired : 4,327 ||
## || Singleton : 1,458 ||
## || Chimeric : 0 ||
## || Unexpected strandness : 114 ||
## || Unexpected fragment length : 180 ||
## || Unexpected read order : 2,575 ||
## || ||
## || Indels : 197 ||
## || ||
## || Running time : 1.0 minutes ||
## || ||
## \\============================================================================//Demultiplexing the BAM file
Next, the BAM file needs to be modified in a way one or two new
columns are generated for the cellular barcode tag and the molecular
barcode (i.e. UMI) tag denoted by CB:Z: and
OX:Z:, respectively.
sorted_tagged_bam <- sc_atac_bam_tagging (inbam = aligned_bam,
output_folder = output_folder,
bam_tags = list(bc="CB", mb="OX"),
nthreads = 12)## Detected bc_len: 16 Detected UMI len: 0
## updating progress every 3 minutes...
## number of read processed:0
## 0 0 0 0
## 70000 reads processed, 70k reads/sec
## number of read processed: 70000
## time elapsed: 1023 milliseconds## Demultiplexed BAM file sorted and indexed ...## Using default value for barcode length (bc_length = 16)## Generating mapping statistics per barcode## Iterating through 5,000,000 reads at a time## chunk1## Completed generating mapping statistics per barcode, saved to /tmp/RtmpPhWjmQ/scPipeATACVignette/scPipe_atac_stats/mapping_stats_per_barcode.csv## The output tagged and sorted BAM file was sent to /tmp/RtmpPhWjmQ/scPipeATACVignetteRemove duplicates
Next, PCR duplicate reads are removed from the BAM file. If
samtools is installed, then
sc_atac_remove_duplicates can be used, with the
installation path of samtools being an optional argument if
a particular version is preferred. However, if samtools
isn’t installed, the PCR duplicate removal should be performed
externally.
Gemerating a fragment file
Next, a fragment file in a .bed format is
generated, where each line represents a unique fragment generated by the
assay. This file is used to generate useful quality control statistics,
as well as for the filter and cellranger cell
calling methods.
## Output folder name is: /tmp/RtmpPhWjmQ/scPipeATACVignette
## Output BED file: /tmp/RtmpPhWjmQ/scPipeATACVignette/fragments.bedPeak calling
Similar to many other tools, the the ability to call peaks on a
pseudo-bulk level on the demultiplexed reads has also been incorporated.
MACS3 is used underneath to achieve this functionality. If
the genome size is not provided then it will be roughly estimated.
Assigning reads to features and feature counting
After the read alignment and BAM file demultiplexing, a feature set
can be used to find the overlap between the aligned reads and the
features using the sc_atac_feature_counting function.
Accepted format of the feature_input should be either a
BED format (i.e. format should have three main columns; chromosome,
feature start and feature end, extension of the file is not considered)
or a genome.fasta file. If using a BED format as the
feature_input, the feature_type should be
either “peak” or “tss”. If using a .fasta for the
feature_input, the feature_type needs to be
defined as genome_bin. This genome.fasta file
will be used within the function to generate a genome_bin
that defines the array of features. The size of the bins can be set
using the bin_size parameter (default: 2000).
Note: avoid using input BAM files larger than 8GB for best performance
features <- file.path(output_folder, "NA_peaks.narrowPeak")
sc_atac_feature_counting (fragment_file = file.path(output_folder, "fragments.bed"),
feature_input = features,
bam_tags = list(bc="CB", mb="OX"),
feature_type = "peak",
organism = "hg38",
yieldsize = 1000000,
exclude_regions = TRUE,
output_folder = output_folder,
fix_chr = "none"
)If genome_bin approach is selected, the following can be used:
reference <- file.path(data.folder, "small_chr21.fa")
sc_atac_feature_counting (fragment_file = file.path(output_folder, "fragments.bed"),
feature_input = reference,
bam_tags = list(bc="CB", mb="OX"),
feature_type = "genome_bin",
organism = "hg38",
cell_calling = FALSE,
yieldsize = 1000000,
exclude_regions = TRUE,
output_folder = output_folder,
fix_chr = "none"
)Call calling is performed on the data prior to finding overlaps. The
cell calling methods available are the emptyDrops method
from DropletUtils,
filter which filters the cells based on various QC
cut-offs, and cellranger which models the signal and noise
as a mixture of two negative binomial distributions, though the
filter method is recommended as it is most commonly used
and the best-suited for scATAC-seq data.
For the filter cell-calling method, there are a variety
of QC metrics that can be used, including:
min_uniq_fragsmax_uniq_fragsmin_frac_peakmin_frac_tssmin_frac_enhancermin_frac_promotermax_frac_mito
sc_atac_feature_counting (fragment_file = file.path(output_folder, "fragments.bed"),
feature_input = features,
bam_tags = list(bc="CB", mb="OX"),
feature_type = "peak",
organism = "hg38",
cell_calling = "filter",
min_uniq_frags = 0,
min_frac_peak = 0,
min_frac_promoter = 0,
yieldsize = 1000000,
exclude_regions = TRUE,
output_folder = output_folder,
fix_chr = "none"
)## `peak`, `tss` or `gene` feature_type is selected for feature input## Creating Galignment object for the feature input ...## hg38is a recognized organism. Using annotation files in repository.## Generating TSS plot data## Generating feature-barcode fragment count matrix## Average no. of fragments per barcode: 1.34885245901639## Raw feature matrix generated: /tmp/RtmpPhWjmQ/scPipeATACVignette/unfiltered_feature_matrix.rds## Calculating TSS enrichment scores## Generating QC metrics for cells## `filter` function is used for cell calling ...## cell called and stored in /tmp/RtmpPhWjmQ/scPipeATACVignette## Cells with less than 200 counts are filtered out.## Number of cells to remove:32## No cells were filtered out since otherwise there would be too few left.## features with less than 10 counts are filtered out.## Number of features to remove:8## No features were filtered out since otherwise there would be too few left.## making sparse## Sparse matrix generated## Sparse count matrix is saved in
## /tmp/RtmpPhWjmQ/scPipeATACVignette/sparse_matrix.mtx## Binary matrix is saved in:
## /tmp/RtmpPhWjmQ/scPipeATACVignette/binary_matrix.rds## Computing feature QC metrics## writing to csv## SCE object is saved in:
## /tmp/RtmpPhWjmQ/scPipeATACVignette/scPipe_atac_SCEobject.rds## sc_atac_feature_counting completed in 6.82266473770142 secondsMapping quality (MAPQ) value can be set to filter data further in
step (default: 30). Additionally, regions that are considered anomalous,
unstructured, or high signal in next-generation sequencing experiments
are excluded using an inbuilt excluded_regions file
(available are for organism hg19,
hg38, and mm10) or a user provided
excluded_regions_filename file in BED format.
The discrepancy of chr between the alignment file and
the feature file/excluded regions file can also be fixed (using
fix_char parameter) if needed by adding the string
chr to the beginning of either the features, and/or
excluded regions. So the options for fix_char are “none”,
“excluded_regions”, “feature”, “both”.
NOTE genome_bin approach may be more
reliable in detecting sensitive features than using a pseudo-bulk peak
calling approach, hence it would make the function slower as well.
The sc_atac_feature_counting function generates a matrix
format of the feature by cell matrix (features as rows, cell barcodes as
columns) in several formats (raw, sparse, binary, jaccard) that can be
used downstream to generate a singleCellExperiment, SCE
object.
feature_matrix <- readRDS(file.path(output_folder, "unfiltered_feature_matrix.rds"))
dplyr::glimpse(feature_matrix)## Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
## ..@ i : int [1:1559] 3 5 3 3 3 3 3 3 3 1 ...
## ..@ p : int [1:1526] 0 1 2 3 4 5 6 7 8 9 ...
## ..@ Dim : int [1:2] 11 1525
## ..@ Dimnames:List of 2
## .. ..$ feature: chr [1:11] "chr21:284477-284941" "chr21:287515-288874" "chr21:325151-325387" "chr21:414105-416460" ...
## .. ..$ barcode: chr [1:1525] "AAACGAACAGTAAGCG" "AAACGAACATGGATGG" "AAACGAATCGGACGAA" "AAACGAATCTGTGTGA" ...
## ..@ x : num [1:1559] 1 1 1 1 1 1 2 1 1 1 ...
## ..@ factors : list()## Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
## ..@ i : int [1:50] 3 3 3 3 3 1 3 3 3 3 ...
## ..@ p : int [1:33] 0 1 2 3 4 5 7 8 9 11 ...
## ..@ Dim : int [1:2] 11 32
## ..@ Dimnames:List of 2
## .. ..$ feature: chr [1:11] "chr21:284477-284941" "chr21:287515-288874" "chr21:325151-325387" "chr21:414105-416460" ...
## .. ..$ barcode: chr [1:32] "AAATGAGCAATCAGGG" "AACTTGGCATGGCCGT" "AATACGCAGTTGGAAT" "ACAGCGCAGATGCGCA" ...
## ..@ x : num [1:50] 7 12 16 3 6 1 21 6 6 9 ...
## ..@ factors : list()Generating the Single-cell Experiment (SCE) object
Easiest way to generate a SCE object is to run
sc_atac_create_sce with the input_folder
parameter. However, if the default dir name was used for previous steps
simply running sc_atac_create_sce() would produce a
sce object in the default location
(i.e. scPipe-atac-output)
sce <- sc_atac_create_sce(input_folder = output_folder,
organism = "hg38",
feature_type = "peak",
pheno_data = NULL,
report = FALSE)## SCE object is saved in:
## /tmp/RtmpPhWjmQ/scPipeATACVignette/scPipe_atac_SCEobject.rdsWe have now completed the preprocessing steps. The feature count
matrix is available as a .rds file in
scPipe-atac-output/feature_matrix.rds and quality control
statistics are saved in the
scPipe-atac-output/scPipe_atac_stats folder. This data is
useful for later quality control as well (QC).
If the report parameter is set to TRUE then an HTML
report with various statistics and plots is generated. A partial
screenshot of the report is shown below.
A convenience function for running the whole pipeline
A function sc_atac_pipeline has been included to make it
easier to run the whole pipeline.
data.folder <- system.file("extdata", package = "scPipe", mustWork = TRUE)
output_folder <- tempdir()
sce <- sc_atac_pipeline(r1 = file.path(data.folder, "small_chr21_R1.fastq.gz"),
r2 = file.path(data.folder, "small_chr21_R3.fastq.gz"),
barcode_fastq = file.path(data.folder, "small_chr21_R2.fastq.gz"),
organism = "hg38",
reference = file.path(data.folder, "small_chr21.fa"),
feature_type = "peak",
remove_duplicates = FALSE,
min_uniq_frags = 0, # set to 0 as the sample dataset only has a small number of reads
min_frac_peak = 0,
min_frac_promoter = 0,
output_folder = output_folder)Downstream analysis
Since the scater and scran packages both use the SingleCellExperiment class, it will be easy to further process this data using these packages for normalization and visualization. Other packages such as SC3 may be useful for clustering and MAST and edgeR for differential expression analysis.
Session Information
## R version 4.4.2 (2024-10-31)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=en_US.UTF-8
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [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
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] scPipe_2.7.1 SingleCellExperiment_1.29.1
## [3] SummarizedExperiment_1.37.0 Biobase_2.67.0
## [5] GenomicRanges_1.59.1 GenomeInfoDb_1.43.2
## [7] IRanges_2.41.2 S4Vectors_0.45.2
## [9] BiocGenerics_0.53.3 generics_0.1.3
## [11] MatrixGenerics_1.19.0 matrixStats_1.4.1
## [13] BiocStyle_2.35.0
##
## loaded via a namespace (and not attached):
## [1] DBI_1.2.3 bitops_1.0-9 httr2_1.0.7
## [4] biomaRt_2.63.0 rlang_1.1.4 magrittr_2.0.3
## [7] compiler_4.4.2 RSQLite_2.3.9 dir.expiry_1.15.0
## [10] png_0.1-8 vctrs_0.6.5 stringr_1.5.1
## [13] pkgconfig_2.0.3 crayon_1.5.3 Rsubread_2.21.1
## [16] fastmap_1.2.0 dbplyr_2.5.0 XVector_0.47.0
## [19] Rsamtools_2.23.1 rmarkdown_2.29 UCSC.utils_1.3.0
## [22] purrr_1.0.2 bit_4.5.0.1 xfun_0.49
## [25] zlibbioc_1.52.0 cachem_1.1.0 jsonlite_1.8.9
## [28] progress_1.2.3 blob_1.2.4 DelayedArray_0.33.3
## [31] reshape_0.8.9 BiocParallel_1.41.0 parallel_4.4.2
## [34] prettyunits_1.2.0 R6_2.5.1 bslib_0.8.0
## [37] stringi_1.8.4 reticulate_1.40.0 rtracklayer_1.67.0
## [40] jquerylib_0.1.4 Rcpp_1.0.13-1 knitr_1.49
## [43] org.Mm.eg.db_3.20.0 R.utils_2.12.3 Matrix_1.7-1
## [46] tidyselect_1.2.1 abind_1.4-8 yaml_2.3.10
## [49] codetools_0.2-20 curl_6.0.1 lattice_0.22-6
## [52] tibble_3.2.1 plyr_1.8.9 withr_3.0.2
## [55] basilisk.utils_1.19.0 KEGGREST_1.47.0 evaluate_1.0.1
## [58] BiocFileCache_2.15.0 xml2_1.3.6 mclust_6.1.1
## [61] Biostrings_2.75.3 pillar_1.10.0 BiocManager_1.30.25
## [64] filelock_1.0.3 RCurl_1.98-1.16 hms_1.1.3
## [67] glue_1.8.0 maketools_1.3.1 tools_4.4.2
## [70] BiocIO_1.17.1 robustbase_0.99-4-1 sys_3.4.3
## [73] data.table_1.16.4 GenomicAlignments_1.43.0 buildtools_1.0.0
## [76] XML_3.99-0.17 grid_4.4.2 tidyr_1.3.1
## [79] AnnotationDbi_1.69.0 GenomeInfoDbData_1.2.13 basilisk_1.19.0
## [82] Rhtslib_3.3.1 restfulr_0.0.15 cli_3.6.3
## [85] rappdirs_0.3.3 S4Arrays_1.7.1 dplyr_1.1.4
## [88] DEoptimR_1.1-3-1 R.methodsS3_1.8.2 sass_0.4.9
## [91] digest_0.6.37 SparseArray_1.7.2 org.Hs.eg.db_3.20.0
## [94] rjson_0.2.23 memoise_2.0.1 htmltools_0.5.8.1
## [97] R.oo_1.27.0 lifecycle_1.0.4 httr_1.4.7
## [100] bit64_4.5.2- Introduction
- Workflow
- Getting started
- Fastq reformatting
- Aligning reads to a reference genome
- Demultiplexing the BAM file
- Remove duplicates
- Gemerating a fragment file
- Peak calling
- Assigning reads to features and feature counting
- Generating the Single-cell Experiment (SCE) object
- A convenience function for running the whole pipeline
- Downstream analysis
- Session Information