How do I increase the sensitivity of Salmon Alevin?

Question

I'm in a little bit of a bind with targeted single-cell sequencing. I'm trying to match up our reads to the targeted amplicon panel (418 targets), and all but one have matched successfully with Salmon Alevin. I have a transcript that produces many 75bp reads that only match the first 28bp of a target, and Alevin is ignoring those reads as a match (the 28bp sequence is genome-unique).

Here's an example read, with matched sequence:

qual FFFFFFFFFFFFFFFFFFFFFFFFFFFF,,,::,,:,:,,,,,,,,,,,,,,,,,:,,,,,,,,:,,,,,:,:
read GCTGCATAATCTCTTCATTCCGAGGAGCACCCCACGCGCGCCGCGCCCCCCCGCCCCCGGGGGGGGGGGGGGG
SNPs ...................................G.G.G..G.GC.....CG...CCGGGGGGGGG..GG..
 ref GCTGCATAATCTCTTCATTCCGAGGAGCACCCCACCCCCACCCCCACCCCCACCCCATTCTTAAATTGTTTGG

As you may be able to see from the above sequence, the reads for this particular transcript get up to 28bp (sometimes with one or two errors), then the wheels start coming off the sequencer, showing an abundance of Gs. This was a NovaSeq run, and a polyG read is the same as no signal using its two-dye process, so I suspect that the sequencing stops after ~35bp, rather than actually synthesising Gs from thin air.

I have lots of these reads; there are about 400k reads that match the reference transcript (Cxcr4, just in case it matters), and about 250k of those reads have a G homopolymer that's at least 8bp in length.

Any ideas on what to do?

FWIW, the above match was taken from a Bowtie2 alignment I did on the reads. I could use Alevin for cell barcode correction only and do the UMI deduplication manually, but it seems a shame to do that when it's only one transcript (from 418) that's spoiling a fully-automated process.

I did try adding in an additional transcript that included the full sequence of a short-match read, but there seem to be too many errors in the non-matching sequence for any reads to count (at least... I assume that's why the counts are not coming through; I'm not quite sure how I can work out what's going on behind the scenes).

Is there any way to tweak Salmon Alevin to get it to accept these shorter matches? If I can change the kmer size, how do I do that?

End trimming is unlikely to help, because that's going to shorten the reads and make them less likely to map to the target sequence. I want a way to include these reads as valid hits to the target sequence.

Answer 1

So, I think there are a few potential options here. Alevin is using selective-alignment internally to determine the mappings for the reads. So, even if you have a k-mer supporting the mapping, if you have gibberish for the rest of the read, the mapping score is going to be very poor and you're not going to recover that mapping locus. There are some parameters you can mess with to try to increase the sensitivity of mapping, though you should always be aware that upping the sensitivity a lot will mean there is less aggressive filtering early in the mapping process which can slow things down. The main parameter you should tweak to allow lower alignment scores is --minScoreFraction; the smaller this value the lower (worse) will be the permitted scores for valid mappings. Additionally, you can combine changes to the --minScoreFraction with the --softclip option to simply score the non-matching end of read as 0 (rather than a long deletion or large set of mismatches); this will lower the effect of bad read ends.

Another alternative you might consider is to try out alevin-fry for processing this data, and mapping the reads with the --sketch flag. This will only check certain structural constraints of the mapping and won't force a particular minimum alignment score. The caveat here is that dropping the requirement for alignment scoring can increase the rate of spurious mappings (reads deriving from elsewhere that spuriously map to some indexed target simply by virtue of sharing some sequence). In transcriptome wide analyses, this can be greatly mitigated by building an appropriate index and using a UMI resolution strategy that complements this expanded index. We've not examined this in detail for targeted single-cell sequencing. However, if you want to give it a try and run into any issues, feel free to drop an issue over on GitHub and we'd be happy to help work through it.

Answer 2

@ATPoint commented that the default kmer size in Salmon indexing is 31, which led me to a solution that involved altering the index generation step:

Instead of:

salmon index -t merged_salmon_${targetName}.fasta -i salmonIndex_${targetName}

I now have:

salmon index -t merged_salmon_${targetName}.fasta -i salmonIndex_k11_${targetName} -k 11

In other words, I changed the kmer size from its default of 31 to 11, which matches the BD-recommended seed substring length for the standard Bowtie2 approach ("-L 10", but I increased that to 11 because Salmon needs odd kmer sizes). I'm pretty sure that my 28bp-mapping reads won't be picked up with a kmer size of 31. Because this is targeted sequencing, I don't expect that a reduced kmer size will be an issue, but I'll keep a watch out for that.

I'm currently re-running the analysis. If it still doesn't work, then I'll see about adding in the tweaks suggested by @nomad.

Update: still not working, so I'll shift onto score tweaks.