Convert Bam To Wig

Question

[Assa Yeroslaviz]

Hi everybody,

I have found an old description of how to convert a bam file into wig using perl.

(How can I convert BAM/SAM to wiggle).

In this example (as far as I understood it) the whole bam file was converted.

samtools pileup bam_small.sorted.bam  |\
perl -ne `'BEGIN{print "track type=wiggle_0 name=fileName description=fileName\n"};($c, $start, undef, $depth) = split; if ($c ne $lastC) { print "variableStep chrom=chr$c\n"; };$lastC=$c;next unless $. % 10 ==0;print "$start\t$depth\n" unless $depth<3;' \
 > fileName1.wig

This command works just fine for single or multiple chromosomes.

Now I wanted to change the command to make a strand-specific wig file for forward and reverse strands, so I tried to tweak the command a bit. AS pileup is deprecated, I used mpileup:

samtools view -h -f 0x10 bam_small.sorted.bam |\
samtools mpileup -S -  |\
perl -ne 'BEGIN{print "track type=wiggle_0 name=fileName description=fileName\n"};($c, $start, undef, $depth) = split; if ($c ne $lastC) { print "variableStep chrom=chr$c\n"; };$lastC=$c;next unless $. % 10 ==0;print "$start\t$depth\n" unless $depth<3;' \
> fileName1.wig

The 0x10 flag is used for the strand of the query, but somehow it doesn't work. I keep getting this error massage:

[bam_header_read] invalid BAM binary header (this is not a BAM file).

Does anyone has an idea how to run this command?

Thanks

Assa

[Pierre Lindenbaum]

what's the difference between your two commands ?

[Assa Yeroslaviz]

sorry I posted the same command twice. Now I corrected it.

Answer 1

[Pierre Lindenbaum]

Try to use

samtools view -f 0x10 -b

(-b is output BAM, -h forces output to SAM), instead of

samtools view -h -f 0x10

[Assa Yeroslaviz]

nice, thank you very much. I do have a followup question.

Why is the number of the inputs for each strand doesn't equal the total number. When I'm using the -f/-F filter option I get almost the same number of input lines as without these parameters. wc -l: -f: 113905 -F: 121720 w.o.: 141199

[John]

It's because reads can have the reverse strand flag set to 1, but not be mapped (in some aligners - most notably BWA).

Use -f 0x10 -F 0x4 to be sure :)

[Pierre Lindenbaum]

I guess it's because some reads are not mapped at all.

[Assa Yeroslaviz]

so do thy come in both strands? Aren't they suppose to be in the -f option only?