--- title: "scPipe: a flexible data preprocessing pipeline for scATAC-seq data" author: "Shani Amarasinghe, Phil Yang" date: "`r Sys.Date()`" output: BiocStyle::html_document: toc_float: true vignette: > %\VignetteIndexEntry{scPipe: a flexible data preprocessing pipeline for scATAC-seq data} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- # 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. ```{r, echo=FALSE, warning=FALSE, message=FALSE, out.width="100%"} knitr::include_graphics(system.file("scATAC-seq_workflow.png", package = "scPipe")) ``` # 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. ```{r, message=FALSE} library(scPipe) data.folder <- system.file("extdata", package = "scPipe", mustWork = TRUE) ``` Specify the output folder. ```{r} output_folder <- paste0(tempdir(), "/scPipeATACVignette") ``` # 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: ```{r} # 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. ```{r} sc_atac_trim_barcode (r1 = r1, r2 = r2, bc_file = barcode_fastq, rmN = FALSE, rmlow = FALSE, output_folder = output_folder) ``` This 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. ```{r,eval=FALSE} 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. ```{r} 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) ``` # 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. ```{r} sorted_tagged_bam <- sc_atac_bam_tagging (inbam = aligned_bam, output_folder = output_folder, bam_tags = list(bc="CB", mb="OX"), nthreads = 12) ``` # Remove 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. ```{r, eval=FALSE} removed <- sc_atac_remove_duplicates(sorted_tagged_bam, output_folder = output_folder) if (!isFALSE(removed)) sorted_tagged_bam <- removed ``` # 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. ```{r} sc_atac_create_fragments(inbam = sorted_tagged_bam, output_folder = output_folder) ``` # Peak 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. ```{r, eval=FALSE} # CHECK IF THE .narrowPeak file is small enough to include in the package, # otherwise, we need to make this basilisk call work?? sc_atac_peak_calling(inbam = sorted_tagged_bam, ref = reference, genome_size = NULL, output_folder = output_folder) ``` # 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 ```{r, eval=FALSE} 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: ```{r,eval=FALSE} 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](https://bioconductor.org/packages/release/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_frags` * `max_uniq_frags` * `min_frac_peak` * `min_frac_tss` * `min_frac_enhancer` * `min_frac_promoter` * `max_frac_mito` ```{r,eval=FALSE} features <- file.path(output_folder, "NA_peaks.narrowPeak") ``` ```{r, eval=TRUE, echo=FALSE} features <- file.path(data.folder, "NA_peaks.narrowPeak") ``` ```{r} 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" ) ``` Mapping 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. ```{r} feature_matrix <- readRDS(file.path(output_folder, "unfiltered_feature_matrix.rds")) dplyr::glimpse(feature_matrix) sparseM <- readRDS(file.path(output_folder, "sparse_matrix.rds")) dplyr::glimpse(sparseM) ``` # 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`) ```{r} sce <- sc_atac_create_sce(input_folder = output_folder, organism = "hg38", feature_type = "peak", pheno_data = NULL, report = FALSE) ``` We 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. ```{r, echo=FALSE, warning=FALSE, message=FALSE, out.width="100%"} knitr::include_graphics(system.file("report_example.png", package = "scPipe")) ``` # 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. ```{r eval = FALSE} 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](https://bioconductor.org/packages/release/scater) and [scran](https://bioconductor.org/packages/release/scran) packages both use the [SingleCellExperiment](https://bioconductor.org/packages/release/SingleCellExperiment) class, it will be easy to further process this data using these packages for normalization and visualization. Other packages such as [SC3](https://bioconductor.org/packages/release/SC3) may be useful for clustering and [MAST](https://bioconductor.org/packages/release/MAST) and [edgeR](https://bioconductor.org/packages/release/edgeR) for differential expression analysis. ```{r, echo=FALSE, results='hide'} # cleanup tempfiles unlink(output_folder, recursive=TRUE) ``` # Session Information {-} ```{r} sessionInfo() ```