Question: Trying to get consensus sequence from paired end alignment
0
gravatar for Sus
18 months ago by
Sus20
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 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 18 months 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 18 months ago • written 18 months ago by finswimmer12k
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: 950 users visited in the last hour