I am working with paired end NGS miseq data from a viral genome with multiple timepoints. I have filtered and trimmed this data for quality and adapter sequences. I have then merged the filtered paired end reads and mapped them to a reference genome using bowtie2. I then used samtools to filter for mapping quality and visualized the resulting filtered bam file on IGV.

My goal is to now scan along the genome using various window sizes (e.g. 300 nt) across the reads that have mapped to the reference. I will then chop up the genome and build phylogenetic trees based on each window.

How can I extract a multiple sequence alignment of these 300 nt windows from a sam or bam file using the coordinates of a reference viral genome? I have tried using samtools view but this does not actually provide the subsection or alignment within the window but rather all reads that overlap the window in some way.

The virus is HIV1 and the immediate goal is to get an alignment of all reads that completely span a user specified coordinate range eg. 1-300 based on some reference sequence. Ideally reads that fit this criteria would be extracted and aligned. After ragged edges are removed, the resulting alignments subjected to phylogenetic analyses.

  • $\begingroup$ this seems a little "chaotic". I would consider doing variant calling and then trying some approach? $\endgroup$
    – Colin D
    Jan 31 at 4:11

1 Answer 1


If need to get msa between every read and ref in the specific region , you can try this https://github.com/orangeSi/bam2msa.git

bam2sam ref.fa input.bam --regions chr1:200-300

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.