SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Lots of unmapped reads - SOLiD bacterial RNA-seq and bowtie mapping Jean RNA Sequencing 10 01-17-2013 11:11 AM
BWA - unmapped reads Adamo Bioinformatics 23 02-22-2012 09:12 PM
BWA output bitwise flag for mapped/unmapped reads wenhuang Bioinformatics 1 08-29-2011 03:54 PM
BWA behaviour with Mate Pair data + Multi read mapping apratap Bioinformatics 5 06-23-2011 04:38 PM
BWA, mostly unmapped reads Michael.James.Clark Bioinformatics 4 03-03-2011 02:20 PM

Reply
 
Thread Tools
Old 01-13-2012, 07:40 AM   #1
kjlee
Member
 
Location: Atlanta, GA

Join Date: Jun 2011
Posts: 12
Default BWA rescue of multi-mapping or unmapped reads

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

Quote:
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
kjlee is offline   Reply With Quote
Old 01-14-2012, 06:02 AM   #2
nilshomer
Nils Homer
 
nilshomer's Avatar
 
Location: Boston, MA, USA

Join Date: Nov 2008
Posts: 1,285
Default

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).
nilshomer is offline   Reply With Quote
Old 01-16-2012, 10:15 AM   #3
kjlee
Member
 
Location: Atlanta, GA

Join Date: Jun 2011
Posts: 12
Default

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.

Quote:
#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 is offline   Reply With Quote
Old 01-16-2012, 10:16 AM   #4
kjlee
Member
 
Location: Atlanta, GA

Join Date: Jun 2011
Posts: 12
Default

Thanks for the advice, Nils.
kjlee is offline   Reply With Quote
Old 01-16-2012, 08:26 PM   #5
nilshomer
Nils Homer
 
nilshomer's Avatar
 
Location: Boston, MA, USA

Join Date: Nov 2008
Posts: 1,285
Default

It would be nice to have them as command line options to see what effect there is when varying these values.
nilshomer is offline   Reply With Quote
Old 04-26-2012, 04:33 PM   #6
vyellapa
Member
 
Location: phoenix

Join Date: Oct 2011
Posts: 59
Default

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.

Code:
Samtools spec says
Added `mpileup -A' to allow to use reads in anomalous pairs in SNP calling.
vyellapa is offline   Reply With Quote
Old 07-13-2012, 10:29 PM   #7
Michael.James.Clark
Senior Member
 
Location: Palo Alto

Join Date: Apr 2009
Posts: 213
Default

Quote:
Originally Posted by nilshomer View Post
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!
__________________
Mendelian Disorder: A blogshare of random useful information for general public consumption. [Blog]
Breakway: A Program to Identify Structural Variations in Genomic Data [Website] [Forum Post]
Projects: U87MG whole genome sequence [Website] [Paper]
Michael.James.Clark is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 09:47 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO