Trying to get consensus sequence from paired end alignment

Question

[Sus]

Hello !

I am trying to get a consensus sequence from my BAM file. For the context, I'm working with paired-end reads in fastq format, so far I've done:

• Building the BWA index for my reference genome
• Alignin with the BWA mem following this command: bwa mem -t 16 ref.fa R1_1.fq R1_2.fq > aln.sam
• Converting SAM to BAM with SAMTOOLS
• Sorting my BAM file and indexing it

I can successfully open my BAM file on IGV to extract the consensus sequence but I want to do include the IUPAC notation for ambiguities.

My issue comes when I try to extract the consensus sequence, the following command I'm using is giving me an empty file at the end and I don't understand why: samtools mpileup -uf ref.fasta aln_sorted.bam | bcftools view -cg - | perl vcfutils.pl vcf2fq > test.fa

Here's the output of the command line above:

[mpileup] 1 samples in 1 input files
<mpileup> Set max per-file depth to 8000
[afs] 0:0.000 1:0.000 2:0.000
Use of uninitialized value $l in numeric lt (<) at vcfutils.pl line 566, <> line 36.
Use of uninitialized value $l in numeric lt (<) at vcfutils.pl line 566, <> line 36.

[finswimmer]

Hello Sus,

check first if the output of bcftools is what you expected before trying to pipe it to vcfutils. Instead of using vcfutils you could also give bcftools consensus a try.

fin swimmer