How do I find the most similar sequences from a large set of short sequences?

Question

I've been asked to evaluate whether the MinION will be sufficient to distinguish between 20bp CRISPR guide RNAs from the GeCKO v2 set. I know that this set has some exactly identical sequences, which wouldn't be distinguishable regardless of the sequencing method used, so I would like to know whether any two non-indentical sequences differ by less than 15% of the 20bp (i.e. only 3 bases different).

There are around 130k 20bp sequences in the dataset, and I'm pretty sure that an all-vs-all approach is not feasible, but perhaps there's some trick that I'm missing.

What I thought of doing was looking at dimer and trimer/codon counts to hunt down sequences that were most similar. I found about 100 sequences that shared the same trimer count signature (but were not identical), and within this group found [only] one pair with 85% identity; all the remainder were 70% or less.

... but I'm having trouble ignoring the doubt I have that these two sequences are the most similar. It might be possible that there are sequences in the dataset that don't have identical trimer count signatures, but differ by less than three bases. I'll give a very obvious example (which isn't in the dataset):

1: AAAAA AAAAA AAAAA AAAAA
2: AAAAC AAAAA AAAAA AAAAA

The trimer counts for these two sequences are different (Sequence 1 -- AAA: 18; Sequence 2 -- AAA: 15; AAC: 1; ACA: 1; CAA: 1), and yet the sequences are only 1 base different (which might be a problem for the MinION). How can I work out if any such sequence pairs exist in the dataset?

Answer 1

130k * 20bp is a small data set. At this scale, SSE2 Smith-Waterman may work well:

git clone https://github.com/attractivechaos/klib
cd klib && gcc -O2 -D_KSW_MAIN ksw.c -o ksw -lz
./ksw -a1 -b1 -q1 -r1 -t14 20bp.fa 20bp.fa > out.tsv

I simulated 130k 20bp reads from E. coli K-12. The last command line takes about ~2 CPU hours based on the partial output.

It is possible to make this much faster by using 5bp seeds. You can also try ssearch and fasta from the fasta aligner package. They are slower, probably because they are not optimized for such input. Computing Hamming distance with SSE2 will be much faster, but is probably of little use as it does not allow gaps. There are also Gene Myers' edit-distance based algorithms, which can be faster than SSE2-SW.

Bowtie2 (as well as most short-read mappers) won't work well at the default setting. Bowtie2 uses long exact seeds. It will miss many 1-mismatch hits, let alone 3-mismatch ones. The number you get from bowtie2 is an underestimate. You might be able to tune bowtie2, but to find 3-mismatch hits, you have to make it a lot slower.

Answer 2

For this application, you could probably also do something like calculate the Hamming distances between all of the strings in an all-vs-all approach (it should not take too long or too much overhead). You could use something like the Hamming distance tools in Julia. Here is an example of what I mean (using Julia):

using StringDistances

k = ["AATTGGCC", "AATTGGCA", "AATTCCGG"]

for s in k
    for y in k
        print(s, " compared to ", y, ":  ", compare(Hamming(), s, y), "\n")
    end
end

The output looks like this:

AATTGGCC compared to AATTGGCC:  1.0
AATTGGCC compared to AATTGGCA:  0.875
AATTGGCC compared to AATTCCGG:  0.5
AATTGGCA compared to AATTGGCC:  0.875
AATTGGCA compared to AATTGGCA:  1.0
AATTGGCA compared to AATTCCGG:  0.5
AATTCCGG compared to AATTGGCC:  0.5
AATTCCGG compared to AATTGGCA:  0.5
AATTCCGG compared to AATTCCGG:  1.0

This is just a simple small example but the core is there. Should not be too bad scaling it up. Hope this helps and provides another helpful option. :)

Answer 3

I'm pretty sure that an all-vs-all approach is not feasible, but perhaps there's some trick that I'm missing.

The "trick" is that I hadn't actually tried this, and Bowtie2 happens to be quite good at doing this. When there are other similar matches, the MAPQ score that Bowtie2 produces are reduced, so all that is needed is to identify those sequences with less than a MAPQ of 42. I've found 2810 such sequences that share some similarity to others (2416 if I exclude reverse-complement matches):

$ bowtie2-build Mouse_GeCKOv2_Library.fasta Mouse_GeCKOv2_Library.fasta
$ bowtie2 -f -x Mouse_GeCKOv2_Library.fasta -U Mouse_GeCKOv2_Library.fasta | 
   samtools view | awk '{if($5 != 42){print ">"$1"\n"$10}}' > similar_seqs.fasta
125795 reads; of these:
  125795 (100.00%) were unpaired; of these:
    1 (0.00%) aligned 0 times
    122982 (97.76%) aligned exactly 1 time
    2812 (2.24%) aligned >1 times
100.00% overall alignment rate