scAN1.0 : A reproducible and standardized pipeline for processing 10X single cell RNAseq data

2.1 Version annotation effect

We first assessed the impact of the Ensembl version of the GTF file on the final output using the human Gonadotroph tumour dataset as an input (see Material and methods). GTF files and their corresponding FASTA files were downloaded via ftp protocol using the –version parameter (see section 4). We set Cellranger as the default mapper and then assessed the impact of 4 different annotation releases (93, 98, 103 and 106) on the number of detected genes (Fig. 2A) and on the count per genes as assessed with the CHGA gene (Fig. 2B).

OPEN IN VIEWER

Fig. 2

Impact of the GTF version. The human gonadotroph tumors dataset was submitted to scAN1.0 with 4 different versions of the GTF for mapping with Cellranger. A. Venn diagram of the number of the detected genes when using the different versions. B. Violin plot representation of the impact of the GTF version on the UMI counts for the CHGA gene.

The overall impact of the GTF version seems relatively modest, especially in regard to the number of UMI counts for CHGA. However, the 106 version allowed to identify a larger number of genes and was kept for the next step.

We then assessed the impact of filtrating the GTF file with the cellranger mkgtf function. This filtration step is intended to remove unwanted genes classified by biotype, such as long non-coding RNAs (https://www.ensembl.org/info/genome/genebuild/biotypes.html). We used the default values of that function that removes biotypes such as gene_biotype:pseudogenes from the GTF annotation file. This step reduces the presence of multi-mapping due to non-regularly processed nucleic acids, removing ambiguity about gene identification and allowing quantification of properly mapped sequences. We observed that this filtration step indeed had a major impact on both the number of genes detected (Fig. 3A) but also on gene counts (Fig. 3B).

OPEN IN VIEWER

Fig. 3

Impact of the GTF filtration. A. Venn diagramm of the number of the detected genes when using either an unfiltered (106) of filtered (filter106) version of the GTF file. B. Violin plot representation of the filtration impact on the UMI counts for the CD68 gene.

We then compared the impact of two popular alignment tools for single-cell RNA sequencing (Kallisto-bustools and Cellranger) using the 10x Genomics pre-built Cellranger reference packages version 2020-A for human. As seen in Fig. 4A, 81% of the genes were identified by both algorithms whereas Kallisto-bustools identified more genes than Cellranger.

OPEN IN VIEWER

Fig. 4

Impact of the mapping tool. A. Venn diagram of the number of the detected genes when using either Cellranger (CR) or Kallisto-bustools (KB) as an alignment tool. B. Violin plot representation on the UMI counts for two genes, CHGA and RBP4.

The impact of the mapper on counts for specific genes seemed to be negligible (Fig. 4B). Therefore this tends to favor Kallisto-bustools for downstream analyses.

We also assessed the impact of the mapper choice on the final clustering step. As seen in Fig. 5, the impact was modest but apparent (e.g. cluster number 1 in the Cellranger dataset was split in two in the Kallisto-Bustools dataset).

OPEN IN VIEWER

Fig. 5

Cluster characterization: UMAP representation (A and B) and Silhouette scores (C and D) of the clusters obtained on data processed with CellRanger (A and C) or Kallisto-bustools (B and D). In E, is shown an alluvial plot highlighting the conservation and differences in cluster composition depending upon the initial mapping method.

Please, note that the different positioning of the clusters on each UMAP is due to known flaws associated with a 2D representation using UMAP (see e.g. [3])

The final combination that was found to be the most effective for our human dataset therefore was using Kallisto-bustools together with the filtered version of the 106 GTF.

We decided to apply scAN1.0 to a published scRNA-seq dataset acquired from murine individual CD8+ T cells throughout the course of their differentiation in response to viral infection [4] and in order to compare the two alignment tools Kallisto-bustools and Cellranger.

As seen in Fig. 6, overall less genes were detected than observed in human cells with 11587 genes identified by both algorithms. Quite interestingly we observed the very same increase (plus 20%) in in the number of genes identified by Kallisto-bustools when compared with Cellranger.

OPEN IN VIEWER

Fig. 6

Murine dataset analysis: Venn diagramm of the number of genes detected using Cellranger (CR) or Kallisto-bustools (KB) as alignment tool the murine CD8 dataset.

4.1 Pipeline description

Figure 1 describes the overall processing of the sequencing files with the ordering of all steps described in subsection 4.5.

4.2 Gonadotroph tumour dataset

Single cell preparation and sequencing. A tumor fragment from a gonadotroph surgically-resected adenoma was collected in Dulbecco’s Modified Eagle Medium (DMEM, cat 41965062; Life Technologies). Single-cell suspension of the resected fragment was obtained through mechanical and enzymatic dissociation (Collagenase P, cat 11213865001) followed by filtration through a 100 μm mesh-strainer (#732-2759, VWR international).

Red blood cells were eliminated using a 10-minute incubation with a commercial red blood cell lysis buffer (eBioscience, cat #00-4300-54). The whole process was achieved within the 2 hours following the surgical resection, cell viability was evaluated to reach at least 70 percent prior encapsulation. Generation of the library was done using a Chromium controller from 10xGenomics. The entire procedure was achieved as recommended by the manufacturer’s for the v3 reagent kit (Chromium Next GEM Single Cell 3’ Kit v3.1 (10xGenomics cat#1000269). Single cell suspension was loaded onto a Chromium Single Cell A Chip, aiming for 5,000 cells in total. The cDNA was amplified after a reverse transcription step, followed by a cleaning using SPRIselect (Beckman Coulter), a quantification step, and enzymatic fragmentation before being sequenced using a NextSeq500 system (Illumina). FASTQ files were generated from sequencer’s base call files (BCLs) for each flowcell directory using cellranger mkfastq.

Study ethic approval. This work is part of the SPACE-PIT study (MR004 n21-5439). It was approved by the Hospices Civils de Lyon ethical committee and registered at the “Centre National Information et Liberté” (CNIL.fr) under the reference 20_098. Informed consent was obtained from the patients.

Processing. The dataset was processed using scAN1.0 with the following parameters:

species = human (default option)

The following options were used to choose either KB or CR:

version = 106

The following options were used for filtrating GTF:

filtergtf

The following options were used to choose GTF version:

version = 106

The number of principal components of the PCA used for UMAP embedding and clustering was set to 10 as determined by the rule of thumb heuristic and the broken stick method [10]. The clustering resolution was set to 0.7.

4.3 CD8 T cells dataset

Single cells preparation and sequencing. We downloaded reads from GEO (GSE131847, https://www-ncbi-nlm-nih-gov.proxy.insermbiblio.inist.fr/geo/) with fastqdump (https://rnnh.github.io/bioinfo-notebook/docs/fastq-dump.html) using the –split-files option. Quality control of reads was performed with fastp (https://github.com/OpenGene/fastp) and resulting FASTQ files were submitted to scAN1.0.

To briefly describe this dataset (see [4] for the full description), P14 transgenic CD8 T cell, which recognise an LCMV epitope, were transferred to histo-compatible hosts that were acutely immunised the day after with the virus. Six days after infection, single responding P14 CD8 T cells from the spleens of immunised hosts were sorted and loaded into Single Cell A chips for partition into Gel Bead In-Emulsions in a Chromium Controller (10x Genomics). Single-cell RNA libraries were prepared according to the 10x Genomics Chromium Single Cell 3 Reagent Kits v2 User Guide and sequenced (paired-end) on a HiSeq 4000 (Illumina).

Processing. The dataset was processed using scAN1.0 with the following parameters:

version = 106

The following options were used to choose either KB or CR:

quantif = cellranger or quantif = kb

4.4 Implementation

The scAN1.0 pipeline was implemented using the reactive workflow manager Nextflow [2], coding with the DSL2 syntax extension. Nextflow simplifies the writing of computational pipelines by making them portable, scalable,parallelizable and ensuring a high level of reproducibility. Nextflow provides native support for container technologies such as Docker or Singularity. Each process in the pipeline will be run in a container. A reproducible container environment is built for each process from Docker images stored on the DockerHub. The analyses can be run on the user’s preferred computing platform. Using the configuration file and corresponding profile, the pipeline can be run on a local computer via Docker or Singularity, as well as on a high-performance computing (HPC) cluster or in cloud-based environments. In the case the workflow cannot complete all steps, a cache-based pipeline resume feature allows to recover the processed steps and the blocked process for debugging before resuming the workflow.

4.5 Input

scAN1.0 takes three mandatory parameters as input: the paired-end FASTQ files from the 10X Chromium sequencing and two genomics files (one FASTA and one GTF).

FASTQ files store the nucleotide sequence and the associated sequencing quality scores. For paired-end sequencing, a FASTQ file is provided in a compressed gzip format for each read (R1 and R2). In the case a sample has been sequenced on several lanes, all R1 and R2 reads are concatenated together.

The mapping steps requires the input of two additional files corresponding to the species of interest:

A FASTA genomic file which stores the raw genome sequence.

4.6 Preprocessing & Mapping

First, FASTQ files are processed and trimmed by using Fastp v0.20.1, an ultra-fast FASTQ preprocessor with useful quality control and data-filtering features [11]. Reads with phred quality (https://en.wikipedia.org/wiki/Phred_quality_score) at least equal to 30 (-q 30) are qualified for the quantification step. Length filtering is disabled (-L) while adapter sequence auto-detection is enabled (–detect_adapter_for_pe). To reduce overlapping annotation, we recommend and implemented (parameter –filtergtf) an optionnal tool of Cellranger (mkgtf function) that allows the GTF filtration based on biotype attributes of sequences. (see: https://support.10xgenomics.com/single-cell-gene-expression/software/release-notes/build and https://www.ensembl.org/info/genome/genebuild/biotypes.html)

The processed files are then mapped to a reference genome in order to quantify gene expression. In scAN1.0, the user can specify two different mappers (see below): Kallisto-bustools v0.26.0 [9] or Cellranger v5.0.1 [12].

Kallisto-bustools is used thought the Python wrapper: kb_python. Starting with FASTA and GTF files, an index of the reference can be built as a colored De Bruij graph with Kallisto via kb ref. with default parameter. Once an index has been generated or downloaded, kb count uses Kallisto to pseudoalign reads and bustools to quantify the data.

4.7 Quality control

Empty Droplets. With droplet-based technologies, most of the barcodes in the matrix correspond to empty droplets (e.g., barcodes with the sum of expressions over all genes being null). They must be removed from Kallisto-bustools generated gene expression matrices. Thus, the Kallisto-bustools outputs were imported into R with a customized R function. The UMI total counts were ranked using barcodeRanks() function from DropletUtils v1.14.2. Empty droplets were removed by selecting the inflection point value on the resulting knee plot (lower cutoff = 10).

On its side, Cellranger handles empty droplets by itself and matrices cleared of empty droplets were directly imported from the standard filtered barcode output.

After importation, either gene expression matrices were converted as Seurat object (Seurat v4.0.4) with CreateSeuratObject function, including features detected in at least 3 cells (min.cells =3).

Low quality cells. On the resulting empty droplet-free matrices, the pipeline computes 3 QC metrics per sample: the number of unique features, the number of UMIs and the percent of mitochondrial gene counts; per cell. These QC metrics are then used to discard three main types of low quality cells [13]:

Cells in apoptosis may exhibit high% mitochondrial and low number of UMIs per cell.

For this, a threshold can be set by the user on scAN1.0 parameters min_feature_RNA, min_ncount_RNA, max_percent_mito, max_feature_RNA, max_ncount_RNA; which default values are 500 and 0 for the first two and “adaptive” for the others. With the “adaptive” configuration, scAN1.0 defines a certain number of median absolute deviations (MADs) away from the median to define maximum values (median + 3*mad) [14].

Furthermore, the pipeline uses the R package DoubletFinder v2.0.3 to detect and remove the potential doublets from the dataset. See [15] and [16].

Non-expressed genes. Genes with sum count along cells equal to 0 (e.g. not-expressed genes) are removed.

4.8 Normalization

After selection of high quality single-cells droplets and significantly expressed genes, the count matrices are then normalized using Sanity, a Bayesian algorithm to infer gene-expression state [17]. We are fully aware that the normalization of single cell transcriptomic data is a research field on its own (see e.g. [18] and citations therein), and we expect this block in the pipeline to be susceptible to be modified in future versions of the pipeline. The modularity of the Nextflow syntax makes it ideal for such additions.

4.9 Clustering and two-dimensional visualization

A final step consists in variable features selection with Seurat::FindVariableFeatures using the vst method (selection.method = “vst") and selecting the 2000 first highly variable genes (nfeatures = 2000), followed by a first linear dimensionality reduction using PCA (Seurat::RunPCA with default parameter). The m first axis of the PCA are then used for the nearest-neighbor graph construction with Seurat::FindNeighbors function (dims=1:m). Cluster determination was performed using the Louvain algorithm [19] run with Seurat::FindClusters function. The resolution parameter set by default to 0.7. The quality of the clustering was assessed using the Silhouette score [20]. Finally, non-linear dimensional representation (using either t-SNE [21] or UMAP [22]) is then performed using Seurat::RunTSNE or Seurat::RunUMAP with default parameters using the same number of dimensions than the nearest-neighbor graph building.

The resolution parameters used for clustering and the number of principal components kept for both clustering and dimension reduction embeddings can be respectively modified by the user at the start of the pipeline by the –resolution_clustering and the –principal_component parameters (default values: 0.7 and 10, respectively).

Alternatively, this last step can be skipped allowing the user to use their own clustering method. Similarly, to avoid the introduction of layers of complexity and simplify the pipeline usage, the automatic annotation of clusters was not introduced. Users can annotate their dataset manually.

4.10 Output

The complete list of process outputs is available on the Readme file of the gitlab repository (see section 5).

By default, all outputs are saved in the results directory, which is structured into subdirectories that are named according to each process. Results and graphics are stored in the subdirectories that correspond to each process. For the post-quantification steps, scAN1.0 provides a Seurat object and an updated matrix for each step. The user can easily manipulate the Seurat object by loading it into R. Since cluster annotation has not been integrated into the pipeline, users can manually annotate the clusters using marker genes or cell type reference signatures.