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.