How do I get the 3bp sequence on either side of a co-ordinate

Question

[nrav58]

Hi all,

After I generate BAM or SAM files from my data, and have the coordinates of the 5' end of each read, how can I get the 3bp sequence on either side of this co-ordinate (total 6bp sequence) from the reference genome?

Answer 1

[Adam]

I recommend bedtools. Start by expanding your 5' coordinates by the 3bp.

bedtools slop -b 3 -i five_prime.bed -g my_genome.fa.fai > five_prime_window.bed

Then, get those sequences from your genome.

bedtools getfasta -fi my_genome.fa -bed five_prime_window.bed