I am not talking about consensus sequence, I know how to get consensus sequence using mpileup in samtools/bcftools. As I understand , SAM/BAM files are basically sequence alignment format so it's natural to expect a straightforward way of conversion between them. Sadly I have yet to find one. Vcf-kit allows for conversion between vcf file and fasta alignment, but the fasta alignment contains only the SNPs. I need the whole alignment. Does samtools/bcftools provide option for such conversion? If not, is there any other tool?
I am not sure what you mean by "fasta alignment file". If you mean a multi-sequence alignment (MSA) in the fasta format, you can't get that because SAM keeps pairwise alignments only and doesn't align inserted sequences. Even if you don't care about inserted sequences, a MSA in fasta is far to big to be practical. Alternatively, by "fasta alignment file", you could mean pairwise alignment, but it is still impractical to output every pairwise alignment in a separate fasta file.
There are tools to convert SAM to BLAST-like format if that is what you want. The following shows one of them. I believe there are other tools easier to use, but I forget what are they.
Anyway, to install:
wget https://raw.githubusercontent.com/lh3/minimap2/master/misc/paftools.js curl -L https://github.com/attractivechaos/k8/releases/download/v0.2.4/k8-0.2.4.tar.bz2 | tar -jxf - cp k8-0.2.4/k8-`uname -s` k8
Then you can convert your
./k8 paftools.js sam2paf -L ALN.sam | ./k8 paftools.js view -l0 - | less -S
Or generate a SAM stream from
ALN.bam and pipe it with:
samtools view -h ALN.bam | ./k8 paftools.js sam2paf -L - | ./k8 paftools.js view -l0 -
Here is an example output:
>MT_orang 16499 0 16025 + MT_human 16569 576 16569 13773 16095 60 tp:A:P mm:i:2150 gn:i:172 go:i:101 Ref+: 577 GTTTATGTAGCTTACCTCCT---CAAAGCAATACACTGAAAATGTTTAGACGGGCTCACATCACCCCATAAACAAATAGGTTTGGTCCTAGCCTTTCTATTAGCTCTTAGTAAGATTACACATGCAAGCATCCCCGTTCCAGTGAGTTCACCCTCTAAATCACCACGATCAAAAGGAACAAGCATCA |||||||||||||| | || ||||||||| |||||||||||| | ||||||| |||| | |||||||||||||||||||||||||||||||||||||||||||||||| || ||||||||||||||||||||| |||||||| || ||||| || |||| || || |||| ||||||||| Qry+: 1 GTTTATGTAGCTTA--TTCTATCCAAAGCAATGCACTGAAAATGTCTCGACGGGCCCACA-CGCCCCATAAACAAATAGGTTTGGTCCTAGCCTTTCTATTAGCTCTTAGTGAGGTTACACATGCAAGCATCCCCGCCCCAGTGAG-TCGCCCTCCAAGTCACTCTGACTAAGAGGAGCAAGCATCA //
You can change
-l80 to print alignments in multiple lines.
InIGV you can visuals the alignments of individual reads to the reference as contained in the .bam file. I guess what you want to is get that information into a file containing all the (paired) read sequences aligned to a particular region, with padding for the missing data. That way you could inspect all SNPs or sequencing errors, check possible haplotypes across a few SNPs, use these for phylogenetic analysis etc.
I am also looking for such tool and best I could find so far is an answer on BioConductor. It refers to some code in R that can be used for this purpose:
MSA fasta file can be generated from SAM, but it is almost impractical to generate this HUGE file. If you still insist to generate the file, then please check this.