Hi, I created a reference of a genome using bowtie-build and i want to map the reads to the genome for which i used bowtie command
bowtie2 -x reference -1 left.fq -2 right.fq -S out.sam
i am interested in the unmapped reads from it.I am looking for a very strict parameter like having 100% alignment match of reads to genome for it to be considered valid mapping. could anyone please suggest which parameters to be used for the same. Looking forward, Thnx
2 answers
In my opinion, in this case, the selection of alignments should be done by filtering existing alignments and not by "forcing" the aligner to behave a certain way. Look into the SAM format specification and note how imperfect alignments are annotated.
For example, the NM tag will contain the value of the so-called edit distance to the reference. When that value is zero the query matches perfectly.
If, as in your definition, every alignment with a change should be called unaligned then filter for the NM tag value that is larger than 0 and select your aligments that way. If retaining a BAM file is essential you can even flip the flag to be unmapped.
Check out the penalty parameters --mp, --rdg and --rfg which you can set to something like 10000, the number of allowed seed mismatches -N which you can set to 0 and the option --very-sensitive to enforce stringent end-to-end alignment (the latter actually already contains seed mismatch of 0). Together these should discard every non-perfect alignment. Additionally you could be stringent with the post-filtering towards MAPQ to only keep very unique reads, e.g. with samtools view -q 30.
Log in to answer this question.
It may also help to add some context for why you want to do this? There may be a better option instead of doing this.