align kmers to reference genome

Question

[glowieelife]

dear all,

I have a fasta file with a list of 25-mers and I am trying to align it to the reference genome ref.fa using bowtie2

I did bowtie2 -x ref -f ref_25mers.fa -S ref_25mers.sam but it gives the result

1 reads; of these:
  1 (100.00%) were unpaired; of these:
    1 (100.00%) aligned 0 times
    0 (0.00%) aligned exactly 1 time
    0 (0.00%) aligned >1 times
0.00% overall alignment rate

which is not supposed to happen as the 25mers come from the reference genome itself... does anyone know what is the correct way to align these 25mers to the reference genome?

thank you so much

ref_25mers.fa:

>ref_25mers
AGCCTCAGGAAGGAGGCAGTGCTGC
GCCTCAGGAAGGAGGCAGTGCTGCC
CCTCAGGAAGGAGGCAGTGCTGCCA
CTCAGGAAGGAGGCAGTGCTGCCAG
TCAGGAAGGAGGCAGTGCTGCCAGC
CAGGAAGGAGGCAGTGCTGCCAGCC
AGGAAGGAGGCAGTGCTGCCAGCCC
GGAAGGAGGCAGTGCTGCCAGCCCT
GAAGGAGGCAGTGCTGCCAGCCCTT
AAGGAGGCAGTGCTGCCAGCCCTTG
AGGAGGCAGTGCTGCCAGCCCTTGG
GGAGGCAGTGCTGCCAGCCCTTGGG
GAGGCAGTGCTGCCAGCCCTTGGGG
AGGCAGTGCTGCCAGCCCTTGGGGA
GGCAGTGCTGCCAGCCCTTGGGGAC
GCAGTGCTGCCAGCCCTTGGGGACA
CAGTGCTGCCAGCCCTTGGGGACAA
AGTGCTGCCAGCCCTTGGGGACAAC
GTGCTGCCAGCCCTTGGGGACAACA
TGCTGCCAGCCCTTGGGGACAACAG
GCTGCCAGCCCTTGGGGACAACAGC
CTGCCAGCCCTTGGGGACAACAGCC
TGCCAGCCCTTGGGGACAACAGCCT
GCCAGCCCTTGGGGACAACAGCCTG
CCAGCCCTTGGGGACAACAGCCTGT
CAGCCCTTGGGGACAACAGCCTGTC
AGCCCTTGGGGACAACAGCCTGTCC
GCCCTTGGGGACAACAGCCTGTCCC

[Pierre Lindenbaum]

>ref_25mers
AGCCTCAGGAAGGAGGCAGTGCTGC
GCCTCAGGAAGGAGGCAGTGCTGCC
CCTCAGGAAGGAGGCAGTGCTGCCA
CTCAGGAAGGAGGCAGTGCTGCCAG
TCAGGAAGGAGGCAGTGCTGCCAGC
CAGGAAGGAGGCAGTGCTGCCAGCC

this is not a multifasta file of N-kmers of 25 bp, this is just ONE big fasta file with lines of 25 characters.

[Pierre Lindenbaum]

if you just want to have the locations, with perfect match, you could could "just" map each of those kmer using a boyer-moore algorithm.

[GenoMax]

You probably want un-gapped alignments so instead of bowtie v.2.x use v.1.x.

Answer 1

[liorglic]

Not sure if that's the best/fastest way, but I use BWA like this:

# map  
bwa aln <genome.fasta> <kmers.fasta> > <result.sai>  
bwa samse <genome.fasta> <result.sai> <kmers.fasta> -f <mapping.sam>  
# filter for perfect mapping  
samtools view -h -e 'flag.unmap != 4 && qlen == 25 && [NM] && [NM] == 0' <mapping.sam> > <mapping.filter.sam>

As noted above, your example input is currently not in valid fasta format - each kmer needs to have its own header line, so before you map you must assign IDs add add the fasta headers.