3
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ I cannot say that I follow why these particular reads are actually 28bp followed by gibberish but the default kmer size in salmon indexing is 31, hence reads shorter than inevitably go unmapped. You could index with a smaller k. Why exactly do these reads have these polyGs? Is that a consequence of how these targeted libraries are made, so is the fragment only 28bp of meaningful content? $\endgroup$
    – ATpoint
    Jun 7 at 20:08
  • $\begingroup$ Since the end is gobbledegook, you should quality trim. $\endgroup$
    – swbarnes2
    Jun 7 at 22:53
4
$\begingroup$

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.

$\endgroup$
2
  • 1
    $\begingroup$ Thank you; tweaking the minScoreFraction did the trick. I went from a max Cxcr4 read count per cell of 1 with default parameters, to 4 with minScoreFraction=0.5, to 26 with minScoreFraction=0.33. $\endgroup$
    – gringer
    Jun 8 at 8:47
  • $\begingroup$ Out of curiosity, did that require changing the kmer length during indexing? $\endgroup$
    – ATpoint
    Jun 10 at 8:50
0
$\begingroup$

@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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.