This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Trying to get consensus sequence from paired end alignment

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.
paired-end alignment bwa samtools
• 2,656 views • link
updated by Biostar Biostar • written by Sus • Sus

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

• 1 views • link
written by finswimmer finswimmer

0 answers

No answers yet.

Log in to answer this question.

More posts like this