Error while trying to get consensus Fastq from BAM

Question

[shreyasibiswas88]

Hi all,

I am trying to get consensus Fastq from BAM file and I am running into this error. I looked for an answer in many posts but could not resolve it. Please help.

command:

/home/samtools-1.2/samtools mpileup -uf Rattus_norvegicus.Rnor_5.0.dna_sm.fa dedupRratUK_recal.bam -d 8000 | /home/bcftools-1.2/bcftools view -cg - | /home/bcftools-1.2/vcfutils.pl vcf2fq -d 10 -D 8000 > RratUK.fastq

Error:

Error: Could not parse --min-ac g
[mpileup] 1 samples in 1 input files
Use of uninitialized value in length at /home/sb47/Programs/bcftools-1.2/vcfutils.pl line 565.
Use of uninitialized value in length at /home/sb47/Programs/bcftools-1.2/vcfutils.pl line 565.

Thank you.

[stolarek.ir]

try running step by step (make sure bam is sorted. Or try vcfutils that is attached to samtools /usr/share/samtools/vcfutils.pl)

Answer 1

[Jorge Amigo]

the commands you're using (pileup generation + variant calling + fastq generation) seem correct, although the version of bcftools you're using doesn't work with the options you've probably seen on many old posts. the problem is that newer bcftools uses call to call variants, rather than the former bcftools view -cg

/home/samtools-1.2/samtools mpileup -uf Rattus_norvegicus.Rnor_5.0.dna_sm.fa dedupRratUK_recal.bam -d 8000 \
| /home/bcftools-1.2/bcftools call -mv \
| /home/bcftools-1.2/vcfutils.pl vcf2fq -d 10 -D 8000 \
> RratUK.fastq

[Lina F]

Just a note, the vcfutils.pl script threw an error for me. However, I was able to use this command to generate fasta consensus sequences:

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

Source: https://github.com/samtools/bcftools/wiki/HOWTOs#consensus-calling

I am using samtools v1.3.1, htslib v1.3.1, and bcftools v1.3.1

[cdastmalchi]

Hi there, when I try running these commands, the line cat e1a.fasta | bcftools consensus calls.vcf.gz > EG28-SEQ-2020129-79.cns.fa gives me a Segmentation fault (core dumped). Any idea how that can be addressed?

Answer 2

[Jorge Amigo]

fastq generation from a bam file is much easier: the reads are just there, and you just have to extract them. several tools may help you with this, but my personal choice is bedtools bamtofastq

bedtools bamtofastq [OPTIONS] -i <BAM> -fq <FASTQ>

although the newer samtools bam2fq option seems to perform fine too.

edit: this answer only applies if all the reads from the bam file are to be extracted in fastq format, but not to generate the consensus fastq requested.

[Istvan Albert]

This is probably not what the OP is asking, the title is misleading. They probably want to get the consensus FASTA file out of the BAM file pileup.

[Jorge Amigo]

you're probably right Istvan. I think it was me that I read it too fast ;)