Finding transcriptome matches for thousands of 21-mers

Question

I am on a project to find sequence matches in mRNA databases across 20+ taxa. I was presented an EPA memorandum that outlines a method I was requested to replicate. I am not sure it is the best method and the memorandum didn't go into much detail. I played around a bit and could use a bit of advice.

Setup: I have an unspecified number (I haven't been told yet) of dsRNA segments approximately 300bp in length. I need to match these to the human transcriptome as well as over 20 other taxa.

Criteria: For each 300bp dsRNA, I am to find mRNA in the taxa that have 14 or more matches within a 21bp window. Then I sort the data by taxon, transcripts matched, and the annotation for the matched mRNA.

Approach (this is where I have questions): The EPA memorandum says the Burrows-Wheeler Aligner (BWA) was used to align a 21-mer sliding window along the target transcriptomes to look for matches of 14 or greater within the window. The PI said to create all 21-mers using a sliding window along the dsRNA sequence. Easy enough.

Here are my questions:

1. Is BWA the best approach to use? I've never used BWA MEM for anything so small.
2. How should I set the parameters for the BWA for this case? The defaults are inadequate, but I'm just taking stabs in the dark to see what falls out. So far, I have adjusted:

• Minimum seed length (-k) down to 3
• band width (-w) down to 7
• ignore alignment scores lower than (-T) range from 1 to 21
• gap open penalty (-O) between 1 and 6
• mismatch penalty (-B) between 1 and 4

3. Why do I see a Bitwise Flag of 0? In adjusting the parameters, the resulting SAM will contain matches where the Bitwise Flag is 0. This seems like nonsense to me, suggesting that I may be on the wrong track.

Sample Execution:

./bwa mem -k 5 -B 1 -O2 -T 5 ../ncbi_dataset/GCF000001405.40.rna.fna ../seqA.fasta | gzip -3 > ../bwa_results/aln_seqA.sam.gz

Bitwise Flag == 0?

seqA_331_352   16      NM_014873.3     6590    0       6S15M   *       0       0       GATCGGTGTAAATCCCATATC   *       NM:i:1  MD:Z:14A0       AS:i:14 XS:i:14
seqA_332_353   16      XM_011510229.4  7276    0       7S14M   *       0       0       TGATCGGTGTAAATCCCATAT   *       NM:i:0  MD:Z:14 AS:i:14 XS:i:14
seqA_333_354   0       XM_017008212.3  1680    0       7S14M   *       0       0       TATGGGATTTACACCGATCAA   *       NM:i:0  MD:Z:14 AS:i:14 XS:i:13
seqA_334_355   0       XM_017008212.3  1680    0       6S15M   *       0       0       ATGGGATTTACACCGATCAAC   *       NM:i:1  MD:Z:14A0       AS:i:14 XS:i:13
seqA_335_356   0       XM_017008212.3  1680    0       5S14M2S *       0       0       TGGGATTTACACCGATCAACT   *       NM:i:0  MD:Z:14 AS:i:14 XS:i:13```

Answer 1

Depending on the input 300bp query sequences (A, B, C and so on) you may get identical 21-kmers derived from the same (say A RNA) or several (A, B) sequences. To skip the repetitive searches one can cluster at 100% identity 21-kmers keeping track of the origins using i.e. vsearch If you must use bwa mem try faster bwamem2. But before plunging into using bwa or blast I suggest doing some test runs with vsearch or LAST to compare the results/do the benchmarks.