What is cellranger doing in comparison to other methods?

Question

I've recently started working with the 10X-Genomics platform with Illumina (MiSeq and HiSeq) for single-cell RNA-Seq. I've been recommended the "cellranger" (version 2.1.0) which I understand handles the barcoding of the platform and performs processing of the data. For example:

1) convert raw Illumina .bcl files into .fastq (demultiplex)

cellranger mkfastq

2) perform alignment, barcode counting, and generate gene-barcode-martices

cellranger count

However, it is unclear which underlying bioinformatics procedure (with which parameters) are being used here. For this pipeline, I wish to know how the trimming, alignment and counts are being performed. This is important to understand whether further normalisation and quality checking is necessary. What steps would I need to do to perform a similar analysis? If for example, I wish to document in more detail what has been performed here in a reproducible manner or reproduce the analysis pipeline on a different platform (such as Nadia) to compare them.

Ideally, I would like to know the analogous pipeline using open-source tools but at the moment I am first seeking to understand what is going on under the hood when I am running the cellranger.

Answer 1

cellranger mkfastq is not necessary anymore. It used to be that the cellranger software wanted the reads to be interleaved, and you could use cellranger to do that for you if you couldn't do it yourself. Newer versions of cellranger will take the fastq files just like Illumina's bcl2fastq makes them.

Cellranger count aligns the reads, filters away duplicates based on UMIs, tries to figure out which cell barcodes really captured cells, and tells you for each cell barcode how many reads hit each gene. You can go through the .bams, and see which reads were counted for what gene and what barcode, based on the flags and other tags in the .bam

If you want another pipeline to compare to, you can use the McCarroll lab's DropSeq pipeline, which is a bunch of little scripts, around the output of STAR aligner (just like 10XGenomics uses). I have put the same data sets through both 10X and Dropseq, and they return about the same thing.

Answer 2

The best approach would be to contact 10x Genomics support: [email protected] . They usually respond within a few hours.

Otherwise, you can check the methods from the first 10x Genomics paper:

The Cell Ranger Single-Cell Software Suite was used to perform sample demultiplexing, barcode processing and single-cell 3′ gene counting (http://software.10xgenomics.com/single-cell/overview/welcome). First, sample demultiplexing was performed based on the 8 bp sample index read to generate FASTQs for the Read1 and Read2 paired-end reads, as well as the 14 bp GemCode barcode. Ten basepair UMI tags were extracted from Read2 (14 libraries were made with 5 bp UMI tags, as noted in Supplementary Table 1, due to an earlier iteration of the methods. For these samples, 5 bp UMI tags were extracted from Read2.). Then, Read1, which contains the cDNA insert, was aligned to an appropriate reference genome using STAR35. For mouse cells, mm10 was used and for human cells, hg19 was used. For samples with mouse and human cell mixtures, the union of hg19 and mm10 were used. For ERCC samples, ERCC reference (https://tools.thermofisher.com/content/sfs/manuals/cms_095047.txt) was used.

Next, GemCode barcodes and UMIs were filtered. All of the known listed of barcodes that are 1-Hamming-distance away from an observed barcode are considered. Then, the posterior probability that the observed barcode was produced by a sequencing error is computed, given the base qualities of the observed barcode and the prior probability of observing the candidate barcode (taken from the overall barcode count distribution). If the posterior probability for any candidate barcode is at least 0.975, then the barcode is corrected to the candidate barcode with the highest posterior probability. If all candidate sequences are equally probable, then the one appearing first by lexical order is picked.

UMIs with sequencing quality score >10 were considered valid if they were not homopolymers. Qual=10 implies 90% base call accuracy. A UMI that is 1-Hamming-distance away from another UMI (with more reads) for the same cell barcode and gene is corrected to the UMI with more reads. This approach is nearly identical to that in Jaitin et al.4, and is similar to that in Klein et al.8 (although Klein et al.8 also used UMIs to resolve multimapped reads, which was not implemented here).

Last, PCR duplicates were marked if two sets of read pairs shared a barcode sequence, a UMI tag, and a gene ID (Ensembl GTFs GRCh37.82, ftp://ftp.ensembl.org/pub/grch37/release-84/gtf/homo_sapiens/Homo_sapiens.GRCh37.82.gtf.gz and GRCm38.84, ftp://ftp.ensembl.org/pub/release-84/gtf/mus_musculus/Mus_musculus.GRCm38.84.gtf.gz, were used). Only confidently mapped (MAPQ=255), non-PCR duplicates with valid barcodes and UMIs were used to generate gene-barcode matrix.

Cell barcodes were determined based on distribution of UMI counts. All top barcodes within the same order of magnitude (>10% of the top nth barcode, where n is 1% of the expected recovered cell count) were considered cell barcodes. Number of reads that provide meaningful information is calculated as the product of four metrics: (1) valid barcodes; (2) valid UMI; (3) associated with a cell barcode; and (4) confidently mapped to exons.

Of course, they have since updated the library prep and the Cell Ranger software, so the exact protocol is somewhat different.