Generating consensus sequence from bam file

Question

[elienesantossilva09]

Hello, I am trying to generate consensus bam file sequences. For this I used the commands recommended by the link below.

https://samtools.github.io/bcftools/howtos/consensus-sequence.html`][

I did a test and managed to generate a consensus in fasta format. but when I used the command below (first command) I received a message that my sample was considered diploid because I did not specify the ploidy. in my case it is mitochondrial dna and i need to specify in first command. How do I do that?

1 command

bcftools mpileup -Ou -f NC_012920.1.fasta 1-321.bam | bcftools call -mv -Oz -o 1-321.vcf.gz
bcftools index 1-321.vcf.gz

other sequential commands:

     (base)liene@gar:~/Desktop/doc_liene/vcfs$ 
     (base)bcftools norm -f NC_012920.1.fasta 1-321.vcf.gz -o 1-321N.bcf
     Lines   total/split/realigned/skipped: 38/0/3/0
     (base) liene@gar:~/Desktop/doc_liene/vcfs$ 
     bcftools filter --IndelGap 5 1-321N.bcf -Ob -o 1-321NI.bcf
     (base) liene@gar:~/Desktop/doc_liene/vcfs$ 
     cat NC_012920.1.fasta |bcftools consensus 1-321.vcf.gz > 1-321MTconsensus.fa
      Note: the --sample option not given, applying all records
      The site NC_012920.1:8270 overlaps with another variant, skipping...
      The site NC_012920.1:8272 overlaps with another variant, skipping...
      (base) liene@gar:~/Desktop/doc_liene/vcfs

Another question I have is whether I made all the other commands right. although it generated bcf files, they were not used in the consensus extract command syntax. I would like to understand that. Thanks

[elienesantossilva09]

Thank you very much. It was of great help.

[ATpoint]

Please use ADD COMMENT for comments. If the answer was helpful, consider to upvote and accept:

[GenoMax]

Did you add this to the original question today?

Another question I have is whether I made all the other commands right. although it generated bcf files, they were not used in the consensus extract command syntax. I would like to understand that.

Adding new questions to original post is not going to get them noticed (we don't have a good way to keep track unless you explicitly add a manual time stamp).

Answer 1

[finswimmer]

Hello,

the ploidy can be set explicit for a region and sex using a ploidy file. That is a tab delimited file with the columns CHROM, FROM, TO, SEX, PLOIDY. If no ploidy file is specified with --ploidy-file bcftools assumes this default:

X 1 60000 M 1
X 2699521 154931043 M 1
Y 1 59373566 M 1
Y 1 59373566 F 0
MT 1 16569 M 1
MT 1 16569 F 1
*  * *     M 2
*  * *     F 2

So if your mitochondrial reference sequence is name MT you don't need to do anything. There are more defaults you can set with --ploidy. See the output of bcftools call --ploidy list for more information.

although it generated bcf files, they were not used in the consensus extract command syntax. I would like to understand that.

You normalize and filter your vcf file after the calling step, but you take the original vcf file as an input for the consensus creation. I guess this is not what you want to do. Just take the filtered bcf file if this is what you want to do.

fin swimmer