Aligning to specific chromosomal regions using BWA.

[oliver_ngs]

Hi,

I have done some Illumina next generation sequencing after SureSelect target enrichment, and I want to align my reads to a specific region of the genome using BWA.

My fasta file consist of the specific 3 Mb targeted. If the name in the fasta file is, example >chr7:15000000-18000000, BWA doesn't seem to recognize this as the target interval and assumes the first base is at position one. The alignments are then in the wrong position when I try to import into IGV, and the gene annotations are all incorrect.

Is there a way of getting BWA to recognize the target interval, so alignment will be given the correct genomic position if I just use part of a chromosome as the reference sequence? Or do I need to align my data to an entire chromosome/genome so that the positions are correct?

Any suggestion would be very much appreciated!

Cheers,
Oliver

[ttnguyen]

I think a simple way to do is aligning your data to entire genome then only select the alignments that belong to your specific regions.

Other way maybe building a new index for your fasta file and then map your reads to this and after mapping step, you need to convert the 'local' coordinate to the 'global' coordinate.

[gaffa]

I doubt that you could make BWA assign coordinates differently. But you could possibly change the coordinates yourself in the resulting SAM file, i.e. add 15000000-1 to every coordinate or whatever.

However, you might want to consider mapping to the whole genome instead, for conceptual reasons. Target enrichment is not 100% efficient, so you are at very high risk of mis-mapping reads if you use only a specific region as your reference.

There has been some discussion of this on Biostar:

http://biostar.stackexchange.com/questions/7572/is-there-any-reference-exome/7617#7617

http://biostar.stackexchange.com/questions/4413

[krobison]

Aligning to a subset of the genome is always a bad idea. Remember that the aligner tries to find the best match between your sequence and the given genome. With a subset, you run a risk of getting an incorrect alignment between your sequence and the restricted reference which would not occur if the aligner had the entire genome (and would therefore find the correct match). The speed improvement you will get will be modest & the downstream informatics more painful.

If you are really, really stubborn about this, then generate masked versions of your chromosomes where everything except your target regions (and a few hundred basepairs of neighboring sequence) are replaced by N. That way your coordinates are preserved.

[oliver_ngs]

Thank you all very much for the advice. I'll try aligning to the entire genome. Cheers, Oliver.