BWA rescue of multi-mapping or unmapped reads

[kjlee]

In reading the BWA manual, I have come across a question that I still cannot answer. It reads:

In paired-end alignment, BWA pairs all hits it found. It further performs Smith-Waterman alignment for unmapped reads with mates mapped to rescue mapped mates, and for high-quality anomalous pairs to fix potential alignment errors.

So BWA will rescue mates presumably if there is a high-quality, uniquely mapping mate pair. My question is: what are the parameters that the unmapped mate have to meet to be considered correctly placed? For instance, if I have 2x100 reads, and one maps uniquely, does the other one have to map with the same stringency as I specified in the aln stage. Or if say, only 30nt of the unmapped read map, and the other 70nt do not (due to say, an inversion), will the unmapped read be mapped and the 70 non-mapping nucleotides just be soft-clipped in the BAM file?

Any help would be greatly appreciated.

Cheers,
Kevin

[nilshomer]

In some cases, the source code is the best documentation. I assume you mean in "bwa sampe", so check out the "bwa_paired_sw" function. It looks like unmapped mates are rescued if the mapping quality of the other end is greater than SW_MIN_MAPQ (set to 17).

[kjlee]

I believe that I found what I was looking for. Per Nils' suggestion, I went to the source code.

My primary query concerned the length of the match (of an unmapped read that was re-mapped using a Smith-Waterman alignment in the vicinity of a uniquely mapped read) that was required for BWA to place it.

#define SW_MIN_MATCH_LEN 20
#define SW_MIN_MAPQ 17
...
bwa_cigar_t *bwa_sw_core(bwtint_t l_pac, const ubyte_t *pacseq, int len, const ubyte_t *seq, int64_t *beg, int reglen, int *n_cigar, uint32_t *_cnt)
{
...
// check whether there are too many N's
if (reglen < SW_MIN_MATCH_LEN || (int64_t)l_pac - *beg < len) return 0;
for (k = 0, x = 0; k < len; ++k)
if (seq[k] >= 4) ++x;
if ((float)x/len >= 0.25 || len - x < SW_MIN_MATCH_LEN) return 0;
...
if (x < SW_MIN_MATCH_LEN || y < SW_MIN_MATCH_LEN) { // not good enough
free(path); free(cigar); free(ref_seq);
*n_cigar = 0;
return 0;

So basically, if there are more then 20 bases that match (in a Smith-Waterman re-alignment) the unmapped (or improperly mapped) mate will be rescued.

Cheers.

[kjlee]

Thanks for the advice, Nils.

[nilshomer]

It would be nice to have them as command line options to see what effect there is when varying these values.

[vyellapa]

What do anomalous pairs mean. Samtools mpileup is removing such reads(when compared to seeing these in IGV) unless -A option is added but Im not sure what it means.

Samtools spec says
Added `mpileup -A' to allow to use reads in anomalous pairs in SNP calling.

[Michael.James.Clark]

It would be nice to have them as command line options to see what effect there is when varying these values.

Indeed it would be!