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: https://stackoverflow.com/questions/52408588/parallelizing-nested-for-loop-with-gnu-parallel
Thanks for your suggestions!
The
GitHublink appears to be incorrect.Thanks, link updated.
Is your script generally applicable/easily extensible to other SPLiT-seq datasets or is it for this specific one you mention above?
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.
Thanks! You may want to mention that in original post/on GitHub.
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 (bcl2fastqstarts from .bcl), not sure how much divergence in functionality there is after that step.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.
To flesh this out more as a fully fledged tool, I would suggest showing some run time comparisons for different data sizes etc.
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.
The 8bit barcode sequences used in the SPLiT-Seq protocol are redundant between the different barcoding Rounds. You can confirm this here https://drive.google.com/file/d/1l3MkxYVWOPTV5KPpq2uXMLzKOKwVeZyw/view?usp=sharing
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.
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.
@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 "splitseqdemultiplexing.sh" 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 "splitseqdemultiplexing.sh" 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.
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.