Bowtie2 + Multiple Sequence References

Question

[elfeu]

Hello guys,

I would like to make a database with multiple contigs to map my reads against this database. However, I don't know how to use bowtie2 and samtools to select reads that mapped uniquely to a specific ORFs and reads which charted in two ORFs. How do I do this filter in samtools?

Thank you.

[Pavel Senin]

I think that you can rely on SAM flag that indicate "primary" (the best) and "secondary" alignments, also, depending on your workflow, you can configure bowtie to report only best alignments, (-M and -k options). edit: or maybe, by knowing ORFs locations you just want to segregate among mapped to them "unique" and other reads?

Answer 1

[Ashutosh Pandey]

Depending on your requirement, I would suggest you to use BWA than bowtie2. Bowtie2 has -k option that fetches the top k mapping positions. But these mapping positions may have different alignment scores. But BWA provides "X0:Number of best hits" tag in the SAM/BAM file that tells you how many positions a read can be aligned with the best alignment score. I think this is more suitable for your problem. For example, if you use -k2 option in Bowtie2, it may give you two positions but the first position may have better alignment score than the second. BWA will have X0:1 for such read as there is only one best alignment position. Also, -k2 will only fetch the first two possible alignments for that read. There can be more alignments possible. In that case, BWA 'XT' tag will give all number of suboptimal alignments.

In short, both these aligners can be used for your problem but the output from BWA will be easy to play with compared to Bowtie2.

[Devon Ryan]

I'll just add that BWA also has more coherent MAPQ scores. With bowtie2, you can actually have lower probability alignments with higher MAPQ scores, which I don't think is the case for BWA.

Answer 2

[Devon Ryan]

What definition of "mapped uniquely" do you want to go by? If by that you mean, "mapped at least slightly better to one specific location than any other", then just use reads with MAPQ scores >= 2 (bowtie2 will asign a MAPQ of 1 or 0 if a read can align equally well to multiple places). If you want a bigger difference between a best and second best alignment to call something unique, then just up that threshold (you might want a higher threshold anyway, depending on what you intend to do with the results). To extract reads above a certain threshold (say, 2) from foo.bam and put them in foo.filtered.bam:

samtools view -q 2 -b foo.bam > foo.filtered.bam

[Ashutosh Pandey]

I would agree with this answer. The solution will depend on the definition of "mapped uniquely". There can be two options: 1) two alignment positions that have the same alignment scores or 2) two alignment positions that have slightly different scores. If it is second, then use the above solution. For the first condition, BWA should help.