Generation of the Count Matrix

This lesson explains how raw single-cell RNA-seq sequencing data is processed into a cell-by-gene count matrix. This includes descriptions of how UMIs are used to identify PCR duplicates and the overall structure of reads to quantify gene expression per cell. The lesson then outlines the standard preprocessing workflow: formatting reads and filtering noisy barcodes, demultiplexing samples, mapping or pseudo-mapping to the transcriptome, collapsing UMIs and quantifying reads to produce the count matrix.

Keywords

Single-cell RNA-seq, Count matrix, UMIs, Cell barcodes, Demultiplexing, Mapping, Cell Ranger

Approximate time: 90 minutes

Learning objectives

In this lesson, we will:

• Describe how raw scRNA-seq reads are processed into a cell-by-gene count matrix.
  
• Explain the roles of cell barcodes, sample indices, and UMIs in assigning and de-duplicating reads.
  
• Outline the main preprocessing steps: formatting and filtering barcodes, alignment, and collapsing UMIs for quantification.

Library preparation

Depending on the library preparation method used, the RNA sequences (also referred to as reads or tags) will be derived either from the 3’ ends (or 5’ ends) of the transcripts (10X Genomics, CEL-seq2, Drop-seq, inDrops) or from full-length transcripts (Smart-seq).

Figure 1: Schematic of different library preparation methods for single-cell RNA-sequencing.
Image credit: Papalexi et al. (2018)

The choice of method depends on the biological question of interest. The following advantages are listed for each approach:

• 3’ (or 5’)-end sequencing:
  • More accurate quantification through the use of unique molecular identifiers distinguishing biological duplicates from amplification (PCR) duplicates
  • Larger number of cells sequenced allows better identity of cell type populations
  • Cheaper per cell cost
  • Best results for experiments with > 10,000 cells
• Full length sequencing:
  • Detection of isoform-level differences in expression
  • Identification of allele-specific differences in expression
  • Deeper sequencing of a smaller number of cells
  • Best for samples with few cells

Many of the same analysis steps need to occur for 3’-end sequencing as for full-length, but 3’ protocols have been increasing in popularity and consist of a few more steps in the analysis. Therefore, our materials are going to detail the analysis of data from these 3’ protocols with a focus on the droplet-based methods (inDrops, Drop-seq, 10X Genomics).

3’-end reads

For the analysis of scRNA-seq data, it is helpful to understand what information is present in each of the reads and how we use it moving forward through the analysis. The following steps pertain to all droplet-based methods.

For the 3’-end sequencing methods, reads originating from different molecules of the same transcript would have originated only from the 3’ end of the transcripts, so would have a high likelihood of having the same sequence. However, the PCR step during library preparation could also generate read duplicates. To determine whether a read is a biological or technical duplicate, these methods use unique molecular identifiers, or UMIs.

• Reads with different UMIs mapping to the same transcript were derived from different molecules and are biological duplicates - each read should be counted.
• Reads with the same UMI + cell barcode and map to the same gene originated from the same molecule. These are technical duplicates (PCR duplicates) and should be collapse and counted as a single read.
• In the image below, the reads for ACTB should be collapsed and counted as a single read, while the reads for ARL1 should each be counted.

Figure 2: Simplified representation of UMIs and cell barcodes are used to identify uniquely sequenced reads.
Image credit: Adapted from Macosko et al. (2015)

So we know that we need to keep track of the UMIs, but what other information do we need to properly quantify the expression in each gene in each of the cells in our samples? Regardless of droplet method, the following are required for proper quantification at the cellular level:

Figure 3: Architecture of the molecule sequenced, showcasing the primer, barcode, UMI, and read.
Image credit: Sarah Boswell

• Sample index: determines which sample the read originated from (red bottom arrow)
  • Added during library preparation - needs to be documented
• Cellular barcode: determines which cell the read originated from (purple top arrow)
  • Each library preparation method has a stock of cellular barcodes used during the library preparation
• Unique molecular identifier (UMI): determines which transcript molecule the read originated from
  • The UMI will be used to collapse PCR duplicates (purple bottom arrow)
• Sequencing read1: the Read1 sequence (red top arrow)
• Sequencing read2: the Read2 sequence (purple bottom arrow)

Single-cell RNA-seq workflow

The scRNA-seq method will determine how to parse the barcodes and UMIs from the sequencing reads. So, although a few of the specific steps will slightly differ, the overall workflow will generally follow the same steps regardless of method. The general workflow is shown below:

Figure 4: Standard computational workflow of a scRNA-seq dataset.
Image credit: Luecken et al. (2019)

The steps of the workflow are:

• Generation of the count matrix (method-specific steps): formatting reads, demultiplexing samples, mapping and quantification
• Quality control of the raw counts: filtering of poor quality cells
• Clustering of filtered counts: clustering cells based on similarities in transcriptional activity (cell types = different clusters)
• Marker identification and cluster annotation: identifying gene markers for each cluster and annotating known cell type clusters
• Optional downstream steps

Regardless of the analysis being done, conclusions about a population based on a single sample per condition are not trustworthy. BIOLOGICAL REPLICATES ARE STILL NEEDED! That is, if you want to make conclusions that correspond to the population and not just the single sample.

Generation of the count matrix

We are going to start by discussing the first part of this workflow, which is generating the count matrix from the raw sequencing data. We will focus on the 3’ end sequencing used by droplet-based methods, such as inDrops, 10X Genomics, and Drop-seq.

Figure 5: Initial steps of a scRNA-seq analysis.

After sequencing, the sequencing facility will either output the raw sequencing data as BCL or FASTQ format or will generate the count matrix. If the reads are in BCL format, then we will need to convert to FASTQ format. There is a useful command-line tool called bcl2fastq that can easily perform this conversion.

Demultiplexing

We do not demultiplex at this step in the workflow. You may have sequenced 6 samples, but the reads for all samples may be present all in the same BCL or FASTQ file.

The generation of the count matrix from the raw sequencing data will go through similar steps for many of the scRNA-seq methods.

Standard steps

Figure 6: Overview of count generation steps.

alevin is a command-line tool that estimates expression of scRNA-seq data for which the 3’ ends of transcripts were sequenced. umi-tools and zUMIs are additional tools that can perform these processes. These tools incorporate collapsing of UMIs to correct for amplification bias. The steps in this process include the following:

1. Formatting reads and filtering noisy cellular barcodes
2. Demultiplexing the samples
3. Mapping/pseudo-mapping to transcriptome
4. Collapsing UMIs and quantification of reads

Cell Ranger pipeline

If using 10X Genomics library preparation method, then the Cell Ranger pipeline would be used for all of the above steps.

1. Formatting reads and filtering noisy cellular barcodes

The FASTQ files can then be used to parse out the cell barcodes, UMIs, and sample barcodes. For droplet-based methods, many of the cellular barcodes will match a low number of reads (< 1000 reads) due to:

• Encapsulation of free floating RNA from dying cells
• Simple cells (RBCs, etc.) expressing few genes
• Cells that failed for some reason

These excess barcodes need to be filtered out of the sequence data prior to read alignment. To do this filtering, the ‘cellular barcode’ and the ‘molecular barcode’ are extracted and saved for each cell. For example, if using umi-tools, the information is added to the header line for each read, with the following format:

@HWI-ST808:130:H0B8YADXX:1:1101:2088:2222:CELL_GGTCCA:UMI_CCCT
AGGAAGATGGAGGAGAGAAGGCGGTGAAAGAGACCTGTAAAAAGCCACCGN
+
@@@DDBD>=AFCF+<CAFHDECII:DGGGHGIGGIIIEHGIIIGIIDHII#

The cellular barcodes are known values from the kit being used. Therefore, any read that does not match this whitelist after correction (allowing for acceptable mismatches due to sequencing mistakes) are removed.

2. Demultiplexing sample reads

The next step of the process is to demultiplex the samples, if sequencing more than a single sample. This is the one step of this process not handled by the umi-tools, but is accomplished by zUMIs. We would need to parse the reads to determine the sample barcode associated with each cell.

3. Mapping/pseudo-mapping to cDNAs

To determine which gene the read originated from, the reads are aligned using traditional (STAR) or light-weight methods (Kallisto/RapMap).

4. Collapsing UMIs and quantification of reads

The duplicate UMIs are collapsed, and only the unique UMIs are quantified using a tool like Kallisto or featureCounts. The resulting output is a cell by gene matrix of counts:

Figure 7: Resultant count matrix after pre-processing steps are completed.
Image credit: Lafzi et al. (2018)

Each value in the matrix represents the number of reads in a cell originating from the corresponding gene. Using the count matrix, we can explore and filter the data, keeping only the higher quality cells.

This is the starting point of most computational analyses!

---

Next Lesson >>

Back to Schedule