Single cell workflow

Identification of cell- and UMI barcodes from single-cell sequencing.

Introduction

This workflow extracts cell barcodes and UMIs from 10x-generated single cell libraries. It was initially created as a Nextflow port of Sockeye.

In brief, the workflow does the following:

• Adapter identification, fused read splitting and stranding.
• Mapping of reads to genomic reference.
• Gene and transcript read assignment.
• Cell barcode and UMI extraction and correction.
• Generation of gene and transcript count matrices for unique UMIs.
• Tagging BAM files with cell barcodes and UMIs.
• Calculation of library saturation.

This workflow supports the following 10x kits:

• 3’: v2/v3 and v4 (GEM-X)
• 5’: v1/v2
• multiome (gene expression only): v1
• visium spatial transcriptomics: v1

The BLAZE preprint provided useful benchmarking of the original sockeye implementation. This assisted in the selection of appropriate thresholds for cell cut-off and for defining the limits of the gene x cell matrix.

The isoform selection procedure used in this workflow was adapted from that found in the FLAMES package.

Compute requirements

Recommended requirements:

• CPUs = 64
• Memory = 256GB

Minimum requirements:

• CPUs = 8
• Memory = 32GB

Approximate run time: Approximately 8h for 120M reads with the recommended requirements.

ARM processor support: False

Install and run

These are instructions to install and run the workflow on command line. You can also access the workflow via the EPI2ME Desktop application.

The workflow uses Nextflow to manage compute and software resources, therefore Nextflow will need to be installed before attempting to run the workflow.

The workflow can currently be run using either Docker or Singularity to provide isolation of the required software. Both methods are automated out-of-the-box provided either Docker or Singularity is installed. This is controlled by the -profile parameter as exemplified below.

It is not required to clone or download the git repository in order to run the workflow. More information on running EPI2ME workflows can be found on our website.

The following command can be used to obtain the workflow. This will pull the repository in to the assets folder of Nextflow and provide a list of all parameters available for the workflow as well as an example command:

nextflow run epi2me-labs/wf-single-cell --help

To update a workflow to the latest version on the command line use the following command:

nextflow pull epi2me-labs/wf-single-cell

A demo dataset is provided for testing of the workflow. It can be downloaded and unpacked using the following commands:

wget https://ont-exd-int-s3-euwst1-epi2me-labs.s3.amazonaws.com/wf-single-cell/wf-single-cell-demo.tar.gz
tar -xzvf wf-single-cell-demo.tar.gz

The workflow can then be run with the downloaded demo data using:

nextflow run epi2me-labs/wf-single-cell \
    --expected_cells 100 \
    --fastq 'wf-single-cell-demo/chr17.fq.gz' \
    --kit '3prime:v3' \
    --ref_genome_dir 'wf-single-cell-demo' \
    --genes_of_interest 'wf-single-cell-demo/umap_plot_genes.csv' \
    -profile standard

For further information about running a workflow on the command line see https://labs.epi2me.io/wfquickstart/

Related protocols

This workflow is designed to take input sequences that have been produced from Oxford Nanopore Technologies devices.

Find related protocols in the Nanopore community.

• Library prep and sequencing protocol for the 10x 5’ kit
• Library prep and sequencing protocol for the 10x 3’ kit

Input example

This workflow accepts either FASTQ or BAM files as input.

The FASTQ or BAM input parameters for this workflow accept one of three cases: (i) the path to a single FASTQ or BAM file; (ii) the path to a top-level directory containing FASTQ or BAM files; (iii) the path to a directory containing one level of sub-directories which in turn contain FASTQ or BAM files. In the first and second cases (i and ii), a sample name can be supplied with --sample. In the last case (iii), the data is assumed to be multiplexed with the names of the sub-directories as barcodes. In this case, a sample sheet can be provided with --sample_sheet.

(i)                     (ii)                 (iii)    
input_reads.fastq   ─── input_directory  ─── input_directory
                        ├── reads0.fastq     ├── barcode01
                        └── reads1.fastq     │   ├── reads0.fastq
                                             │   └── reads1.fastq
                                             ├── barcode02
                                             │   ├── reads0.fastq
                                             │   ├── reads1.fastq
                                             │   └── reads2.fastq
                                             └── barcode03
                                              └── reads0.fastq

Input parameters

Input Options

Sample Options

Output Options

Advanced options

Outputs

Output files may be aggregated including information for all samples or provided per sample. Per-sample files will be prefixed with respective aliases and represented below as {{ alias }}.

Pipeline overview

The following section details the main steps of the workflow.

1. Concatenates input files and generate per read stats.

The fastcat/bamstats tool is used to concatenate multi-file samples to be processed by the workflow. It will also output per read stats including average read lengths and qualities.

2. Stranding and identification of full-length reads

Reads derived from a 10x Genomics library will ideally be flanked by two different adapter sequences. These reads are more likely to represent full length mRNA molecules, although that isn’t guaranteed as some cDNAs may have been created by internal priming or represent other cDNA synthesis artifacts. See this 10x Genomics Technical Note for more information.

This stage of the workflow assigns an adapter configuration to each read. This assigned configurations allows the stranding and orientating the reads. The following schematic shows an example read structure from 10x Genomics 3′ library .

Adapters are located within the reads using vsearch (Read1 and TSO in Fig.1 in the case of the 3prime kit). The table below details the adapter sequences for each of the 10x Genomics kits, along with links to the relevant user guides.

Concatenated reads are identified by adapter configuration and are split into individual subreads and reorientated if required. The following table details the various configurations and the actions taken for each.

Adapter configuration summaries can be found in the output file {{ alias }}/{{ alias }}.config_stats.json"

To only process full length reads the option --full_length_only should be set to true (default: true). If set to false, reads with only a single adapter or other non-full-length adapter configurations will also be processed.

3. Extract cell barcodes and UMIs

The next step is to extract 10x Genomics barcodes and UMI sequences from the stranded and trimmed reads.

In order to do this, the first 100bp of each read are aligned to a reference probe using parasail. This probe contains a suffix of the adapter1 sequence, some ambiguities (“Ns”) representing the barcode and UMI, and a polyT tract.

In this case, the first 16bp of the read that aligns to the N probe region is where the barcode is expected to be, so this sequence is assigned as uncorrected barcode. The following 12 bp are where the UMI is expected, and is assigned to uncorrected umi.

The probe sequence will vary depending on the size of the expected barcode and UMI, which can be found in the following table.

Workflow options:

• The size of the adapter1 suffix can be specified with: barcode_adapter1_suff_length

Once the barcode and UMI has been extracted, these are then trimmed from the reads along with the adapter1 sequence. These trimmed reads are then used in the alignment step.

4. Aligning reads to genome

The next stage is to align the preprocessed reads to the reference genome. This enables gene and transcript read assignment in downstream steps.

The stranded and trimmed FASTQ reads are mapped to the reference genome using minimap2.
The parameter resources_mm2_max_threads int controls the threads given to an alignment process. Other optional parameters can be supplied to minimap2 using resources_mm2_flags (for example --resources_mm2_flags="-I 16GB").

5. Barcode correction

The aim of this stage is to correct errors present in the previously extracted barcodes. 10x Genomics barcodes are not random; all possible barcode sequences can be found in a whitelist of known barcodes. The appropriate whitelist for each kit and version is chosen automatically by the workflow. The whitelist is used to generate a shortlist of high quality barcodes present in the sample, which is then used to correct barcode errors.

The correction proceeds as follows:

• A shortlist of all high quality barcodes present in our sample is generated. This is done by adding an uncorrected barcode to the shortlist if it:
  • has a min quality > barcode_min_quality (default 15)

In each cell library there are expected to be some low quality cells and empty droplets that can be identified by their low number of reads. To remove these cells, the shortlist is filtered with a quantile based threshold. This threshold is determined by ranking the cells by read count and taking the top n cells (n = expected_cells). The read count 95th percentile / 20 is the threshold used. This threshold can be visualised in the knee plots generated by the workflow.

note: As visium barcodes do not represent cells, but rather tissue coordinates, shortlist cell count thresholding is not performed for visium analysis.

For uncorrected barcodes not present in the shortlist, they are cross-referenced against the shortlist, and are assigned a barcode from this list if the candidate barcode meets the following criteria:

• the query and closest-matched shortlist barcode have an edit distance <= 2
• The edit distance difference between the top and hit and second top hit in the shortlist >= 2

note: Edit distance here refers specifically to the Levenshtein distance, which is the minimum number of single base changes it would take to transform one UMI into another (including insertions, substitutions or deletions).

Workflow options:

• barcode_max_ed: Max edit distance between the uncorrected barcode and the matching whitelist barcode (default 2).
• barcode_min_ed_diff: Min difference in edit distance between (1) the uncorrected barcode vs top hit and (2) uncorrected barcode vs runner-up hit (default 2).
• expected_cells: Number of expected cells. Enter the number of expected cells in the sample.

6. Gene and transcript assignment

Now that barcodes and UMIs have been assigned, the next stage is to assign a gene and transcript to each read.

The first step is to create a transcriptome using stringtie2 (using long read mode -L). The minimum required transcripts for a transcript to be called is 2. This can be changed, along with other stringtie settings using the stringtie_opts options (default “-c 2”). The resulting query transcriptome is then annotated with the supplied reference genome annotation using gffcompare, annotating query transcripts with reference transcript and gene metadata.

The original reads are then mapped to this transcriptome, producing one or more alignments per read.

Transcript assignment

To assign a read to a transcript, we follow a similar procedure to that used by FLAMES.

• Transcripts that map to intronic query transcripts (those mapping to query transcripts with gffcompare classes [‘i’, ‘y’, ‘p’, ‘s’]) are set to unknown (-). This does not affect the gene assignment.
• Any read that uniquely maps to the annotated transcriptome is assigned that transcript.
• Any read that does not uniquely map to a transcript is processed with the following procedure:
  • Order alignments by alignment score (AS)
  • If the top alignment has a better AS or query coverage than the second, and has a minimum transcript coverage of 0.4, assign to this transcript, else set as unknown (-)
  • If the AS is the same between the 1st and 2nd alignment the choose the alignment that has the highest transcript coverage, if the coverage is > 0.8

Gene assignment

For each read, the alignment with the highest genomic MAPQ score is used for gene assignemnt if this score > gene_assigns_minqv (default 30). The gene ID is then derived from the transcript assigned in the transcript assignment step.

7. UMI correction

The next step is to correct UMI errors. The previous barcode correction step leveraged a whitelist of all potential barcodes, drastically narrowing the search space. However, UMIs are random sequences and so there is no way of knowing beforehand which UMIs will be in our library. To reduce the UMI search space, the reads are initially clustered by assigned gene, or if no gene is assigned, by genomic interval. This also has the effect of reducing the likelihood of UMI collisions which occur due to the number of total reads in a library usually being significantly higher than the total possible number of UMIS (16.7 million combinations for a 12bp UMI).

UMI clusters are then generated from these gene-clustered reads using UMI-tools with a slightly modified version of the directional method. For this clustering, Levenshtein distance threshold is used (default 2) instead of Hamming distance, in order to account for UMI indels. The selected node of each cluster is the one with the highest number of reads and is used to assign a corrected_umi to the rest of the reads in that cluster.

8. Make expression matrices

The gene x cell and transcript x cell expression matrices are the main outputs of the workflow and can be used in further analysis of the single cell experiment.

The expression matrices are generated by collapsing the corrected UMIs (counting each unique UMI once) and summing the counts of features (gene or transcript) per cell to give a feature x cell expression matrix and are output as a folder of files in Market Exchange (MEX) format (see the output docs

The expression count matrices are further processed in the following way to give gene x cell processed matrices and are also output in Market Exchange (MEX) format.

• Cells are dropped that contain less than matrix_min_genes genes or transcripts (default 200)
• Genes are dropped which are present in fewer than matrix_min_cells (default 3)
• Cells where mitochondrial genes make up more than matrix_max_mito (default 20%) are dropped
• Counts are normalized to matrix_norm_count (default 10,000) reads/cell
• Normalized counts are finally log10 transformed

9. Tagging bam files

BAM files generated from aligning reads to the reference genome are now tagged with the following information from the workflow. By default, the BAMs are output per chromosome, but can be concatenated into a single file per sampe using merge_bam. The new bam will contain the following tags:

• CB: corrected cell barcode sequence
• CR: uncorrected cell barcode sequence
• CY: Phred quality scores of the uncorrected cell barcode sequence
• UB: corrected UMI sequence
• UR: uncorrected UMI sequence
• UY: Phred quality scores of the uncorrected UMI sequence

10. Calculate library saturation

Sequencing saturation is an estimate of the library complexity that has been captured in a sequencing run. As read depth per cell increases, the number of genes and distinct UMIs identified will increase at a rate that is dependent on the complexity of the input library.

Reads which have been assigned a corrected barcode, corrected UMI and gene are subsampled by varying degrees and the resulting median reads/cell is plotted against either:

• median genes per cell: this gives an indication of the gene complexity captured and whether increasing read depth would lead to significantly more genes being identified.
• median UMIs per cell: this gives an indication of the UMI complexity captured. If this plot plateaus out, it may indicate that PCR duplicates represent a high proportion of the reads at higher sequencing depth.
• Sequencing saturation: calculated as 1 - (number of unique UMIs across all cells / number of reads across all cells). Values near 0 indicate that very few duplicate UMIs are being identified, whereas values nearer to 1 indicate a higher proportion of duplicate UMIs are being captured. 1 / (1 - sequencing saturation) can be used can be used to estimate how many additional reads would be required to identify a new UMI. For example, if the sequencing saturation is 0.5 (50%) then for each two new reads, one of those should represent a new UMI. If the sequencing saturation is lower, at 0.2 for example, then on average 1.25 reads would need to be sequenced to obtain a new UMI.

11. Make UMAP plots

UMAP (Uniform Manifold Approximation and Projection) is a data visualisation algorithm used for dimensionality reduction. It aims to preserve the local structure and relationships in high-dimensional data (in this case a gene x cell count matrix) when projecting it into a two-dimensional space. This allows structure within the data, such as cell type and state, to be visualised, and can be useful for quality control and gaining an initial view of the data. To generate UMAP plots use plot_umaps. The data used for the UMAP generation are the processed expression matrices

Several UMAP plots are created:

• gene expression plots show the gene expression UMAPs with each cell annotated with the mean gene expression count
• mitochondrial expression plots are the same but annotated with mean mitochondrial expression
• single gene expression UMAPs again show the gene expression UMAPs, but each plot is overlaid with a gene count data from a single gene. These genes of interest can be specified by adding gene names to a file with the path spcified by umap_plot_genes path/to/gene+names.csv
• transcript expression plots shows the transcript expression UMAP plots with each cell annotated with the mean gene transcript count per cell

The UMAP algorithm is stochastic, therefore analysing the same data multiple times, using identical parameters, can lead to visually different projections. In order to have some confidence in the observed results, it can be useful to run the projection multiple times. The number of repeated projections can be set with umap_n_repeats (default 3)

Troubleshooting

• If the workflow fails please run it with the demo data set to ensure the workflow itself is working. This will help us determine if the issue is related to the environment, input parameters or a bug.
• See how to interpret some common nextflow exit codes here.

When using singularity the following error may occur:

RuntimeError: cannot cache function 'rdist': no locator available for file '/home/epi2melabs/...'

If you receive this error, you may need to set the numba cache directory to a location that is writable by the singularity container. To do this, add the following contents to a file (numba.config for example) and use -c numba.config in the workflow command to apply it.

profiles {
    singularity {
        singularity {
            enabled = true
            autoMounts = true
            runOptions = '--writable-tmpfs'
        }
    }
}
env {
    NUMBA_CACHE_DIR = "${launchDir}/numba_cache"
}

FAQ’s

If your question is not answered here, please report any issues or suggestions on the github issues page or start a discussion on the community.

Related blog posts

Related blog posts

See the EPI2ME website for lots of other resources and blog posts.