Tool: SPLiT-Seq Demultiplexing a bash tool for extraction of barcoded single cells
gravatar for paulranum11
9 months ago by
paulranum1130 wrote:

Hi Everyone,

I am interested in the recently published SPLiT-Seq single cell/nucleus RNA-Seq paper from Rosenberg et al that came out in Science

In this paper the authors describe a custom iterative barcoding strategy to index single cells for single cell RNA-Seq without droplets or cell isolation. They have archived their data on GEO: GSE110823. However I was unable to find any scripts or bioinformatics tools that they used to process the raw data.

I have attempted to throw together a SPLiT-Seq demultiplexing strategy using bash. With the goal of making it more robust, accurate, and functional I wanted to share it with the bioinformatics community. Any comments, criticisms, bug reports, or pointers to any official SPLiT-Seq tools are much appreciated.

Here is a link to the Github repository.

sequencing rna-seq tool • 1.5k views
ADD COMMENTlink modified 6 months ago by rainwashautumn0 • written 9 months ago by paulranum1130

To flesh this out more as a fully fledged tool, I would suggest showing some run time comparisons for different data sizes etc.

ADD REPLYlink written 9 months ago by jrj.healey13k

The GitHub link appears to be incorrect.

ADD REPLYlink written 9 months ago by genomax69k

Thanks, link updated.

ADD REPLYlink written 9 months ago by paulranum1130

Is your script generally applicable/easily extensible to other SPLiT-seq datasets or is it for this specific one you mention above?

ADD REPLYlink written 9 months ago by genomax69k

It should be directly applicable to any SPLiT-Seq dataset using the published indexes. If new indexes are used a user can simply change the barcodes.txt files to reflect their index sequences.

ADD REPLYlink written 9 months ago by paulranum1130

Thanks! You may want to mention that in original post/on GitHub.

ADD REPLYlink written 9 months ago by genomax69k

Can you describe how the function of this is different from the standard bcl2fastq? From the looks of your script, you are starting with .fastq files (bcl2fastq starts from .bcl), not sure how much divergence in functionality there is after that step.

ADD REPLYlink written 9 months ago by steve2.0k

I am not an expert on bcl2fastq but I understand that it is used as the first demultiplexing step after Illumina sequencing taking the .bcl files generated by the sequencer then both demultiplexing and converting them to standard .fastq format for each indexed sample. To do this bcl2fastq is capable extracting samples indexed with standard Illumina barcodes from one of their index kits. Because SPLiT-Seq uses a custom 3 stage ligation based barcoding design in addition to standard illumina indexes I assume that bcl2fastq is not capable of demultiplexing the single cells based on this non-illumina indexing strategy. As an example the archived SPLiT-Seq datasets on GEO appear to have had the illumina indexes removed while the split-seq barcodes 1-3 are still there.

I am sure there are similarities as both tools try to demulitplex samples... so I would be very interested to know if bcl2fastq could be made to work on SPLiT-Seq data. Perhaps other people here may have more insight.

ADD REPLYlink modified 9 months ago • written 9 months ago by paulranum1130

Hello, firstly I would like to appreciate for your work. But there is some doubts about the Round*_barcodes_new.txt files. I am curious about why round1 barcode is 13 bits and 2 is 13 bits, 3 is 14or 15 bits. We know that the round1 unique barcode is 8bits, yet you add 5 redundant bits which is totally the same in all barcodes.

ADD REPLYlink written 6 months ago by rainwashautumn0

The 8bit barcode sequences used in the SPLiT-Seq protocol are redundant between the different barcoding Rounds. You can confirm this here

The barcodes can be redundant because they have flanking sequences used for the ligation that are distinct. As a result the position of each barcode and its flanking sequence in the resulting read is unique. The several bases of flanking sequence added to the 8bit barcode are included in the "barcodes_new3.txt" files to allow the script to detect a match unique to the intended barcode position. For example, a Round2 barcode will not detect a barcode added during Round 1 or Round 3 because they lack the Round2 flanking sequence. If you look at an example read sequence and the barcode sequences that I linked to above this logic will become clear.

ADD REPLYlink modified 6 months ago • written 6 months ago by paulranum1130

Thanks a lot for your meticulously reply. Now I have noticed that the full dataset which is used for benchmarking contains only 131 nuclei (according to the paper )while your script generated 14,423 .fastq files (cells) as output. I am new to bioinformatics and hope my question do not bother you a lot.

ADD REPLYlink written 5 months ago by rainwashautumn0

@rainwashautumn I have also noticed this discrepancy with GSM3017260 (SRR6750041). I emailed the authors of the paper months ago about this point but never received a response. I think the issue may be one of two things 1. the dataset is mislabeled in the GEO archive or 2. They performed additional filtering to somehow arrive at 131 nuclei... eliminating a very large number of barcode combinations... From the testing we have done on "" we know that it performs as expected on datasets containing a computationally predefined number of barcode combinations. We also have confidence that each of the barcode combinations identified exist in the .fastq file downloaded from SRR6750041 because "" records each combination found and we can confirm its presence in the raw reads.
Finally, we seem to get results much closer to those reported in the Split-seq paper if we process other files for example GSM3017260 (SRR6750042). Anyway you point out a very interesting discrepancy that I have not been able to wrap my head around. If you figure it out please share your findings with us here.

ADD REPLYlink written 5 months ago by paulranum1130

HUGE UPDATE to SPLiT-Seq_Demultiplexing released on Jan_16_2019. An enormous increase in speed was achieved through rewriting steps 1 and 2 in python to make use of hash data structures. Support for a user defined number of mismatches was added. Suggestions, comments, and criticisms welcome.

ADD REPLYlink modified 6 months ago • written 6 months ago by paulranum1130
gravatar for RamRS
9 months ago by
Houston, TX
RamRS22k wrote:


I have a few questions and some feedback on the code:


  • What is a "round"?
  • Why are there 3 rounds?
  • Is that the typical use case or can be there more (or fewer) "rounds"? (I see that the paper answers these questions)
  • Why use cut -b 1- to read the files into an array?

Code feedback:

  • I'd recommend switching to files passed as command line arguments and processed using getopts, instead of hard-coding filenames in scripts. That way, you can have your input file wherever
  • Use set -eo pipefail so the script stops if any step fails
  • Ensure log and out files can be created (i.e. dir has write permissions), or better yet, just write to STDOUT and STDERR as appropriate
  • Do not implicitly rm. Instead, create directories of the format prefix.XXX using mktemp [-d]. See here
  • There are a lot of loops and conditionals, maybe restructure the code so you need fewer loops? Better yet, switch to python. Shell loops are bad (Advanced topic)
  • Unless you're using regular expressions, grep -F (or grep -Fx if you're matching an entire line) will reduce the load on grep.
  • You're exporting and parallelizing the grep function, why not do that all over the script instead of just in one place?
ADD COMMENTlink modified 9 months ago • written 9 months ago by RamRS22k

Ram, thanks very much for your feedback I appreciate it.

Explanation of "rounds": The SPLiT-Seq science paper describes 3 successive "rounds" of barcoding. I am borrowing the round1, round2, round3, terminology from their publication. Each round of barcoding done in the wet lab is performed through the addition of a specific set of barcode Oligos. The typical use case is 3 rounds of barcodes. Three rounds of barcoding are used to achieve a high number of possible barcode combinations. Round1 (48) x round2 (96) x round3 (96) = 442,368 possible combinations.

Code feedback: Thanks for your detailed comments. I will work on revising the code.

One area that you identify for revision (avoiding shell loops) is an area I have looked at but struggled to revise particularly in the step1 (demultiplexing) part of the script. I created a stack overflow question on this topic here:

Thanks for your suggestions!

ADD REPLYlink modified 9 months ago • written 9 months ago by paulranum1130

And you got a reply from the author of GNU parallel ;)

Ole won’t lead you wrong!

ADD REPLYlink written 9 months ago by jrj.healey13k

The person on SO has parallelized the innermost loop, I'm guessing you could try applying that logic (abstracting the logic as a function, using input files are arguments and parallelizing the function) to the outer loops as well, although I'm not sure how the task distribution would work in such a implementation.

ADD REPLYlink written 9 months ago by RamRS22k
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1738 users visited in the last hour