The data

I have a large number of reads from sequences that were generated by randomly shuffling regions of two parent sequences together. See the following image:


The regions that are shuffled are around 100-200 bp in length. The full shuffled sequence is around 1000 bp long. In general, the regions will be in the same order as the original sequence. The shuffled sequences are fragmented and sequenced into 150bp reads. One read can contain more than two borders between shuffle regions from a different parent sequence.

The sequences are very similar (95%). Because of this high similarity you cannot pinpoint the border between two regions. The border will be a range between a variant that is unique to source1 and a variant that is unique to source2. That is ok because I am not interested in knowing from which source an identical region came.

The question

I would like to identify which part of the read originates from parent sequence one or two. What would be a suitable tool to do this?

My attempts

I tried to look at some multiple sequence alignment tools like MAFFT or Clustal but I am not interested in aligning all reads to each other. They also don't identify the shuffle regions well in my testing.


1 Answer 1


This is the sort of thing that LAST has been specifically designed to work well for:

LAST is designed for moderately large data (e.g. genomes, DNA reads, proteomes). It's especially good at finding rearrangements and recombinations: we believe last-split does that more rigorously than anything else.

Here's my usual approach for mapping with LAST:

lastdb -uRY4 -R01 reference_sequences.fa reference_sequences.fa
# training is probably not needed if your sub-sequences have no
# variation from the parental sequences
last-train reference_sequences.fa input_sequences.fa > trained.mat
lastal -P 10 -p trained.mat --split reference_sequences.fa input_sequences.fa > aligned.maf

LAST is very customisable. If you want to adjust the mapping or split sensitivity, have a look at the lastal documentation here, and the last-split documentation here.

For example, there's a mismap probability adjustment which excludes alignments with a mismap probability above a certain threshold:

# only form alignments if the mismap probability is 0.1% or less
lastal --split-m 0.001 ...
  • 1
    $\begingroup$ I tried a lot of the parameters. The mismap threshold did not have an impact for me. I found that lowering the minimum alignment score with the parameter lastal -e to 50 worked well to increase the sensitivity. It is now working perfectly. $\endgroup$
    – Yano
    Nov 3, 2022 at 16:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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