Case Study: Preprocessing CEL-seq2 data

CEL-seq2 is a modified version of the CEL-seq method with higher sensitivity and lower cost. It uses a polyT primer with index sequences to capture a cell’s mRNA. The index sequence consists of the a cell specific index, which is designed, and a molecular index, also called unique molecular identifier (UMI), which is a random sequence that should be different for each primer (the number of distinct UMIs will depend on its length). The ouput fastq files from a CEL-seq2 experiment is paired-ended where the first read contains the transcript and second read contains the cellindex+UMI with some polyA sequence.

Getting started

We begin by creating a folder to store the results.

library(scPipe)
library(SingleCellExperiment)
data_dir = tempdir()

Organising the files required

To process the data, we need the genome fasta file, gff3 exon annotation and a cell barcode annotation. The barcode annotation should be a .csv file with at least two columns, where the first column has the cell id and the second column contains the barcode sequence. We use ERCC spike-in genes for this demo. All files can be found in the extdata folder of the scPipe package:

# file path:
ERCCfa_fn = system.file("extdata", "ERCC92.fa", package = "scPipe")
ERCCanno_fn = system.file("extdata", "ERCC92_anno.gff3", package = "scPipe")
barcode_annotation_fn = system.file("extdata", "barcode_anno.csv", package = "scPipe")

The read structure for this example is paired-ended, with the first longer read mapping to transcripts and second shorter read consisting of 6bp UMIs followed by 8bp cell barcodes. NOTE: by convention, scPipe always assumes read1 refers to the read with the transcript sequence, which is usually the longer read. These data are also available in the in extdata folder:

fq_R1 = system.file("extdata", "simu_R1.fastq.gz", package = "scPipe")
fq_R2 = system.file("extdata", "simu_R2.fastq.gz", package = "scPipe")

Data Preprocessing

Fastq reformatting

The pipeline starts with fastq file reformatting. We move the barcode and UMI sequences to the read name and leave the transcript sequence as is. This outputs a read name that looks like @[barcode_sequence]*[UMI_sequence]#[readname] … The read structure of read 2 in our example dataset has the 8 bp long cell barcode starting at position 6 and the 6 bp long UMI sequence starting at the first position. So the read structure will be : list(bs1=-1, bl1=0, bs2=6, bl2=8, us=0, ul=6). bs1=-1, bl1=0 means we don’t have an index in read 1 so we set a negative value as its start position and give it zero length. bs2=6, bl2=8 means we have an index in read two which starts at position 6 in the read and is 8 bases in length. us=0, ul=6 means we have a UMI at the start of read 2 which is 6 bases long. NOTE: we use a zero based index system, so the indexing of the sequence starts at zero.

sc_trim_barcode(file.path(data_dir, "combined.fastq.gz"),
                fq_R1,
                fq_R2,
                read_structure = list(bs1=-1, bl1=0, bs2=6, bl2=8, us=0, ul=6))

## trimming fastq file...
## pass QC: 46009
## removed_have_N: 0
## removed_low_qual: 0
## time elapsed: 117 milliseconds

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. The Rsubread package is not available on Windows, so the following alignment step will only run on Mac OSX or Linux.

Rsubread::buildindex(basename=file.path(data_dir, "ERCC_index"), reference=ERCCfa_fn)

## 
##         ==========     _____ _    _ ____  _____  ______          _____  
##         =====         / ____| |  | |  _ \|  __ \|  ____|   /\   |  __ \ 
##           =====      | (___ | |  | | |_) | |__) | |__     /  \  | |  | |
##             ====      \___ \| |  | |  _ <|  _  /|  __|   / /\ \ | |  | |
##               ====    ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
##         ==========   |_____/ \____/|____/|_|  \_\______/_/    \_\_____/
##        Rsubread 2.16.0
## 
## //================================= setting ==================================\\
## ||                                                                            ||
## ||                Index name : ERCC_index                                     ||
## ||               Index space : base space                                     ||
## ||               Index split : no-split                                       ||
## ||          Repeat threshold : 100 repeats                                    ||
## ||              Gapped index : no                                             ||
## ||                                                                            ||
## ||       Free / total memory : 78.3GB / 125.4GB                               ||
## ||                                                                            ||
## ||               Input files : 1 file in total                                ||
## ||                             o ERCC92.fa                                    ||
## ||                                                                            ||
## \\============================================================================//
## 
## //================================= Running ==================================\\
## ||                                                                            ||
## || Check the integrity of provided reference sequences ...                    ||
## || There were 1 notes for reference sequences.                                ||
## || The notes can be found in the log file, '/tmp/Rtmpwxfbe7/ERCC_index.log'.  ||
## || Scan uninformative subreads in reference sequences ...                     ||
## || 1 uninformative subreads were found.                                       ||
## || These subreads were excluded from index building.                          ||
## || Estimate the index size...                                                 ||
## ||    8%,   0 mins elapsed, rate=126.4k bps/s                                 ||
## ||   49%,   0 mins elapsed, rate=731.1k bps/s                                 ||
## ||   66%,   0 mins elapsed, rate=960.4k bps/s                                 ||
## || 3.0 GB of memory is needed for index building.                             ||
## || Build the index...                                                         ||
## ||    8%,   0 mins elapsed, rate=12.6k bps/s                                  ||
## ||   49%,   0 mins elapsed, rate=74.7k bps/s                                  ||
## ||   66%,   0 mins elapsed, rate=99.3k 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.2 minutes.                      ||
## ||          Index /tmp/Rtmpwxfbe7/ERCC_index was successfully built.          ||
## ||                                                                            ||
## \\============================================================================//

Rsubread::align(index=file.path(data_dir, "ERCC_index"),
  readfile1=file.path(data_dir, "combined.fastq.gz"),
  output_file=file.path(data_dir, "out.aln.bam"), phredOffset=64)

## 
##         ==========     _____ _    _ ____  _____  ______          _____  
##         =====         / ____| |  | |  _ \|  __ \|  ____|   /\   |  __ \ 
##           =====      | (___ | |  | | |_) | |__) | |__     /  \  | |  | |
##             ====      \___ \| |  | |  _ <|  _  /|  __|   / /\ \ | |  | |
##               ====    ____) | |__| | |_) | | \ \| |____ / ____ \| |__| |
##         ==========   |_____/ \____/|____/|_|  \_\______/_/    \_\_____/
##        Rsubread 2.16.0
## 
## //================================= setting ==================================\\
## ||                                                                            ||
## || Function      : Read alignment (RNA-Seq)                                   ||
## || Input file    : combined.fastq.gz                                          ||
## || Output file   : out.aln.bam (BAM)                                          ||
## || Index name    : ERCC_index                                                 ||
## ||                                                                            ||
## ||                    ------------------------------------                    ||
## ||                                                                            ||
## ||                               Threads : 1                                  ||
## ||                          Phred offset : 64                                 ||
## ||                             Min votes : 3 / 10                             ||
## ||                        Max mismatches : 3                                  ||
## ||                      Max indel length : 5                                  ||
## ||            Report multi-mapping reads : yes                                ||
## || Max alignments per multi-mapping read : 1                                  ||
## ||                                                                            ||
## \\============================================================================//
## 
## //=============== Running (24-Oct-2023 19:01:13, pid=3734380) ================\\
## ||                                                                            ||
## || 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 [1,1]                    ||
## || 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.                                               ||
## ||   23% completed, 0.5 mins elapsed, rate=119.8k reads per second            ||
## ||   43% completed, 0.5 mins elapsed, rate=144.8k reads per second            ||
## ||   62% completed, 0.5 mins elapsed, rate=158.1k reads per second            ||
## ||   82% completed, 0.5 mins elapsed, rate=162.9k reads per second            ||
## ||   72% completed, 0.5 mins elapsed, rate=0.7k reads per second              ||
## ||   78% completed, 0.5 mins elapsed, rate=0.8k reads per second              ||
## ||   83% completed, 0.5 mins elapsed, rate=0.9k reads per second              ||
## ||   88% completed, 0.5 mins elapsed, rate=0.9k reads per second              ||
## ||   92% completed, 0.5 mins elapsed, rate=1.0k reads per second              ||
## ||   97% completed, 0.5 mins elapsed, rate=1.0k reads per second              ||
## ||  102% completed, 0.5 mins elapsed, rate=1.1k reads per second              ||
## ||  107% completed, 0.5 mins elapsed, rate=1.1k reads per second              ||
## ||  113% completed, 0.5 mins elapsed, rate=1.2k reads per second              ||
## ||                                                                            ||
## ||                           Completed successfully.                          ||
## ||                                                                            ||
## \\====================================    ====================================//
## 
## //================================   Summary =================================\\
## ||                                                                            ||
## ||                 Total reads : 46,009                                       ||
## ||                      Mapped : 46,009 (100.0%)                              ||
## ||             Uniquely mapped : 46,009                                       ||
## ||               Multi-mapping : 0                                            ||
## ||                                                                            ||
## ||                    Unmapped : 0                                            ||
## ||                                                                            ||
## ||                      Indels : 0                                            ||
## ||                                                                            ||
## ||                Running time : 0.5 minutes                                  ||
## ||                                                                            ||
## \\============================================================================//

##                       out.aln.bam
## Total_reads                 46009
## Mapped_reads                46009
## Uniquely_mapped_reads       46009
## Multi_mapping_reads             0
## Unmapped_reads                  0
## Indels                          0

Assigning reads to annotated exons

After the read alignment, we assign reads to exons based on the annotation provided using the sc_exon_mapping function. Currently scPipe only supports annotation in gff3 or bed format.

sc_exon_mapping(file.path(data_dir, "out.aln.bam"),
            file.path(data_dir, "out.map.bam"),
            ERCCanno_fn)

## adding annotation files...
## time elapsed: 0 milliseconds
## 
## annotating exon features...
## updating progress every 3 minutes...
## 46009 reads processed, 46k reads/sec
## number of read processed: 46009
## unique map to exon: 46009 (100.00%)
## ambiguous map to multiple exon: 0 (0.00%)
## map to intron: 0 (0.00%)
## not mapped: 0 (0.00%)
## unaligned: 0 (0.00%)

De-multiplexing data and counting genes

Next we use the sc_demultiplex function to split the reads into separate .csv files for each cell in the /count subfolder. Each file contains three columns and each row corresponds to a distinct read. The first column contains the gene ID that the read maps to, the second column gives the UMI sequence for that read and the third column is the distance to the transcript end position (TES) in bp. This file will be used for UMI deduplication and to generate a gene count matrix by calling the sc_gene_counting function.

sc_demultiplex(file.path(data_dir, "out.map.bam"), data_dir, barcode_annotation_fn, has_UMI=FALSE)

## demultiplexing reads by barcode...
## Warning: mitochondrial chromosome not found using chromosome name `MT`.
## time elapsed: 72 milliseconds

sc_gene_counting(data_dir, barcode_annotation_fn)

## summarising gene counts...
## time elapsed: 62 milliseconds

We have now completed the preprocessing steps. The gene count matrix is available as a .csv file in data_dir/gene_count.csv and quality control statistics are saved in the data_dir/stat folder. These data are useful for later quality control (QC).

Preprocessing data generated by other protocols

The greatest difference between scRNA-seq protocols is in their read structure, which is specified by the read_structure argument of the sc_trim_barcode function in scPipe. The major difference between CEL-seq and Drop-seq/Chromium 10x is that the cell barcode is unknown in advance for the latter protocols, so an extra step sc_detect_bc is required before running sc_demultiplex, to identify and generate the cell barcode annotation. scPipe is not optimized for the full-lengh protocols like SMART-seq, but it can process the data generated by such protocols by setting the has_UMI to FALSE in the sc_demultiplex function.

Quality Control

The easiest way to create a SingleCellExperiment object from the output of scPipe preprocessing is using create_sce_by_dir function, that will read in the gene count matrix together with the QC information available in the stat folder.

sce = create_sce_by_dir(data_dir)
dim(sce)

## [1] 92 10

The dataset we analysed above used ERCC simulated reads from 10 cells of perfect quality. In order to demonstrate QC on a more reaslitic example, we have included an experimental dataset in the scPipe software generated by Dr Christine Biben. This experiment profiles gene expression in 383 blood cells for a subset of 1000 genes (to save space). We first create a SingleCellExperiment object directly without using the create_sce_by_dir function:

data("sc_sample_data")
data("sc_sample_qc")
sce = SingleCellExperiment(assays = list(counts = as.matrix(sc_sample_data))) # generate new sce with gene count matrix
QC_metrics(sce) = sc_sample_qc
demultiplex_info(sce) = cell_barcode_matching

## organism/gene_id_type not provided. Make a guess:mmusculus_gene_ensembl/ensembl_gene_id

UMI_dup_info(sce) = UMI_duplication

There are several plots we can create to assess the overall quality of the experiment. The first is to look at the cell barcode demultiplexing statistics. To do this we generate a bar plot that shows the percentage of reads that uniquely match to the cell barcodes, as well as the unmatched proportion of reads and their alignment rates to introns and exons. If we observe a large proportion of unmatched reads that map to exons, this indicates a failure of the cell barcode demultiplexing.

plot_demultiplex(sce)

A second plot shows the duplication rate which can be used to evaluate read depth. UMIs are routinely used to mark individual molecules and after PCR amplification, the same molecule with have multiple copys, which can be identified and removed if they are observed to have the same UMI sequence. Therefore, the copy number of UMIs is an indication of the PCR amplification rate.

plot_UMI_dup(sce)

## `geom_smooth()` using formula = 'y ~ x'

Next we calculate QC metrics and use the detect_outlier function to identify poor quality cells. The detect_outlier function has argument comp to define the maximum component of the gaussian mixture model. Using the default value of 1 is generally sufficient, but in cases where the data are heterogeneous in terms of quality control metrics, setting this value to 2 or 3 can give better results. This function will remove low quality cells if type="low", large cells if type="high" or both when type="both". The conf argument specifies the lower and upper confidence intervals for outliers and detect_outlier is insensistive to the interval values.

sce = calculate_QC_metrics(sce)
sce = detect_outlier(sce)

We can plot the alignment statistics per sample as a barplot using

plot_mapping(sce, percentage = TRUE, dataname = "sc_sample_data")

and generate pairwise plots of the QC metrics as follows:

plot_QC_pairs(sce)

The final step will be to remove the low quality cells by the remove_outliers function.

sce = remove_outliers(sce)
dim(sce)

## [1] 1000  360

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.

sessionInfo()

## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.18-bioc/R/lib/libRblas.so 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=en_US.UTF-8    
##  [3] LC_TIME=en_GB              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: America/New_York
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] Matrix_1.6-1.1              XVector_0.42.0             
##  [3] scPipe_2.2.0                SingleCellExperiment_1.24.0
##  [5] SummarizedExperiment_1.32.0 Biobase_2.62.0             
##  [7] GenomicRanges_1.54.0        GenomeInfoDb_1.38.0        
##  [9] IRanges_2.36.0              S4Vectors_0.40.0           
## [11] BiocGenerics_0.48.0         MatrixGenerics_1.14.0      
## [13] matrixStats_1.0.0           BiocStyle_2.30.0           
## 
## loaded via a namespace (and not attached):
##   [1] DBI_1.1.3                bitops_1.0-7             biomaRt_2.58.0          
##   [4] rlang_1.1.1              magrittr_2.0.3           compiler_4.3.1          
##   [7] RSQLite_2.3.1            mgcv_1.9-0               dir.expiry_1.10.0       
##  [10] png_0.1-8                vctrs_0.6.4              stringr_1.5.0           
##  [13] pkgconfig_2.0.3          crayon_1.5.2             Rsubread_2.16.0         
##  [16] fastmap_1.1.1            dbplyr_2.3.4             labeling_0.4.3          
##  [19] utf8_1.2.4               Rsamtools_2.18.0         rmarkdown_2.25          
##  [22] purrr_1.0.2              bit_4.0.5                xfun_0.40               
##  [25] zlibbioc_1.48.0          cachem_1.0.8             jsonlite_1.8.7          
##  [28] progress_1.2.2           blob_1.2.4               DelayedArray_0.28.0     
##  [31] reshape_0.8.9            BiocParallel_1.36.0      parallel_4.3.1          
##  [34] prettyunits_1.2.0        R6_2.5.1                 RColorBrewer_1.1-3      
##  [37] bslib_0.5.1              stringi_1.7.12           GGally_2.1.2            
##  [40] reticulate_1.34.0        rtracklayer_1.62.0       jquerylib_0.1.4         
##  [43] Rcpp_1.0.11              bookdown_0.36            knitr_1.44              
##  [46] org.Mm.eg.db_3.18.0      R.utils_2.12.2           splines_4.3.1           
##  [49] tidyselect_1.2.0         abind_1.4-5              yaml_2.3.7              
##  [52] codetools_0.2-19         curl_5.1.0               lattice_0.22-5          
##  [55] tibble_3.2.1             plyr_1.8.9               withr_2.5.1             
##  [58] basilisk.utils_1.14.0    KEGGREST_1.42.0          evaluate_0.22           
##  [61] BiocFileCache_2.10.0     xml2_1.3.5               mclust_6.0.0            
##  [64] Biostrings_2.70.0        pillar_1.9.0             BiocManager_1.30.22     
##  [67] filelock_1.0.2           generics_0.1.3           RCurl_1.98-1.12         
##  [70] ggplot2_3.4.4            hms_1.1.3                munsell_0.5.0           
##  [73] scales_1.2.1             glue_1.6.2               tools_4.3.1             
##  [76] BiocIO_1.12.0            robustbase_0.99-0        data.table_1.14.8       
##  [79] GenomicAlignments_1.38.0 XML_3.99-0.14            grid_4.3.1              
##  [82] tidyr_1.3.0              colorspace_2.1-0         AnnotationDbi_1.64.0    
##  [85] nlme_3.1-163             GenomeInfoDbData_1.2.11  basilisk_1.14.0         
##  [88] Rhtslib_2.4.0            restfulr_0.0.15          cli_3.6.1               
##  [91] rappdirs_0.3.3           fansi_1.0.5              S4Arrays_1.2.0          
##  [94] dplyr_1.1.3              gtable_0.3.4             DEoptimR_1.1-3          
##  [97] hash_2.2.6.3             R.methodsS3_1.8.2        sass_0.4.7              
## [100] digest_0.6.33            SparseArray_1.2.0        farver_2.1.1            
## [103] org.Hs.eg.db_3.18.0      rjson_0.2.21             memoise_2.0.1           
## [106] htmltools_0.5.6.1        R.oo_1.25.0              lifecycle_1.0.3         
## [109] httr_1.4.7               MASS_7.3-60              bit64_4.0.5

scPipe: a flexible data preprocessing pipeline for 3’ end scRNA-seq data

Luyi Tian

2023-10-24

Introduction