It's probably the same question over and over but I can't find the solutions.
I mapped some contigs on a Plant reference genome with BWA, sorted the file with samtools and got the coverage by typing:
samtools depth my.bam > qry-depth> wc -l qry-depth
The results of this command / sequence legth * 100 to have the % genome covered.
Then, I wanted to extract the consensus sequence from the bam file:
samtools mpileup -uf reference.fasta my.bam | bcftools call -mv -Oz -o calls.vcf.gz>
tabix calls.vcf.gz> cat reference.fasta | bcftools consensus calls.vcf.gz > con.fasta
And here I am. Now I have a problem and a general question.
For two files I got this error message after the first mpileup commad:
The QS annotation not present at gi|33590464|gb|AY244516.1|:2123 (NCBI header)
So, I can't get the consensus.
For other samples, I know the contigs don't cover the entire genome, so how does sometools export the consensus? By keeping in the not covered region the reference sequence? Is there a way to visualize this?
I hope I'm clear and not repetitive,
Thanks for your help :)