Hello,
I am currently using the bwasw aligner to map long 454 rna-seq reads on the human genome. Due to the 400-600bp length of our reads, many on them span several exon/exon junctions, which causes some problems with processed pseudogenes and paralogs. In addition, 454 sequencers tend to produce erroneous indels around homopolymers. I am looking for a reasonable setting for the scoring parameters a, b, q and r.
Using the default values, anytime a given sequence can be aligned on both a spliced gene (several small high scoring segment pairs across different exons with no or few mismatches) and an unspliced read (a single, long high scoring segment pair with many mismatches), bwasw always chooses the unspliced alignment, even if it implies an important number of mismatches. Since bwasw only reports the best alignment, it is not possible to search for alternative alignments in the resulting BAM file. By contrast, in this particular case, BLAT almost always finds the correct alignment, because it handles introns, and reports a single alignment with the highest score, consisting of several distant blocks. However, it is untractable to align complete transcriptomes with BLAT.
Masking pseudogenes is not a completely satisfactory option because their detection rests on somewhat arbitrary thresholds and relevant databases are fairly dissimilar.
I tried to increase the mismatch penalty (with the -b flag) of bwasw, but the alignment is still incorrect, and raising this parameter too much makes bwasw replace mismatches with pairs of indels because, in this case, the penalty induced by an insertion and a deletion is still inferior to the mismatch penalty. Here is an example of mismatch turning into an indel:
reference: CCCGATT
query: CCCGGTT
alignment with default parameters (-b 3) --> a A/G mismatch appears in the alignment:
CCCGATT reference
CCCGGTT query
alignment with increased mismatch penalty (-b 15) --> the A/G mismatch is removed, but both an insertion and a deletion are being added:
CCC-GATT reference
CCCGG-TT query
There is a discussion in another thread between aleferna and lh3 about BLAT sensitivity with custom -z and -fastmap settings, but is it relevant to use these settings for 454 rna-seq reads ?
Any advice/info, about:
-specific settings,
-a good strategy to find one,
-relevant information on bwasw scoring
-ways to tweak BWA to get all alignments, not only the best one
-alternative alignment algorithms
is welcome.
Best regards
David
I am currently using the bwasw aligner to map long 454 rna-seq reads on the human genome. Due to the 400-600bp length of our reads, many on them span several exon/exon junctions, which causes some problems with processed pseudogenes and paralogs. In addition, 454 sequencers tend to produce erroneous indels around homopolymers. I am looking for a reasonable setting for the scoring parameters a, b, q and r.
Using the default values, anytime a given sequence can be aligned on both a spliced gene (several small high scoring segment pairs across different exons with no or few mismatches) and an unspliced read (a single, long high scoring segment pair with many mismatches), bwasw always chooses the unspliced alignment, even if it implies an important number of mismatches. Since bwasw only reports the best alignment, it is not possible to search for alternative alignments in the resulting BAM file. By contrast, in this particular case, BLAT almost always finds the correct alignment, because it handles introns, and reports a single alignment with the highest score, consisting of several distant blocks. However, it is untractable to align complete transcriptomes with BLAT.
Masking pseudogenes is not a completely satisfactory option because their detection rests on somewhat arbitrary thresholds and relevant databases are fairly dissimilar.
I tried to increase the mismatch penalty (with the -b flag) of bwasw, but the alignment is still incorrect, and raising this parameter too much makes bwasw replace mismatches with pairs of indels because, in this case, the penalty induced by an insertion and a deletion is still inferior to the mismatch penalty. Here is an example of mismatch turning into an indel:
reference: CCCGATT
query: CCCGGTT
alignment with default parameters (-b 3) --> a A/G mismatch appears in the alignment:
CCCGATT reference
CCCGGTT query
alignment with increased mismatch penalty (-b 15) --> the A/G mismatch is removed, but both an insertion and a deletion are being added:
CCC-GATT reference
CCCGG-TT query
There is a discussion in another thread between aleferna and lh3 about BLAT sensitivity with custom -z and -fastmap settings, but is it relevant to use these settings for 454 rna-seq reads ?
Any advice/info, about:
-specific settings,
-a good strategy to find one,
-relevant information on bwasw scoring
-ways to tweak BWA to get all alignments, not only the best one
-alternative alignment algorithms
is welcome.
Best regards
David