Can alevin-fry be used for data from the BC Rhapsody?

[stela2502]

I mean data coming from this tool: https://www.bdbiosciences.com/en-se/products/instruments/single-cell-multiomics-systems

I assume the fastq reads have a different primer structure. Has anybody worked on that?
You help would very much be appreciated!

[rob-p]

Cc @k3yavi and @jeremymsimon; I recall the generalized geometry support upstream in salmon allowing processing this technology, but I don't recall all of the details.

[rob-p]

Thanks, @stela2502, we'll check this out and figure out what's going on. You van use the view command to see a text dump of the RAD file, but I doubt that will point directly at the issue. We'll debug the backtrace.

[stela2502]

Hi Thank you! I have compiled alevin-fry myself and after uncommenting the line

    println!("{:?}", hdr.ref_names);

I got this info:

["seq_1", "seq_2", "seq_3", "seq_4", "seq_5", "seq_6", "seq_7", "seq_8", "seq_9", "seq_10", "seq_11", "seq_12", "seq_13", "seq_14", "seq_15", "seq_16", "seq_17", "seq_18", "seq_19", "seq_20", "seq_21", "seq_22", "seq_23", "seq_24", "seq_25", "seq_26", "seq_27", "seq_28", "seq_29", "seq_30", "seq_31", "seq_32", "seq_33", "seq_34", "seq_35", "seq_36", "seq_37", "seq_38", "seq_39", "seq_40", "seq_41", "seq_42", "seq_43", "seq_44", "seq_45", "seq_46", "seq_47", "seq_48", "seq_49", "seq_50", "seq_51", "seq_52", "seq_53", "seq_54", "seq_55", "seq_56", "seq_57", "seq_58", "seq_59", "seq_60", "seq_61", "seq_62", "seq_63", "seq_64", "seq_65", "seq_66", "seq_67", "seq_68", "seq_69", "seq_70", "seq_71", "seq_72", "seq_73", "seq_74", "seq_75", "seq_76", "seq_77", "seq_78", "seq_79", "seq_80", "seq_81", "seq_82", "seq_83", "seq_84", "seq_85", "seq_86", "seq_87", "seq_88", "seq_89", "seq_90", "seq_91", "seq_92", "seq_93", "seq_94", "seq_95", "seq_96", "seq_97", "seq_98", "seq_99", "seq_100", "seq_101", "seq_102", "seq_103"
2022-09-27 13:32:27 INFO, "seq_104", "seq_105", "seq_106", "seq_107", "seq_108", "seq_109", "seq_110", "seq_111", "seq_112", "seq_113", "seq_114", "seq_115", "seq_116", "seq_117", "seq_118", "seq_119", "seq_120", "seq_121", "seq_122", "seq_123", "seq_124", "seq_125", "seq_126", "seq_127", "seq_128", "seq_129", "seq_130", "seq_131", "seq_132", "seq_133", "seq_134", "seq_135", "seq_136", "seq_137", "seq_138", "seq_139", "seq_140", "seq_141", "seq_142", "seq_143", "seq_144", "seq_145", "seq_146", "seq_147", "seq_148", "seq_149", "seq_150", "seq_151", "seq_152", "seq_153", "seq_154", "seq_155", "seq_156", "seq_157", "seq_158", "seq_159", "seq_160", "seq_161", "seq_162", "seq_163", "seq_164", "seq_165", "seq_166", "seq_167", "seq_168", "seq_169", "seq_170", "seq_171", "seq_172", "seq_173", "seq_174", "seq_175", "seq_176", "seq_177", "seq_178", "seq_179", "seq_180", "seq_181", "seq_182", "seq_183", "seq_184", "seq_185", "seq_186", "seq_187", "seq_188", "seq_189", "seq_190", "seq_191", "seq_192", "seq_193", "seq_194", "seq_195", "seq_196" paired : false, ref_count : 466, num_chunks : 29,650
, "seq_197", "seq_198", "seq_199", "seq_200", "seq_201", "seq_202", "seq_203", "seq_204", "seq_205", "seq_206", "seq_207", "seq_208", "seq_209", "seq_210", "seq_211", "seq_212", "seq_213", "seq_214", "seq_215", "seq_216", "seq_217", "seq_218", "seq_219", "seq_220", "seq_221", "seq_222", "seq_223", "seq_224", "seq_225", "seq_226", "seq_227", "seq_228", "seq_229", "seq_230", "seq_231", "seq_232", "seq_233", "seq_234", "seq_235", "seq_236", "seq_237", "seq_238", "seq_239", "seq_240", "seq_241", "seq_242", "seq_243", "seq_244", "seq_245", "seq_246", "seq_247", "seq_248", "seq_249", "seq_250", "seq_251", "seq_252", "seq_253", "seq_254", "seq_255", "seq_256", "seq_257", "seq_258", "seq_259", "seq_260", "seq_261", "seq_262", "seq_263", "seq_264", "seq_265", "seq_266", "seq_267", "seq_268", "seq_269", "seq_270", "seq_271", "seq_272", "seq_273", "seq_274", "seq_275", "seq_276", "seq_277", "seq_278", "seq_279", "seq_280", "seq_281", "seq_282", "seq_283", "seq_284", "seq_285", "seq_286", "seq_287", "seq_288", "seq_289", "seq_290", "seq_291", "seq_292", "seq_293", "seq_294", "seq_295", "seq_296", "seq_297", "seq_298", "seq_299", "seq_300", "seq_301", "seq_302", "seq_303", "seq_304", "seq_305", "seq_306", "seq_307", "seq_308", "seq_309", "seq_310", "seq_311", "seq_312", "seq_313", "seq_314", "seq_315", "seq_316", "seq_317", "seq_318", "seq_319", "seq_320", "seq_321", "seq_322", "seq_323", "seq_324", "seq_325", "seq_326", "seq_327", "seq_328", "seq_329", "seq_330", "seq_331", "seq_332", "seq_333", "seq_334", "seq_335", "seq_336", "seq_337", "seq_338", "seq_339", "seq_340", "seq_341", "seq_342", "seq_343", "seq_344", "seq_345", "seq_346", "seq_347", "seq_348", "seq_349", "seq_350", "seq_351", "seq_352", "seq_353", "seq_354", "seq_355", "seq_356", "seq_357", "seq_358", "seq_359", "seq_360", "seq_361", "seq_362", "seq_363", "seq_364", "seq_365", "seq_366", "seq_367", "seq_368", "seq_369", "seq_370", "seq_371", "seq_372", "seq_373", "seq_374", "seq_375", "seq_376", "seq_377", "seq_378", "seq_379", "seq_380", "seq_381", "seq_382", "seq_383", "seq_384", "seq_385", "seq_386", "seq_387", "seq_388", "seq_389", "seq_390", "seq_391", "seq_392", "seq_393", "seq_394", "seq_395", "seq_396", "seq_397", "seq_398", "seq_399", "seq_400", "seq_401", "seq_402", "seq_403", "seq_404", "seq_405", "seq_406", "seq_407", "seq_408", "seq_409", "seq_410", "seq_411", "seq_412", "seq_413", "seq_414", "seq_415", "seq_416", "seq_417", "seq_418", "seq_419", "seq_420", "seq_421", "seq_422", "seq_423", "seq_424", "seq_425", "seq_426", "seq_427", "seq_428", "seq_429", "seq_430", "seq_431", "seq_432", "seq_433", "seq_434", "seq_435", "seq_436", "seq_437", "seq_438", "seq_439", "seq_440", "seq_441", "seq_442", "seq_443", "seq_444", "seq_445", "seq_446", "seq_447", "seq_448", "seq_449", "seq_450", "seq_451", "seq_452", "seq_453", "seq_454", "seq_455", "seq_456", "seq_457", "seq_458", "seq_459", "seq_460", "seq_461", "seq_462", "seq_463", "seq_464", "seq_465", "seq_466"]

Seams as if my 'make up the database' hit an error. If what I assume jas happened, then I should be able to simply change the geneid_to_name.txt file. I'll try that right away.

[stela2502]

YES !!! It is working! Crap that was a hard one :-D - at least for me not having worked with rust at all...

[stela2502]

And I had of cause modified the do_quantify function. Line 398 in the file src/quant.rs...

[stela2502]

Hi rob-p. BD-Rhapsody seams to 'hide' the sample ids from the user. They have to be somewhere in the R1, but I can not look for them as the --bc-geometry seams to be rather picky on what it accepts.
This seams to work, but I can not identify the sample IDs:

--umi-geometry '1[53-60]' --bc-geometry '1[1-9,22-30,44-52]' --read-geometry '2[1-end]'

And neither this "--bc-geometry '1[1-15,16-30,44-52]'" not that --bc-geometry '1[1-52]' will work. Any idea of how to tackle this problem?

Of cause I first need to identify the sample ids in the data...

[rob-p]

Hi @stela2502,

I'm glad you got it to at least run; that's progress. I think one issue with that barcode geometry is that the resulting barcodes will be very long If you do e.g. 1[1-52] or some such. Right now, it will be difficult to deal with barcodes > 32 nucleotides long in total. We do have a path to support those, but it's already a substantial ceiling above what most common protocols that I'm aware of use.

My apologies for not being conversant in all of the different protocols, but is there a document that we can refer to regarding the barcode layout for BD-Rhapsody?

Thanks!
Rob

[stela2502]

I have contacted BD support and hope they can provide that in the future.
From what I have seen they might consider this a company secret. At least they do not seam to be very open with this.
I have not looked deeper into the program they provide to analyze the data as I have better things to do right now.
I really hope they come up with the required info.

[rob-p]

Me too! Perhaps @k3yavi and I can also ask around to folks we know using this tech to see if they have any further insight.

[stela2502]

Hi, sorry for the late answer. We have analyzed the same data using the 7 bridges 'tool'. But this does not even report the sequence behind the sample ids (just integers). So I assume it is a 'company secret'. I have not been in contact with the BD people, but I might try that tomorrow.
So now I am trying to decode the read and check whether I somehow can find 6 samples in our fastq files.
I assume that there are some smaller barcodes hidden in the areas 1[10-21,31-43] I am digging into this right now and if I find the data I actually do not necessarily have support in alevin-fry for that. If need be I can run a perl script and extract all this info one more times, export it as csv and read this info in later on. The issue is that the cell barcode would likely be merged to the sample barcode if alevin-fry could handle these long barcodes. So in a way the external approach might be more straight forward. And as the 7 bridges program takes about a day to finish and alevin-fry more like 10 minutes we have plenty of time to waste on finding the sample ids.

[stela2502]

Hi Rob,
from what I read about the BD sample IDs the R2 is mainly sample ID. It also contains a part of the primers, but that can easily be cut away. Together with the list of maximal 12 samples I have written a small R script that processes R2 and uses the Text::LevenshteinXS qw(distance); function to find the 'closest' sample ID. At the moment the distance info is not used in any way.
To use this with alevin-fry I would need to split the fastqs before running the mapping - right? Are you planning to add an internal sample splitting to alevin-fry to support natively BD_Rhapsody in the near future? Otherwise I go on splitting the fastqs ;-)

Thank you for being so responsive! Really great help I got from you up to now!

[stela2502]

extractCellBarcodes_from_R2.zip
The pl script was not supported by github ...

[stela2502]

Hello and sorry for the long wait. I finally could identify the most problematic part - to identify the samples.
At the moment alevin-fry can not identify the sample based on the R2 read. And I think that would likely be overkill anyhow.
What I thought might be a way out is to open both the R1 and R2 files and split both on the sample id in R2.

I implemented that in a very crude way using Perl. And it does take forever.
Given the fact that alevin-fry is extremely good in mapping fastq files to a genome/transcriptome, would it be possible to use it or the downstream logic to split two fastq files into sample specific fastq files?
Or do you know of any other efficient way to doing that?

Thank you for your help so far!

[stela2502]

Would likely make sense to to move the discussion there?
Thank you for your help! I will also try to ask the needletail people for help.

[stela2502]

onecodex/needletail#60

[stela2502]

Hi Rob. Now I 'fixed' the N problem by simply not using any [u8] that contains them....

And after fixing an internal problem with searching the wrong file I finally got it to work.
Quite crucial in the end was a try_finish call on the GzEncoder objects.

THANK YOU ONCE MORE!

I can't wait to apply this to a real dataset! :-D

[rob-p]

That's awesome! The culprit with the N's is likely our kmer library, which expects the strings to be sanitized (i.e. Ns) removed before it packs them in 2-bit encoding. I think the current solution is fine as a quick hack, the other solution would be to just replace the N characters with some other base either deterministically, or pseudorandomly.

I'm happy to help and, as I said, would be happy to help improve and refine the program as you encounter issues or as there is more data to evaluate with :).

[stela2502]

Yes I think you are right with the kmer library. But for the programs logic I can simply remove these kmers. They would not match anyhow as there are no N's in the sample ids :-D. So in the end it should make the program even faster :-D.

I have changed some other problems. E.g. rust did not like two command line options starting with the same letter (read_file1 and read_file2). I remember you had a lock on the out stream in splitp to speed the process up. Is that something that could also be applied to the GzEncoder?

At the moment it looks like there are a lot of R2 reads that do not even match to more that 2 9bp kmers of the sample ids. So there seams to be something off with BD rhapsody's data quality or the whole chemistry logic. But that is a little early to draw this conclusion.

[stela2502]

Another question: As the split2sample mainly decodes a gz file to then encode the same info and save it into another gz file - would there be an option to 'collect' all the gzipped data for a fastq entry and then write that into the output stream? Like an additional variable in a fastq entry - say "gz" that just harbors the gzipped whole fastq entry.
I am asking as just to gzip a 89Gb big fastq file took almost an our on our cluster. If I think that the split2sequence would unzip twice that info just to re-zip it a little later. Who knows - probably we could save two hours calculation time?

[rob-p]

Hi @stela2502,

I just wanted to say I ran across a BC Rhapsody-related tweet from David Eccles, another salmon/alevin user and it occurred to me that he had looked at this a while back in another context. I reached out to him on Twitter and he replied with what look to be some nice resources (https://twitter.com/gringene_bio/status/1585339519532163076?s=20&t=HNah1vFHPXv1gvESm8E7qA)! Perhaps we can see what functionality exists and what duplicate effort can be avoided.

Best,
Rob

[gringer]

In addition to the cell squish code I mention here, last year BD provided me with python code including the cell index lookups. I converted it to R and added an inverse lookup to allow both forward and reverse lookups from the cell barcode sequence to the ID used by the SevenBridges workflow (see attached, and now included in my bioinfscripts repository here).

These barcode sequences seem to be identical with both the short and long barcode formats that I mentioned in the salmon issue, so I expect the ID/sequence mapping will remain the same.

Rhapsody_cell_keys.r.txt

[stela2502]

Thank you very much! This makes sense and I can think of how to use this. Thank you!
@rob-p I assume the sample identification needs to become a two step approach. First step would be to collect cell ids per sample and the second to split the fastqs into sample specific ones using the sample->cellID info.
With this R script and your splitp program I should be able to convert the R1 sequence into a numeric cellid. Would it make sense to modify the R1 for alevin_fry? I mean during the split we could fix 1bp differences to the original cell ids by simply re-creating the R1 from the numeric cell ID :-D (and the UMI). At least fix the parts where the cell id is stored.
That will be a lot of fun to make this work!

[gringer]

Yeah, removing the 0-3bp variable-length prefix in R1 seems to help a lot for salmon alevin: my mapping rates have doubled.

For what it's worth, I've discovered that the prefixes are attached to specific BC1 barcodes, rather than being randomly distributed across all reads. I haven't added those prefixes into my barcode list, because I couldn't work out a quick way to easily incorporate variable-length barcodes into my workflow. It was easier to trim off and discard the prefix, and feed the result into my existing workflow.