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.