Convert Bam To Wig
1
3
Entering edit mode
11.2 years ago
Assa Yeroslaviz ★ 1.7k

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

bam conversion wiggle samtools • 12k views
ADD COMMENT
0
Entering edit mode

what's the difference between your two commands ?

ADD REPLY
0
Entering edit mode

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

ADD REPLY
7
Entering edit mode
11.2 years ago

Try to use

samtools view -f 0x10 -b

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

samtools view -h -f 0x10
ADD COMMENT
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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 :)

ADD REPLY
0
Entering edit mode

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

ADD REPLY
0
Entering edit mode

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

ADD REPLY

Login before adding your answer.

Traffic: 784 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6