rob-p
commented
3 years ago Hi @jeremymsimon,
Somehow, the notification for this in my e-mail got classified as SPAM. Anyway, thank you for the detailed description! I'm going to ping @Gaura here. @Gaura β this is the alternative protocol I was discussing with you yesterday. As you can see, the main issue here is the "noisy" barcodes. Let me know what you think would be necessary to add support for this, and I'm happy to schedule a technical discussion if you want to discuss some options.
jeremymsimon
Gaura
Hey @rob-p and @k3yavi! Avi and I have chatted a bit about this offline, and I've commented about this in the past, but I'm currently attempting to extend alevin's use to the newer combinatorial indexing experiments like SPLiT-seq or the commercial solution, ParseBio (the biotech formerly known as SplitBio). Now that you support different barcode geometries and allow the user to specify where these barcodes should be expected, this should be doable.
I've downloaded the raw data from the original SPLiT-seq paper, and we have some ParseBio data of our own; my goal is to help you all write a vignette on this once we get the details worked out.
In short, the R1 file contains cDNA sequence, and R2 contains a 10bp UMI followed by 3 different barcodes in the following arrangements (note: ParseBio seems to shift BC1 up a few bp, to 79-86):
The intervening bases are a combination of common and variable sequences.
The other wrinkle with these approaches is that the cDNA can be amplified using either oligo-dT or random hexamer barcodes. The net effect of this is that BC1 can be one of two options, and some barcode combinations need to be pooled to represent the same cell. For example:
AACGTGAT-CTGTAGCC-ACACAGAA GATAGACA-CTGTAGCC-ACACAGAA
may be the same cell, amplified by two different means.
We know how the BC1s are paired, so Avi suggested pre-processing the R2 FASTQ file such that we locate random hexamer BC1 sequences and modify them to the matched oligo-dT sequence. In the above example, we'd correct some R2 reads that contain "GATAGACA" in the 87-94 position to "AACGTGAT". I wrote a script (in perl because I'm old-school) that does this and allows for mismatches within 1 hamming distance. I'm sure someone could port this to python, rust, or C/C++ and speed this step up substantially, but it seems to work fine.
The next step is to run alevin. I've done so using salmon v1.5.2 with the following parameters (for a ParseBio run), reversing R1/R2 as per Avi's suggestion:
I know that these libraries are stranded, but I doubt
ISRis the correct architecture, especially given it's not true PE sequencing. Notably, running alevin this way gets me only a 4% mapping rate, and most of my cells get filtered out downstream due to having really low depth.Below is my
salmon_quant.log:And here is my
lib_format_counts.json, which has 0s for the different potential libTypes and shows a very high rate of inconsistent/orphan mappings:And here is the
alevin_meta_info.json, which shows the huge attrition in recognized CBs and poor mapping rate:This, in conjunction with the fact that running with
-l Asuggests incorrectly the data are unstranded (IU), led @mikelove to suspect alevin may be skipping the libType guesswork that salmon typically performs.I know these data are of reasonably good quality, because I have instead processed them using zUMIs, which supports SPLiT-seq data. A side-by-side of the same FASTQs processed this way using alevin or using zUMIs gave me ~270 cells that pass filters (alevin) vs ~50k cells that pass filters (zUMIs), so something is definitely up here.
I'd like to confirm a few details and ask for some guidance on how to move forward:
I've also attached here the top 1000 lines of each R1/R2-corrected FASTQ (ParseBio) in the hopes that was somehow helpful in diagnosis, but I can share more reads or the full files some other way if needed
head_R1.fastq.txt head_R2_corrected.fastq.txt
Thanks so much!
Jeremy