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?
$\begingroup$
$\endgroup$
Add a comment
|
3 Answers
$\begingroup$
$\endgroup$
6
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 ALN.sam with:
./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 -l0 to -l80 to print alignments in multiple lines.
answered Feb 17, 2019 at 15:04
-
$\begingroup$ I meant "aligned fasta format" which I think is MSA. Beside general curiosity, my immediate need for the format is as the input file for gubbins. $\endgroup$ Feb 17, 2019 at 15:18
-
1$\begingroup$ gubbins seems to be a program for the fiel of phylogenetics. So the input file must be an alignment file between your species and not the alignments of the reads for one sample. $\endgroup$ Feb 17, 2019 at 15:54
-
$\begingroup$ @AhmedAbdullah you should probably ask a new question, explaining your final objective. Going from reads aligned to a genome (bam) to a multiple sequence alignment you want to use for phylogenetic analysis makes no sense. If you explain what you are trying to achieve in a new question, we should be able to help. $\endgroup$– terdon ♦Feb 17, 2019 at 19:10
-
$\begingroup$ @terdon I am trying whether this works or not. I have 48 SAM files corresponding to 48 samples. I converted them to 48 ordered BAM file. and the BAM fiels to a single vcf file using : bcftools mpileup -d 1000 -Bf R.fa *.bam | bcftools call --skip-variants indels --multiallelic-caller --variants-only --samples-file samples.txt -O z -o aln.vcf.gz. vcf-kit can generate aligned fasta for only SNPs from this vcf file (vk phylo fast <vcf fiie>. The problem is alignment is only of SNPs, I was wondering whether one can generate full alignment from the (SAM files or the subsequent vcf file) $\endgroup$ Feb 18, 2019 at 9:10
-
$\begingroup$ @AhmedAbdullah you are almost certainly approaching the problem in the wrong way. Please ask a new question explaining exactly what you need. What you're describing here makes no sense. $\endgroup$– terdon ♦Feb 18, 2019 at 9:13
$\begingroup$
$\endgroup$
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: genomicalignments::stackStringsFromBam().
$\begingroup$
$\endgroup$
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.
