How to extract a regional consensus sequence of reads from BAM file

[jeffseq]

Deal all,

How do I get a regional consensus sequence, say ch01: 77295611-77295961, from a BAM file by Samtools?

I extracted bam file for this region: ch01: 77295611-77295961 (region.bam) from all_aln.bam and used the recommanded command:

samtools mpileup -uf ref.fa region.bam | bcftools view -cg - | vcfutils.pl vcf2fq > cns.fq

However, the cns.fq seems contain sequence for the whole reference (took longer time), not just the region ch01: 77295611-77295961.

How to fix this?

Thanks.

Jeff

[jeffseq]

anyone knows or do I need to write a script to parse it out?

[swbarnes2]

Looking at the command line options for mpileup, why didn't you use -r there to specify a region?

Usage: samtools mpileup [options] in1.bam [in2.bam [...]]

Input options:

-6 assume the quality is in the Illumina-1.3+ encoding
-A count anomalous read pairs
-B disable BAQ computation
-b FILE list of input BAM files [null]
-C INT parameter for adjusting mapQ; 0 to disable [0]
-d INT max per-BAM depth to avoid excessive memory usage [250]
-E extended BAQ for higher sensitivity but lower specificity
-f FILE faidx indexed reference sequence file [null]
-G FILE exclude read groups listed in FILE [null]
-l FILE list of positions (chr pos) or regions (BED) [null]
-M INT cap mapping quality at INT [60]
-r STR region in which pileup is generated [null]
-R ignore RG tags
-q INT skip alignments with mapQ smaller than INT [0]
-Q INT skip bases with baseQ/BAQ smaller than INT [13]

[jeffseq]

Thanks. Tried:

samtools mpileup -uf ref.fa region.bam -r ch01: 77295611-77295961 | bcftools view -cg - | vcfutils.pl vcf2fq > cns.fq

the consensus is indeed calculated for region 77295611-77295961 on ch01, however, the cns.fq contains sequence sprend over all ch01, in the form of NNNNN-consensus-NNNNN format. (Need to parse out by writing a script).

Jeff

[kjlee]

I believe that your problem is the space between the "ch01:" and "77295611-77295961". Perhaps try it without the space.

[mmmm]

which command is used to get the complete consensus from the bam file?

[adamr]

Hi all,

I tried this and ran into the same problem. When I broke up the pipe, it looked like bcftools was giving the expected output, but vcfutils then tried to recreate the entire chromosome.

My solution was to make a small modification to vcfutils so that it does not add a "gap" at the beginning of the sequence.

Line 500 now reads:
if (($t[1] - $last_pos > 1) && ($last_pos > 0)) {

That seems to do it...

[nazen]

to adamr:
Thanks a lot it works!

[bwubb]

Ive been trying to work with this pipe in order to generate short sequences.

samtools mpileup -uf ref.fa region.bam | bcftools view -cg - | vcfutils.pl vcf2fq > cns.fq

Ive even included adamr's changes to vcfutils.pl. The above code works but only for a single region. If I try to give a bed file there are issues. The sequences are not named based on the BED format column.

Sample Bed:

1	1158620	1158631	"SDF4.1:1158631"
1	1650793	1650804	"CDK11A,CDK11B;CDK11B.1:1650797"

Also if just does not work for the multiple regions. I will get something like

@1
CATTTTGCnnnnnnnnnnnnnnnnnnnnnnnnnnnn....
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!....
~~~~~~~~~~~~~~

At the very least I was hoping for a multiple sequence fq file. Does anyone have any insight?