How to obtain a full length sequence from bam file

[kamilo889]

I have an BAM file from a single chromosome, and I am interested in obtaining a piece of this sequence in FASTA format as consensus.

I've already try to obtain the sequence and I am interested in put it in another separated file, but this is not a continuous sequence, so I am looking for a way to obtain this piece of sequence in FASTA format. I have tried mpileup, but it does not work (it gives me empty files).

Does anyone knows how it could work?

Thank you.

[deyler]

It's not clear to me what you mean by consensus. Samtools doesn't generate a consensus sequence; it calls variants relative to the reference sequence. So it seems to me you have two possible workflows:

1) use samtools to call variants in the region of interest (ROI)
1a) extract the reference sequence for your ROI from the reference fasta file with samtools faidx
1b) select the variants you want reflected in your consensus sequence based on whatever criteria you have in mind
1c) change the sequence of your ROI to match the variants from 1b

or

2) use samtools view to extract just the reads from your ROI
2a) trim the ends of reads that extend outside the ends of your ROI with a custom script
2b) output the reads in some multiple-sequence alignment format like CLUSTAL
2c) take the multiple-sequence alignment file and put it into some program that generates consensus sequences

CAVEATS:
1) I am not sure I really understand what you are trying to do
2) These instructions pertain to samtools v 0.19+, not the most recent version

[GenoMax]

It is possible to get a consensus fastq sequence using samtools/bcftools.

http://samtools.sourceforge.net/mpileup.shtml
https://www.biostars.org/p/1388/
https://www.biostars.org/p/70585/

[GenoMax]

@kamilo889: What version of samtools did you use to generate the original BAM? Are you using the same version to try and generate the consensus (you likely can't mix and match old/new versions). Unless you have coverage across the entire chromosome I am not sure you will get "full length" sequence.

[deyler]

What I dislike about the samtools-oriented approach is that it somewhat obscures what is really going on. If samtools has no data at a position, it assumes it's reference (which may not be valid). I would prefer a consensus sequence to reflect the absence of information, rather than to show me a consensus base call based on no data. Furthermore it seems like variant calls are not filtered when you pipe them straight into vcfutils, so you might be feeding in garbage. #1 from my post above gives you the opportunity to filter variant calls; option #2 lets you choose a consensus sequence algorithm other than samtools.

It all depends on how good of a consensus sequence you need, I suppose, and what you're using it for. I wouldn't use #1 (mpileup | bcftools | vcfutils) to describe the consensus sequence for a binding site. But for generating an improved reference sequence for alignment, it may be fine.