Can anyone suggest the best way to convert bam to fasta in an aligned format?
I have a data from ddRad sequencing. Initially, I want to count the number of restriction sites obtained. So, I already got suggestions to get it mapped to the genome and count the sites using available python code.
Following are the commands I have used for fasta conversion.
samtools mpileup -uf ref.fa aln.bam | bcftools call -mv -Oz -o calls.vcf.gz tabix calls.vcf.gz cat ref.fa | bcftools consensus calls.vcf.gz > cns.fa
Could you please have look at this and please tell me know if this is the right format where I could get a reliable number of restriction sites?
Other pipeline suggestions are also most welcome.