<p>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. <a href="https://github.com/torognes/vsearch" rel="nofollow noreferrer">vsearch</a>
If you must use <code>bwa mem</code> try faster <code>bwamem2</code>.
But before plunging into using <code>bwa</code> or <code>blast</code> I suggest doing some test runs with <code>vsearch</code> or <a href="https://gitlab.com/mcfrith/last" rel="nofollow noreferrer">LAST</a> to compare the results/do the benchmarks.</p>

CC BY-SA 4.0
2024-09-25T17:03:57.380Z
2024-09-25T17:03:57.380Z
darked89
<p>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.</p>
<p><strong>Setup:</strong>
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.</p>
<p><strong>Criteria:</strong>
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.</p>
<p><strong>Approach (this is where I have questions):</strong>
The EPA memorandum says the <em>Burrows-Wheeler Aligner</em> (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.</p>
<p><strong>Here are my questions:</strong></p>
<ol>
<li><em>Is BWA the best approach to use?</em> I've never used BWA MEM for anything so small.</li>
<li><em>How should I set the parameters for the BWA for this case?</em> The defaults are inadequate, but I'm just taking stabs in the dark to see what falls out. So far, I have adjusted:</li>
</ol>
<ul>
<li><em>Minimum seed length (-k)</em> down to 3</li>
<li><em>band width (-w)</em> down to 7</li>
<li><em>ignore alignment scores lower than (-T)</em> range from 1 to 21</li>
<li><em>gap open penalty (-O)</em> between 1 and 6</li>
<li><em>mismatch penalty (-B)</em> between 1 and 4</li>
</ul>
<ol start="3">
<li><em>Why do I see a Bitwise Flag of 0?</em> In adjusting the parameters, the resulting SAM will contain matches where the <em>Bitwise Flag</em> is 0. This seems like nonsense to me, suggesting that I may be on the wrong track.</li>
</ol>
<p><strong>Sample Execution:</strong></p>
<p><code>./bwa mem -k 5 -B 1 -O2 -T 5 ../ncbi_dataset/GCF000001405.40.rna.fna ../seqA.fasta | gzip -3 &gt; ../bwa_results/aln_seqA.sam.gz</code></p>
<p><strong>Bitwise Flag == 0?</strong></p>
<pre><code>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```
</code></pre>

CC BY-SA 4.0
2024-09-25T12:27:03.560Z
What is it you consider small for bwa? If anything, 300bp is on the long side. And you only have single, individual sequences, right? Not multiple reads with the same sequence? What am I missing? And have you considered using blast instead? It's the kind of thing blast is made for.
terdon
CC BY-SA 4.0
2024-09-26T16:46:57.757Z
The 300bp sequence isn't fed into BWA. Rather, 21-mers are generated using a sliding window. It is the 21-mers that are put into BWA.
Greymouse451
<transcriptome><bwa><k-mer>
Finding transcriptome matches for thousands of 21-mers
# Finding transcriptome matches for thousands of 21-mers

(Asked by: username_0 on 2024-09-24)

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:

Is BWA the best approach to use? I've never used BWA MEM for anything so small.
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


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```

**Tags:** transcriptome, bwa, k-mer
> **Comments:**
> **username_1**: What is it you consider small for bwa? If anything, 300bp is on the long side. And you only have single, individual sequences, right? Not multiple reads with the same sequence? What am I missing? And have you considered using blast instead? It's the kind of thing blast is made for.  
> **username_0**: The 300bp sequence isn't fed into BWA. Rather, 21-mers are generated using a sliding window. It is the 21-mers that are put into BWA.  

---

## 1 answers

### Answer by username_2 on 2024-09-25. Score: 0

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](https://github.com/torognes/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](https://gitlab.com/mcfrith/last) to compare the results/do the benchmarks.


---

