Question: Trying to get consensus sequence from paired end alignment
gravatar for Sus
2.1 years ago by
Sus20 wrote:

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 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 line 566, <> line 36.
Use of uninitialized value $l in numeric lt (<) at line 566, <> line 36.
bwa samtools alignment paired-end • 1.1k views
ADD COMMENTlink written 2.1 years ago by Sus20

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

ADD REPLYlink modified 2.1 years ago • written 2.1 years ago by finswimmer13k
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 626 users visited in the last hour