How can I filter out reference sequences that get multiple hits in BWA MEM ?

[Chrys]

Hi everyone !

I have exported a small number (~500) of 145 pb sequences from a Stacks analyse.

I then used it as my "reference genome" (actually a very fragmented and unordered one, consequently !) to align a .fastq.gz file containing several millions of reads 1. My objective is to filter out, from my "reference", every single 145 pb sequence from which an aligned read actually has multiple hits.

e.g if a given read has a primary alignment on Seq A from reference, but also a secondary alignment on Seq B from reference, I want to filter out both Seq A and B. In the end, I hope to obtain a smaller "reference" with very specific sequences that I can amplify unambiguously by PCR.

I used BWA MEM to do the alignment, with parameters -a and -M turned on, but I would like to know if there would be another way (than making a script that works directly on the SAM file) to detect all sequences that I want to get rid of.

This document about SAM format and SAMtools manual are both very helpful, and I have read a lot of threads on forums including this one but I couldn't find any solution for my issue. Would somebody have any suggestion :-) ?

Thanks a lot !

Chrys

[HESmith]

Point of clarification: assuming that your reads are shorter than 145bp, some of your reads may uniquely match a portion of reference read 'A' while other reads multi-map to references 'A' plus 'B'. Do you want to keep or discard reference 'A'?

[Chrys]

Hello and thank you so much for answering me !

In this case, I would discard reference A. Indeed my ultimate goal is to avoid as much as possible cross amplification when performing PCR, and to avoid the risk to pick up a reference that has a portion highly similar to another.

I know this is not very subtle, but I'm still discovering sequence alignment and I want to start on an easy level

[HESmith]

The simplest approach is probably to output the multi-mapping reads using SAMtools 'view' to filter (an explanation of SAM flags can be found here), then use BEDtools 'intersect -v' to output the reference elements that are not contained in your multi-mapper BAM.