Hi all,
In one (or actually two) of our BAM files we came across something unexpected. We mapped a pair of FASTQ's to two different small references. The references had highly homologous regions in them, and this is about a read mapping to one of those homologous bits.
We found that some reads that mapped to one of the highly homologous regions did show up in the BAM for one reference, but not the other. The aligned segment mapped identically well to both references after blasting.
While inspecting the aligned section in the BAM where it did show up, I noticed that the alignment score (AS:i) was 30, while I had mapped with -T 33 (bwa mem), this struck me as odd; the alignment should never have passed filtering.
When then trying to reproduce and isolating the identified read, mapping it to both references, the read suddenly behaved identical for both references. Creating an alignment when I set -T to 29, and not creating an alignment when -T was set to 33, identical for both references.
When then thinking it was maybe a hiccup or bit-flip and just rerunning the original job (complete fastq to both references), the read behaved differently again; creating an alignment in one reference with AS 30, but not in the other (in both cases I'm sure I put -T to 33).
The flags for the BAM where it the alignment showed: 163/83 for the mates of the pair The flags for the BAM where the alignment did not show: 181/121 for the mates of the pair (the mate got mapped in all cases with an AS of 37)
What could have caused bwa to make different filtering decisions between the mappings? Why did it allow an alignment with AS30 while I had set -T 33?
Any insights or ideas on what to look into further get to the core might be appreciated!
Curious if the data from a larger reference i.e. are you aligning to a reduced representation of the original genome.
Yes, we are mapping to a subset of the original genome where the reads came from.
The aligned segment in question does not align to another part of the complete reference.