Question: Trying to get consensus sequence from paired end alignment
0
gravatar for Sus
12 months ago by
Sus10
Sus10 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 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.
ADD COMMENTlink written 12 months ago by Sus10

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 12 months ago • written 12 months ago by finswimmer11k
Please log in to add an answer.

Help
Access

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