BAM format file to FASTA alignment file

Question

[ben.ward]

Hi,

I've been given a reference genome and BAM file(s) which are reads aligned to that reference. This will be the first time I've done anything with BAM/SAM formats. I'd like to know how I would go about generating for each BAM file, the full sequence, using the reads in the BAM file and the reference they are aligned to, essentially resulting in a FASTA format alignment of the reference, and then each sequence represented by each BAM file. What are the considerations? For example at the top of my head, whether reads are paired end and their strandedness matters.

I'd like to know if this can be done in R / Bioconductor as I have some familiarity with it - mostly using Biostrings.

Thanks,

Ben.

[Lee Katz]

I re-read your question and I think you actually want a fasta format alignment. You actually probably don't want that. You are new to BAM/SAM and so you might not realize yet that it would be pretty awful having it in a different alignment format. Take a look at the actual file or the first several lines and just see what the file means before moving onto something more familiar like fasta.

e.g., samtools view file.bam | head -n 999 | column -t | less -S

Answer 1

[Devon Ryan]

The simplest method for generating the consensus sequence is to simply use bcftools consensus. An example is given in the documentation. Obviously, you'll need to have called the variants. GATK has a similar tool.

[Lee Katz]

In samtools v1.1, you can use samtools bam2fq.

Usage:   samtools bam2fq [-nO] [-s <outSE.fq>] <in.bam>

Options: -n        don't append /1 and /2 to the read name
         -O        output quality in the OQ tag if present
         -s FILE   write singleton reads to FILE [assume single-end]

[Devon Ryan]

I've interpreted OP as requesting the consensus sequence rather than reproducing fastq files.

[al-ash]

Bcftools consensus command does:

Create consensus sequence by applying VCF variants to a reference fasta file

so I think that it is far from being straightforward to use this tool for obtaining the consensus of mapped reads which does not contain bases from the reference sequence - which I believe is what OP was looking for. I'm now solving the same task (to reconstruct consensus of mapped reads and not reference sequence modified by VCF variants). I'm using samtools mpileup + bcftools:

samtools mpileup -ABuf reference.nt mapped.bam | bcftools call -cOz --pval-threshold 0.99 > mapped.vcf.gz
tabix mapped.vcf.gz 
cat reference.nt | bcftools consensus mapped.vcf.gz > mapped.fastq

... and the only consensus I can get still contains nucleotides from reference at some positions. In case you could suggest a specific commands or parameters how to use bcftools consensus to generate read consensus that would be really helpful.