BWA - read identity

[Antony03]

Hi,

I'm using BWA to map my short sequencing reads (Illumina MiSeq) on a reference sequence. Is this possible to set the minimum % of identity with BWA? For exemple, if I set 90%, then all reads under this threshold will be discarded. I know the parameter "-n" can do something like this but the manual is not really clear (for me). Is "-n" the minimal absolute mismatch in term of base pairs?

Thank you!

Antony

[Brian Bushnell]

Since nobody else has answered this, I will mention that I do not know of a way to generically do that in BWA. The manual says:

-n NUM Maximum edit distance if the value is INT, or the fraction of missing alignments given 2% uniform base error rate if FLOAT. In the latter case, the maximum edit distance is automatically chosen for different read lengths. [0.04]

...which is clear in the case of int but unclear in the case of float. If all your reads are the same length you can calculate the edit distance that would correspond to that identity, though: e.g. 150bp * (1-0.90) = 15.

Alternately, you can use BBMap, which has an "idfilter" flag - e.g. "idfilter=0.90", which will discard reads that don't map with at least 90% identity, regardless of their length. It also has an "idtag" flag, which will annotate all the alignments in the sam file with the percent identity, so you can manually filter them if you want.

[dariober]

[This is a hack. I'm putting it here more for curiosity than anything else]

On the aligned bam file you could capture the read length and the NM tag ("Edit distance to the reference, including ambiguous bases but excluding clipping") and print the alignment if L/NM < threshold. With samtools/awk:

samtools view -h aln.bam \
| awk -v maxpct=0.05 -v OFS='\t' -v FS='\t' '{
    if($0 ~ /^@/){
        print $0
    } else {
        xnm=gensub(/.*\tNM:i:/, "", "g", $0); 
        nm=gensub(/\t.*/, "", "g", xnm); 
        pct=nm/length($10);
        if(pct < maxpct) {
            print $0, "XP:f:"pct
        }
    }
}' | less

(The L/NM goes in tag XP)